source/unres/src_MD/xdrf_em64/libxdrf.m4

   1 /*____________________________________________________________________________
   2  |
   3  | libxdrf - portable fortran interface to xdr. some xdr routines
   4  |           are C routines for compressed coordinates
   5  |
   6  | version 1.1
   7  |
   8  | This collection of routines is intended to write and read
   9  | data in a portable way to a file, so data written on one type
  10  | of machine can be read back on a different type.
  11  |
  12  | all fortran routines use an integer 'xdrid', which is an id to the
  13  | current xdr file, and is set by xdrfopen.
  14  | most routines have in integer 'ret' which is the return value.
  15  | The value of 'ret' is zero on failure, and most of the time one
  16  | on succes.
  17  |
  18  | There are three routines useful for C users:
  19  |  xdropen(), xdrclose(), xdr3dfcoord().
  20  | The first two replace xdrstdio_create and xdr_destroy, and *must* be
  21  | used when you plan to use xdr3dfcoord(). (they are also a bit
  22  | easier to interface). For writing data other than compressed coordinates
  23  | you should use the standard C xdr routines (see xdr man page)
  24  |
  25  | xdrfopen(xdrid, filename, mode, ret)
  26  |      character *(*) filename
  27  |      character *(*) mode
  28  |
  29  |      this will open the file with the given filename (string)
  30  |      and the given mode, it returns an id in xdrid, which is
  31  |      to be used in all other calls to xdrf routines.
  32  |      mode is 'w' to create, or update an file, for all other
  33  |      values of mode the file is opened for reading
  34  |
  35  |      you need to call xdrfclose to flush the output and close
  36  |      the file.
  37  |      Note that you should not use xdrstdio_create, which comes with the
  38  |      standard xdr library
  39  |
  40  | xdrfclose(xdrid, ret)
  41  |      flush the data to the file, and closes the file;
  42  |      You should not use xdr_destroy (which comes standard with
  43  |      the xdr libraries.
  44  |
  45  | xdrfbool(xdrid, bp, ret)
  46  |      integer pb
  47  |
  48  |      This filter produces values of either 1 or 0
  49  |
  50  | xdrfchar(xdrid, cp, ret)
  51  |      character cp
  52  |
  53  |      filter that translate between characters and their xdr representation
  54  |      Note that the characters in not compressed and occupies 4 bytes.
  55  |
  56  | xdrfdouble(xdrid, dp, ret)
  57  |      double dp
  58  |
  59  |      read/write a double.
  60  |
  61  | xdrffloat(xdrid, fp, ret)
  62  |      float fp
  63  |
  64  |      read/write a float.
  65  |
  66  | xdrfint(xdrid, ip, ret)
  67  |      integer ip
  68  |
  69  |      read/write integer.
  70  |
  71  | xdrflong(xdrid, lp, ret)
  72  |      integer lp
  73  |
  74  |      this routine has a possible portablility problem due to 64 bits longs.
  75  |
  76  | xdrfshort(xdrid, sp, ret)
  77  |      integer *2 sp
  78  |
  79  | xdrfstring(xdrid, sp, maxsize, ret)
  80  |      character *(*)
  81  |      integer maxsize
  82  |
  83  |      read/write a string, with maximum length given by maxsize
  84  |
  85  | xdrfwrapstring(xdris, sp, ret)
  86  |      character *(*)
  87  |
  88  |      read/write a string (it is the same as xdrfstring accept that it finds
  89  |      the stringlength itself.
  90  |
  91  | xdrfvector(xdrid, cp, size, xdrfproc, ret)
  92  |      character *(*)
  93  |      integer size
  94  |      external xdrfproc
  95  |
  96  |      read/write an array pointed to by cp, with number of elements
  97  |      defined by 'size'. the routine 'xdrfproc' is the name
  98  |      of one of the above routines to read/write data (like xdrfdouble)
  99  |      In contrast with the c-version you don't need to specify the
 100  |      byte size of an element.
 101  |      xdrfstring is not allowed here (it is in the c version)
 102  |
 103  | xdrf3dfcoord(xdrid, fp, size, precision, ret)
 104  |      real (*) fp
 105  |      real precision
 106  |      integer size
 107  |
 108  |      this is *NOT* a standard xdr routine. I named it this way, because
 109  |      it invites people to use the other xdr routines.
 110  |      It is introduced to store specifically 3d coordinates of molecules
 111  |      (as found in molecular dynamics) and it writes it in a compressed way.
 112  |      It starts by multiplying all numbers by precision and
 113  |      rounding the result to integer. effectively converting
 114  |      all floating point numbers to fixed point.
 115  |      it uses an algorithm for compression that is optimized for
 116  |      molecular data, but could be used for other 3d coordinates
 117  |      as well. There is subtantial overhead involved, so call this
 118  |      routine only if you have a large number of coordinates to read/write
 119  |
 120  | ________________________________________________________________________
 121  |
 122  | Below are the routines to be used by C programmers. Use the 'normal'
 123  | xdr routines to write integers, floats, etc (see man xdr)
 124  |
 125  | int xdropen(XDR *xdrs, const char *filename, const char *type)
 126  |      This will open the file with the given filename and the
 127  |      given mode. You should pass it an allocated XDR struct
 128  |      in xdrs, to be used in all other calls to xdr routines.
 129  |      Mode is 'w' to create, or update an file, and for all
 130  |      other values of mode the file is opened for reading.
 131  |      You need to call xdrclose to flush the output and close
 132  |      the file.
 133  |
 134  |      Note that you should not use xdrstdio_create, which
 135  |      comes with the standard xdr library.
 136  |
 137  | int xdrclose(XDR *xdrs)
 138  |      Flush the data to the file, and close the file;
 139  |      You should not use xdr_destroy (which comes standard
 140  |      with the xdr libraries).
 141  |
 142  | int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
 143  |      This is \fInot\fR a standard xdr routine. I named it this
 144  |      way, because it invites people to use the other xdr
 145  |      routines.
 146  |
 147  |      (c) 1995 Frans van Hoesel, hoesel@chem.rug.nl
 148 */
 149
 150
 151 #include <limits.h>
 152 #include <malloc.h>
 153 #include <math.h>
 154 /* #include <rpc/rpc.h>
 155 #include <rpc/xdr.h> */
 156 #include "xdr.h"
 157 #include <stdio.h>
 158 #include <stdlib.h>
 159 #include "xdrf.h"
 160
 161 int ftocstr(char *, int, char *, int);
 162 int ctofstr(char *, int, char *);
 163
 164 #define MAXID 20
 165 static FILE *xdrfiles[MAXID];
 166 static XDR *xdridptr[MAXID];
 167 static char xdrmodes[MAXID];
 168 static unsigned int cnt;
 169
 170 typedef void (* FUNCTION(xdrfproc)) (int *, void *, int *);
 171
 172 void
 173 FUNCTION(xdrfbool) ARGS(`xdrid, pb, ret')
 174 int *xdrid, *ret;
 175 int *pb;
 176 {
 177         *ret = xdr_bool(xdridptr[*xdrid], (bool_t *) pb);
 178         cnt += sizeof(int);
 179 }
 180
 181 void
 182 FUNCTION(xdrfchar) ARGS(`xdrid, cp, ret')
 183 int *xdrid, *ret;
 184 char *cp;
 185 {
 186         *ret = xdr_char(xdridptr[*xdrid], cp);
 187         cnt += sizeof(char);
 188 }
 189
 190 void
 191 FUNCTION(xdrfdouble) ARGS(`xdrid, dp, ret')
 192 int *xdrid, *ret;
 193 double *dp;
 194 {
 195         *ret = xdr_double(xdridptr[*xdrid], dp);
 196         cnt += sizeof(double);
 197 }
 198
 199 void
 200 FUNCTION(xdrffloat) ARGS(`xdrid, fp, ret')
 201 int *xdrid, *ret;
 202 float *fp;
 203 {
 204         *ret = xdr_float(xdridptr[*xdrid], fp);
 205         cnt += sizeof(float);
 206 }
 207
 208 void
 209 FUNCTION(xdrfint) ARGS(`xdrid, ip, ret')
 210 int *xdrid, *ret;
 211 int *ip;
 212 {
 213         *ret = xdr_int(xdridptr[*xdrid], ip);
 214         cnt += sizeof(int);
 215 }
 216
 217 void
 218 FUNCTION(xdrflong) ARGS(`xdrid, lp, ret')
 219 int *xdrid, *ret;
 220 long *lp;
 221 {
 222         *ret = xdr_long(xdridptr[*xdrid], lp);
 223         cnt += sizeof(long);
 224 }
 225
 226 void
 227 FUNCTION(xdrfshort) ARGS(`xdrid, sp, ret')
 228 int *xdrid, *ret;
 229 short *sp;
 230 {
 231         *ret = xdr_short(xdridptr[*xdrid], sp);
 232         cnt += sizeof(sp);
 233 }
 234
 235 void
 236 FUNCTION(xdrfuchar) ARGS(`xdrid, ucp, ret')
 237 int *xdrid, *ret;
 238 char *ucp;
 239 {
 240         *ret = xdr_u_char(xdridptr[*xdrid], ucp);
 241         cnt += sizeof(char);
 242 }
 243
 244 void
 245 FUNCTION(xdrfulong) ARGS(`xdrid, ulp, ret')
 246 int *xdrid, *ret;
 247 unsigned long *ulp;
 248 {
 249         *ret = xdr_u_long(xdridptr[*xdrid], ulp);
 250         cnt += sizeof(unsigned long);
 251 }
 252
 253 void
 254 FUNCTION(xdrfushort) ARGS(`xdrid, usp, ret')
 255 int *xdrid, *ret;
 256 unsigned short *usp;
 257 {
 258         *ret = xdr_u_short(xdridptr[*xdrid], usp);
 259         cnt += sizeof(unsigned short);
 260 }
 261
 262 void
 263 FUNCTION(xdrf3dfcoord) ARGS(`xdrid, fp, size, precision, ret')
 264 int *xdrid, *ret;
 265 float *fp;
 266 int *size;
 267 float *precision;
 268 {
 269         *ret = xdr3dfcoord(xdridptr[*xdrid], fp, size, precision);
 270 }
 271
 272 void
 273 FUNCTION(xdrfstring) ARGS(`xdrid, STRING_ARG(sp), maxsize, ret')
 274 int *xdrid, *ret;
 275 STRING_ARG_DECL(sp);
 276 int *maxsize;
 277 {
 278         char *tsp;
 279
 280         tsp = (char*) malloc(((STRING_LEN(sp)) + 1) * sizeof(char));
 281         if (tsp == NULL) {
 282             *ret = -1;
 283             return;
 284         }
 285         if (ftocstr(tsp, *maxsize+1, STRING_PTR(sp), STRING_LEN(sp))) {
 286             *ret = -1;
 287             free(tsp);
 288             return;
 289         }
 290         *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int) *maxsize);
 291         ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
 292         cnt += *maxsize;
 293         free(tsp);
 294 }
 295
 296 void
 297 FUNCTION(xdrfwrapstring) ARGS(`xdrid,  STRING_ARG(sp), ret')
 298 int *xdrid, *ret;
 299 STRING_ARG_DECL(sp);
 300 {
 301         char *tsp;
 302         int maxsize;
 303         maxsize = (STRING_LEN(sp)) + 1;
 304         tsp = (char*) malloc(maxsize * sizeof(char));
 305         if (tsp == NULL) {
 306             *ret = -1;
 307             return;
 308         }
 309         if (ftocstr(tsp, maxsize, STRING_PTR(sp), STRING_LEN(sp))) {
 310             *ret = -1;
 311             free(tsp);
 312             return;
 313         }
 314         *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int)maxsize);
 315         ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
 316         cnt += maxsize;
 317         free(tsp);
 318 }
 319
 320 void
 321 FUNCTION(xdrfopaque) ARGS(`xdrid, cp, ccnt, ret')
 322 int *xdrid, *ret;
 323 caddr_t *cp;
 324 int *ccnt;
 325 {
 326         *ret = xdr_opaque(xdridptr[*xdrid], (caddr_t)*cp, (u_int)*ccnt);
 327         cnt += *ccnt;
 328 }
 329
 330 void
 331 FUNCTION(xdrfsetpos) ARGS(`xdrid, pos, ret')
 332 int *xdrid, *ret;
 333 int *pos;
 334 {
 335         *ret = xdr_setpos(xdridptr[*xdrid], (u_int) *pos);
 336 }
 337
 338 void
 339 FUNCTION(xdrf) ARGS(`xdrid, pos')
 340 int *xdrid, *pos;
 341 {
 342         *pos = xdr_getpos(xdridptr[*xdrid]);
 343 }
 344
 345 void
 346 FUNCTION(xdrfvector) ARGS(`xdrid, cp, size, elproc, ret')
 347 int *xdrid, *ret;
 348 char *cp;
 349 int *size;
 350 FUNCTION(xdrfproc) elproc;
 351 {
 352         int lcnt;
 353         cnt = 0;
 354         for (lcnt = 0; lcnt < *size; lcnt++) {
 355                 elproc(xdrid, (cp+cnt) , ret);
 356         }
 357 }
 358
 359
 360 void
 361 FUNCTION(xdrfclose) ARGS(`xdrid, ret')
 362 int *xdrid;
 363 int *ret;
 364 {
 365         *ret = xdrclose(xdridptr[*xdrid]);
 366         cnt = 0;
 367 }
 368
 369 void
 370 FUNCTION(xdrfopen) ARGS(`xdrid,  STRING_ARG(fp), STRING_ARG(mode), ret')
 371 int *xdrid;
 372 STRING_ARG_DECL(fp);
 373 STRING_ARG_DECL(mode);
 374 int *ret;
 375 {
 376         char fname[512];
 377         char fmode[3];
 378
 379         if (ftocstr(fname, sizeof(fname), STRING_PTR(fp), STRING_LEN(fp))) {
 380                 *ret = 0;
 381         }
 382         if (ftocstr(fmode, sizeof(fmode), STRING_PTR(mode),
 383                         STRING_LEN(mode))) {
 384                 *ret = 0;
 385         }
 386
 387         *xdrid = xdropen(NULL, fname, fmode);
 388         if (*xdrid == 0)
 389                 *ret = 0;
 390         else
 391                 *ret = 1;
 392 }
 393
 394 /*___________________________________________________________________________
 395  |
 396  | what follows are the C routines for opening, closing xdr streams
 397  | and the routine to read/write compressed coordinates together
 398  | with some routines to assist in this task (those are marked
 399  | static and cannot be called from user programs)
 400 */
 401 #define MAXABS INT_MAX-2
 402
 403 #ifndef MIN
 404 #define MIN(x,y) ((x) < (y) ? (x):(y))
 405 #endif
 406 #ifndef MAX
 407 #define MAX(x,y) ((x) > (y) ? (x):(y))
 408 #endif
 409 #ifndef SQR
 410 #define SQR(x) ((x)*(x))
 411 #endif
 412 static int magicints[] = {
 413     0, 0, 0, 0, 0, 0, 0, 0, 0,
 414     8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
 415     80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
 416     812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
 417     8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
 418     82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
 419     832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
 420     8388607, 10568983, 13316085, 16777216 };
 421
 422 #define FIRSTIDX 9
 423 /* note that magicints[FIRSTIDX-1] == 0 */
 424 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
 425
 426
 427 /*__________________________________________________________________________
 428  |
 429  | xdropen - open xdr file
 430  |
 431  | This versions differs from xdrstdio_create, because I need to know
 432  | the state of the file (read or write) so I can use xdr3dfcoord
 433  | in eigther read or write mode, and the file descriptor
 434  | so I can close the file (something xdr_destroy doesn't do).
 435  |
 436 */
 437
 438 int xdropen(XDR *xdrs, const char *filename, const char *type) {
 439     static int init_done = 0;
 440     enum xdr_op lmode;
 441     const char *type1;
 442     int xdrid;
 443
 444     if (init_done == 0) {
 445         for (xdrid = 1; xdrid < MAXID; xdrid++) {
 446             xdridptr[xdrid] = NULL;
 447         }
 448         init_done = 1;
 449     }
 450     xdrid = 1;
 451     while (xdrid < MAXID && xdridptr[xdrid] != NULL) {
 452         xdrid++;
 453     }
 454     if (xdrid == MAXID) {
 455         return 0;
 456     }
 457     if (*type == 'w' || *type == 'W') {
 458             type = "w+";
 459             type1 = "w+";
 460             lmode = XDR_ENCODE;
 461     } else if (*type == 'a' || *type == 'A') {
 462             type = "w+";
 463             type1 = "a+";
 464             lmode = XDR_ENCODE;
 465     } else {
 466             type = "r";
 467             type1 = "r";
 468             lmode = XDR_DECODE;
 469     }
 470     xdrfiles[xdrid] = fopen(filename, type1);
 471     if (xdrfiles[xdrid] == NULL) {
 472         xdrs = NULL;
 473         return 0;
 474     }
 475     xdrmodes[xdrid] = *type;
 476     /* next test isn't usefull in the case of C language
 477      * but is used for the Fortran interface
 478      * (C users are expected to pass the address of an already allocated
 479      * XDR staructure)
 480      */
 481     if (xdrs == NULL) {
 482         xdridptr[xdrid] = (XDR *) malloc(sizeof(XDR));
 483         xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode);
 484     } else {
 485         xdridptr[xdrid] = xdrs;
 486         xdrstdio_create(xdrs, xdrfiles[xdrid], lmode);
 487     }
 488     return xdrid;
 489 }
 490
 491 /*_________________________________________________________________________
 492  |
 493  | xdrclose - close a xdr file
 494  |
 495  | This will flush the xdr buffers, and destroy the xdr stream.
 496  | It also closes the associated file descriptor (this is *not*
 497  | done by xdr_destroy).
 498  |
 499 */
 500
 501 int xdrclose(XDR *xdrs) {
 502     int xdrid;
 503
 504     if (xdrs == NULL) {
 505         fprintf(stderr, "xdrclose: passed a NULL pointer\n");
 506         exit(1);
 507     }
 508     for (xdrid = 1; xdrid < MAXID; xdrid++) {
 509         if (xdridptr[xdrid] == xdrs) {
 510
 511             xdr_destroy(xdrs);
 512             fclose(xdrfiles[xdrid]);
 513             xdridptr[xdrid] = NULL;
 514             return 1;
 515         }
 516     }
 517     fprintf(stderr, "xdrclose: no such open xdr file\n");
 518     exit(1);
 519
 520 }
 521
 522 /*____________________________________________________________________________
 523  |
 524  | sendbits - encode num into buf using the specified number of bits
 525  |
 526  | This routines appends the value of num to the bits already present in
 527  | the array buf. You need to give it the number of bits to use and you
 528  | better make sure that this number of bits is enough to hold the value
 529  | Also num must be positive.
 530  |
 531 */
 532
 533 static void sendbits(int buf[], int num_of_bits, int num) {
 534
 535     unsigned int cnt, lastbyte;
 536     int lastbits;
 537     unsigned char * cbuf;
 538
 539     cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
 540     cnt = (unsigned int) buf[0];
 541     lastbits = buf[1];
 542     lastbyte =(unsigned int) buf[2];
 543     while (num_of_bits >= 8) {
 544         lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
 545         cbuf[cnt++] = lastbyte >> lastbits;
 546         num_of_bits -= 8;
 547     }
 548     if (num_of_bits > 0) {
 549         lastbyte = (lastbyte << num_of_bits) | num;
 550         lastbits += num_of_bits;
 551         if (lastbits >= 8) {
 552             lastbits -= 8;
 553             cbuf[cnt++] = lastbyte >> lastbits;
 554         }
 555     }
 556     buf[0] = cnt;
 557     buf[1] = lastbits;
 558     buf[2] = lastbyte;
 559     if (lastbits>0) {
 560         cbuf[cnt] = lastbyte << (8 - lastbits);
 561     }
 562 }
 563
 564 /*_________________________________________________________________________
 565  |
 566  | sizeofint - calculate bitsize of an integer
 567  |
 568  | return the number of bits needed to store an integer with given max size
 569  |
 570 */
 571
 572 static int sizeofint(const int size) {
 573     unsigned int num = 1;
 574     int num_of_bits = 0;
 575
 576     while (size >= num && num_of_bits < 32) {
 577         num_of_bits++;
 578         num <<= 1;
 579     }
 580     return num_of_bits;
 581 }
 582
 583 /*___________________________________________________________________________
 584  |
 585  | sizeofints - calculate 'bitsize' of compressed ints
 586  |
 587  | given the number of small unsigned integers and the maximum value
 588  | return the number of bits needed to read or write them with the
 589  | routines receiveints and sendints. You need this parameter when
 590  | calling these routines. Note that for many calls I can use
 591  | the variable 'smallidx' which is exactly the number of bits, and
 592  | So I don't need to call 'sizeofints for those calls.
 593 */
 594
 595 static int sizeofints( const int num_of_ints, unsigned int sizes[]) {
 596     int i, num;
 597     unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
 598     num_of_bytes = 1;
 599     bytes[0] = 1;
 600     num_of_bits = 0;
 601     for (i=0; i < num_of_ints; i++) {
 602         tmp = 0;
 603         for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
 604             tmp = bytes[bytecnt] * sizes[i] + tmp;
 605             bytes[bytecnt] = tmp & 0xff;
 606             tmp >>= 8;
 607         }
 608         while (tmp != 0) {
 609             bytes[bytecnt++] = tmp & 0xff;
 610             tmp >>= 8;
 611         }
 612         num_of_bytes = bytecnt;
 613     }
 614     num = 1;
 615     num_of_bytes--;
 616     while (bytes[num_of_bytes] >= num) {
 617         num_of_bits++;
 618         num *= 2;
 619     }
 620     return num_of_bits + num_of_bytes * 8;
 621
 622 }
 623
 624 /*____________________________________________________________________________
 625  |
 626  | sendints - send a small set of small integers in compressed format
 627  |
 628  | this routine is used internally by xdr3dfcoord, to send a set of
 629  | small integers to the buffer.
 630  | Multiplication with fixed (specified maximum ) sizes is used to get
 631  | to one big, multibyte integer. Allthough the routine could be
 632  | modified to handle sizes bigger than 16777216, or more than just
 633  | a few integers, this is not done, because the gain in compression
 634  | isn't worth the effort. Note that overflowing the multiplication
 635  | or the byte buffer (32 bytes) is unchecked and causes bad results.
 636  |
 637  */
 638
 639 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
 640         unsigned int sizes[], unsigned int nums[]) {
 641
 642     int i;
 643     unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
 644
 645     tmp = nums[0];
 646     num_of_bytes = 0;
 647     do {
 648         bytes[num_of_bytes++] = tmp & 0xff;
 649         tmp >>= 8;
 650     } while (tmp != 0);
 651
 652     for (i = 1; i < num_of_ints; i++) {
 653         if (nums[i] >= sizes[i]) {
 654             fprintf(stderr,"major breakdown in sendints num %d doesn't "
 655                     "match size %d\n", nums[i], sizes[i]);
 656             exit(1);
 657         }
 658         /* use one step multiply */
 659         tmp = nums[i];
 660         for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
 661             tmp = bytes[bytecnt] * sizes[i] + tmp;
 662             bytes[bytecnt] = tmp & 0xff;
 663             tmp >>= 8;
 664         }
 665         while (tmp != 0) {
 666             bytes[bytecnt++] = tmp & 0xff;
 667             tmp >>= 8;
 668         }
 669         num_of_bytes = bytecnt;
 670     }
 671     if (num_of_bits >= num_of_bytes * 8) {
 672         for (i = 0; i < num_of_bytes; i++) {
 673             sendbits(buf, 8, bytes[i]);
 674         }
 675         sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
 676     } else {
 677         for (i = 0; i < num_of_bytes-1; i++) {
 678             sendbits(buf, 8, bytes[i]);
 679         }
 680         sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
 681     }
 682 }
 683
 684
 685 /*___________________________________________________________________________
 686  |
 687  | receivebits - decode number from buf using specified number of bits
 688  |
 689  | extract the number of bits from the array buf and construct an integer
 690  | from it. Return that value.
 691  |
 692 */
 693
 694 static int receivebits(int buf[], int num_of_bits) {
 695
 696     int cnt, num;
 697     unsigned int lastbits, lastbyte;
 698     unsigned char * cbuf;
 699     int mask = (1 << num_of_bits) -1;
 700
 701     cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
 702     cnt = buf[0];
 703     lastbits = (unsigned int) buf[1];
 704     lastbyte = (unsigned int) buf[2];
 705
 706     num = 0;
 707     while (num_of_bits >= 8) {
 708         lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
 709         num |=  (lastbyte >> lastbits) << (num_of_bits - 8);
 710         num_of_bits -=8;
 711     }
 712     if (num_of_bits > 0) {
 713         if (lastbits < num_of_bits) {
 714             lastbits += 8;
 715             lastbyte = (lastbyte << 8) | cbuf[cnt++];
 716         }
 717         lastbits -= num_of_bits;
 718         num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
 719     }
 720     num &= mask;
 721     buf[0] = cnt;
 722     buf[1] = lastbits;
 723     buf[2] = lastbyte;
 724     return num;
 725 }
 726
 727 /*____________________________________________________________________________
 728  |
 729  | receiveints - decode 'small' integers from the buf array
 730  |
 731  | this routine is the inverse from sendints() and decodes the small integers
 732  | written to buf by calculating the remainder and doing divisions with
 733  | the given sizes[]. You need to specify the total number of bits to be
 734  | used from buf in num_of_bits.
 735  |
 736 */
 737
 738 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
 739         unsigned int sizes[], int nums[]) {
 740     int bytes[32];
 741     int i, j, num_of_bytes, p, num;
 742
 743     bytes[1] = bytes[2] = bytes[3] = 0;
 744     num_of_bytes = 0;
 745     while (num_of_bits > 8) {
 746         bytes[num_of_bytes++] = receivebits(buf, 8);
 747         num_of_bits -= 8;
 748     }
 749     if (num_of_bits > 0) {
 750         bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
 751     }
 752     for (i = num_of_ints-1; i > 0; i--) {
 753         num = 0;
 754         for (j = num_of_bytes-1; j >=0; j--) {
 755             num = (num << 8) | bytes[j];
 756             p = num / sizes[i];
 757             bytes[j] = p;
 758             num = num - p * sizes[i];
 759         }
 760         nums[i] = num;
 761     }
 762     nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
 763 }
 764
 765 /*____________________________________________________________________________
 766  |
 767  | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
 768  |
 769  | this routine reads or writes (depending on how you opened the file with
 770  | xdropen() ) a large number of 3d coordinates (stored in *fp).
 771  | The number of coordinates triplets to write is given by *size. On
 772  | read this number may be zero, in which case it reads as many as were written
 773  | or it may specify the number if triplets to read (which should match the
 774  | number written).
 775  | Compression is achieved by first converting all floating numbers to integer
 776  | using multiplication by *precision and rounding to the nearest integer.
 777  | Then the minimum and maximum value are calculated to determine the range.
 778  | The limited range of integers so found, is used to compress the coordinates.
 779  | In addition the differences between succesive coordinates is calculated.
 780  | If the difference happens to be 'small' then only the difference is saved,
 781  | compressing the data even more. The notion of 'small' is changed dynamically
 782  | and is enlarged or reduced whenever needed or possible.
 783  | Extra compression is achieved in the case of GROMOS and coordinates of
 784  | water molecules. GROMOS first writes out the Oxygen position, followed by
 785  | the two hydrogens. In order to make the differences smaller (and thereby
 786  | compression the data better) the order is changed into first one hydrogen
 787  | then the oxygen, followed by the other hydrogen. This is rather special, but
 788  | it shouldn't harm in the general case.
 789  |
 790  */
 791
 792 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
 793
 794
 795     static int *ip = NULL;
 796     static int oldsize;
 797     static int *buf;
 798
 799     int minint[3], maxint[3], mindiff, *lip, diff;
 800     int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
 801     int minidx, maxidx;
 802     unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
 803     int flag, k;
 804     int small, smaller, larger, i, is_small, is_smaller, run, prevrun;
 805     float *lfp, lf;
 806     int tmp, *thiscoord,  prevcoord[3];
 807     unsigned int tmpcoord[30];
 808
 809     int bufsize, xdrid, lsize;
 810     unsigned int bitsize;
 811     float inv_precision;
 812     int errval = 1;
 813
 814     /* find out if xdrs is opened for reading or for writing */
 815     xdrid = 0;
 816     while (xdridptr[xdrid] != xdrs) {
 817         xdrid++;
 818         if (xdrid >= MAXID) {
 819             fprintf(stderr, "xdr error. no open xdr stream\n");
 820             exit (1);
 821         }
 822     }
 823     if (xdrmodes[xdrid] == 'w') {
 824
 825         /* xdrs is open for writing */
 826
 827         if (xdr_int(xdrs, size) == 0)
 828             return 0;
 829         size3 = *size * 3;
 830         /* when the number of coordinates is small, don't try to compress; just
 831          * write them as floats using xdr_vector
 832          */
 833         if (*size <= 9 ) {
 834             return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
 835                 (xdrproc_t)xdr_float));
 836         }
 837
 838         xdr_float(xdrs, precision);
 839         if (ip == NULL) {
 840             ip = (int *)malloc(size3 * sizeof(*ip));
 841             if (ip == NULL) {
 842                 fprintf(stderr,"malloc failed\n");
 843                 exit(1);
 844             }
 845             bufsize = size3 * 1.2;
 846             buf = (int *)malloc(bufsize * sizeof(*buf));
 847             if (buf == NULL) {
 848                 fprintf(stderr,"malloc failed\n");
 849                 exit(1);
 850             }
 851             oldsize = *size;
 852         } else if (*size > oldsize) {
 853             ip = (int *)realloc(ip, size3 * sizeof(*ip));
 854             if (ip == NULL) {
 855                 fprintf(stderr,"malloc failed\n");
 856                 exit(1);
 857             }
 858             bufsize = size3 * 1.2;
 859             buf = (int *)realloc(buf, bufsize * sizeof(*buf));
 860             if (buf == NULL) {
 861                 fprintf(stderr,"malloc failed\n");
 862                 exit(1);
 863             }
 864             oldsize = *size;
 865         }
 866         /* buf[0-2] are special and do not contain actual data */
 867         buf[0] = buf[1] = buf[2] = 0;
 868         minint[0] = minint[1] = minint[2] = INT_MAX;
 869         maxint[0] = maxint[1] = maxint[2] = INT_MIN;
 870         prevrun = -1;
 871         lfp = fp;
 872         lip = ip;
 873         mindiff = INT_MAX;
 874         oldlint1 = oldlint2 = oldlint3 = 0;
 875         while(lfp < fp + size3 ) {
 876             /* find nearest integer */
 877             if (*lfp >= 0.0)
 878                 lf = *lfp * *precision + 0.5;
 879             else
 880                 lf = *lfp * *precision - 0.5;
 881             if (fabs(lf) > MAXABS) {
 882                 /* scaling would cause overflow */
 883                 errval = 0;
 884             }
 885             lint1 = lf;
 886             if (lint1 < minint[0]) minint[0] = lint1;
 887             if (lint1 > maxint[0]) maxint[0] = lint1;
 888             *lip++ = lint1;
 889             lfp++;
 890             if (*lfp >= 0.0)
 891                 lf = *lfp * *precision + 0.5;
 892             else
 893                 lf = *lfp * *precision - 0.5;
 894             if (fabs(lf) > MAXABS) {
 895                 /* scaling would cause overflow */
 896                 errval = 0;
 897             }
 898             lint2 = lf;
 899             if (lint2 < minint[1]) minint[1] = lint2;
 900             if (lint2 > maxint[1]) maxint[1] = lint2;
 901             *lip++ = lint2;
 902             lfp++;
 903             if (*lfp >= 0.0)
 904                 lf = *lfp * *precision + 0.5;
 905             else
 906                 lf = *lfp * *precision - 0.5;
 907             if (fabs(lf) > MAXABS) {
 908                 /* scaling would cause overflow */
 909                 errval = 0;
 910             }
 911             lint3 = lf;
 912             if (lint3 < minint[2]) minint[2] = lint3;
 913             if (lint3 > maxint[2]) maxint[2] = lint3;
 914             *lip++ = lint3;
 915             lfp++;
 916             diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
 917             if (diff < mindiff && lfp > fp + 3)
 918                 mindiff = diff;
 919             oldlint1 = lint1;
 920             oldlint2 = lint2;
 921             oldlint3 = lint3;
 922         }
 923         xdr_int(xdrs, &(minint[0]));
 924         xdr_int(xdrs, &(minint[1]));
 925         xdr_int(xdrs, &(minint[2]));
 926
 927         xdr_int(xdrs, &(maxint[0]));
 928         xdr_int(xdrs, &(maxint[1]));
 929         xdr_int(xdrs, &(maxint[2]));
 930
 931         if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
 932                 (float)maxint[1] - (float)minint[1] >= MAXABS ||
 933                 (float)maxint[2] - (float)minint[2] >= MAXABS) {
 934             /* turning value in unsigned by subtracting minint
 935              * would cause overflow
 936              */
 937             errval = 0;
 938         }
 939         sizeint[0] = maxint[0] - minint[0]+1;
 940         sizeint[1] = maxint[1] - minint[1]+1;
 941         sizeint[2] = maxint[2] - minint[2]+1;
 942
 943         /* check if one of the sizes is to big to be multiplied */
 944         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
 945             bitsizeint[0] = sizeofint(sizeint[0]);
 946             bitsizeint[1] = sizeofint(sizeint[1]);
 947             bitsizeint[2] = sizeofint(sizeint[2]);
 948             bitsize = 0; /* flag the use of large sizes */
 949         } else {
 950             bitsize = sizeofints(3, sizeint);
 951         }
 952         lip = ip;
 953         luip = (unsigned int *) ip;
 954         smallidx = FIRSTIDX;
 955         while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
 956             smallidx++;
 957         }
 958         xdr_int(xdrs, &smallidx);
 959         maxidx = MIN(LASTIDX, smallidx + 8) ;
 960         minidx = maxidx - 8; /* often this equal smallidx */
 961         smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
 962         small = magicints[smallidx] / 2;
 963         sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
 964         larger = magicints[maxidx] / 2;
 965         i = 0;
 966         while (i < *size) {
 967             is_small = 0;
 968             thiscoord = (int *)(luip) + i * 3;
 969             if (smallidx < maxidx && i >= 1 &&
 970                     abs(thiscoord[0] - prevcoord[0]) < larger &&
 971                     abs(thiscoord[1] - prevcoord[1]) < larger &&
 972                     abs(thiscoord[2] - prevcoord[2]) < larger) {
 973                 is_smaller = 1;
 974             } else if (smallidx > minidx) {
 975                 is_smaller = -1;
 976             } else {
 977                 is_smaller = 0;
 978             }
 979             if (i + 1 < *size) {
 980                 if (abs(thiscoord[0] - thiscoord[3]) < small &&
 981                         abs(thiscoord[1] - thiscoord[4]) < small &&
 982                         abs(thiscoord[2] - thiscoord[5]) < small) {
 983                     /* interchange first with second atom for better
 984                      * compression of water molecules
 985                      */
 986                     tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
 987                         thiscoord[3] = tmp;
 988                     tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
 989                         thiscoord[4] = tmp;
 990                     tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
 991                         thiscoord[5] = tmp;
 992                     is_small = 1;
 993                 }
 994
 995             }
 996             tmpcoord[0] = thiscoord[0] - minint[0];
 997             tmpcoord[1] = thiscoord[1] - minint[1];
 998             tmpcoord[2] = thiscoord[2] - minint[2];
 999             if (bitsize == 0) {
1000                 sendbits(buf, bitsizeint[0], tmpcoord[0]);
1001                 sendbits(buf, bitsizeint[1], tmpcoord[1]);
1002                 sendbits(buf, bitsizeint[2], tmpcoord[2]);
1003             } else {
1004                 sendints(buf, 3, bitsize, sizeint, tmpcoord);
1005             }
1006             prevcoord[0] = thiscoord[0];
1007             prevcoord[1] = thiscoord[1];
1008             prevcoord[2] = thiscoord[2];
1009             thiscoord = thiscoord + 3;
1010             i++;
1011
1012             run = 0;
1013             if (is_small == 0 && is_smaller == -1)
1014                 is_smaller = 0;
1015             while (is_small && run < 8*3) {
1016                 if (is_smaller == -1 && (
1017                         SQR(thiscoord[0] - prevcoord[0]) +
1018                         SQR(thiscoord[1] - prevcoord[1]) +
1019                         SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) {
1020                     is_smaller = 0;
1021                 }
1022
1023                 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small;
1024                 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small;
1025                 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small;
1026
1027                 prevcoord[0] = thiscoord[0];
1028                 prevcoord[1] = thiscoord[1];
1029                 prevcoord[2] = thiscoord[2];
1030
1031                 i++;
1032                 thiscoord = thiscoord + 3;
1033                 is_small = 0;
1034                 if (i < *size &&
1035                         abs(thiscoord[0] - prevcoord[0]) < small &&
1036                         abs(thiscoord[1] - prevcoord[1]) < small &&
1037                         abs(thiscoord[2] - prevcoord[2]) < small) {
1038                     is_small = 1;
1039                 }
1040             }
1041             if (run != prevrun || is_smaller != 0) {
1042                 prevrun = run;
1043                 sendbits(buf, 1, 1); /* flag the change in run-length */
1044                 sendbits(buf, 5, run+is_smaller+1);
1045             } else {
1046                 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
1047             }
1048             for (k=0; k < run; k+=3) {
1049                 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
1050             }
1051             if (is_smaller != 0) {
1052                 smallidx += is_smaller;
1053                 if (is_smaller < 0) {
1054                     small = smaller;
1055                     smaller = magicints[smallidx-1] / 2;
1056                 } else {
1057                     smaller = small;
1058                     small = magicints[smallidx] / 2;
1059                 }
1060                 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1061             }
1062         }
1063         if (buf[1] != 0) buf[0]++;;
1064         xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */
1065         return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]));
1066     } else {
1067
1068         /* xdrs is open for reading */
1069
1070         if (xdr_int(xdrs, &lsize) == 0)
1071             return 0;
1072         if (*size != 0 && lsize != *size) {
1073             fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; "
1074                     "%d arg vs %d in file", *size, lsize);
1075         }
1076         *size = lsize;
1077         size3 = *size * 3;
1078         if (*size <= 9) {
1079             return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
1080                 (xdrproc_t)xdr_float));
1081         }
1082         xdr_float(xdrs, precision);
1083         if (ip == NULL) {
1084             ip = (int *)malloc(size3 * sizeof(*ip));
1085             if (ip == NULL) {
1086                 fprintf(stderr,"malloc failed\n");
1087                 exit(1);
1088             }
1089             bufsize = size3 * 1.2;
1090             buf = (int *)malloc(bufsize * sizeof(*buf));
1091             if (buf == NULL) {
1092                 fprintf(stderr,"malloc failed\n");
1093                 exit(1);
1094             }
1095             oldsize = *size;
1096         } else if (*size > oldsize) {
1097             ip = (int *)realloc(ip, size3 * sizeof(*ip));
1098             if (ip == NULL) {
1099                 fprintf(stderr,"malloc failed\n");
1100                 exit(1);
1101             }
1102             bufsize = size3 * 1.2;
1103             buf = (int *)realloc(buf, bufsize * sizeof(*buf));
1104             if (buf == NULL) {
1105                 fprintf(stderr,"malloc failed\n");
1106                 exit(1);
1107             }
1108             oldsize = *size;
1109         }
1110         buf[0] = buf[1] = buf[2] = 0;
1111
1112         xdr_int(xdrs, &(minint[0]));
1113         xdr_int(xdrs, &(minint[1]));
1114         xdr_int(xdrs, &(minint[2]));
1115
1116         xdr_int(xdrs, &(maxint[0]));
1117         xdr_int(xdrs, &(maxint[1]));
1118         xdr_int(xdrs, &(maxint[2]));
1119
1120         sizeint[0] = maxint[0] - minint[0]+1;
1121         sizeint[1] = maxint[1] - minint[1]+1;
1122         sizeint[2] = maxint[2] - minint[2]+1;
1123
1124         /* check if one of the sizes is to big to be multiplied */
1125         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
1126             bitsizeint[0] = sizeofint(sizeint[0]);
1127             bitsizeint[1] = sizeofint(sizeint[1]);
1128             bitsizeint[2] = sizeofint(sizeint[2]);
1129             bitsize = 0; /* flag the use of large sizes */
1130         } else {
1131             bitsize = sizeofints(3, sizeint);
1132         }
1133
1134         xdr_int(xdrs, &smallidx);
1135         maxidx = MIN(LASTIDX, smallidx + 8) ;
1136         minidx = maxidx - 8; /* often this equal smallidx */
1137         smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
1138         small = magicints[smallidx] / 2;
1139         sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1140         larger = magicints[maxidx];
1141
1142         /* buf[0] holds the length in bytes */
1143
1144         if (xdr_int(xdrs, &(buf[0])) == 0)
1145             return 0;
1146         if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0)
1147             return 0;
1148         buf[0] = buf[1] = buf[2] = 0;
1149
1150         lfp = fp;
1151         inv_precision = 1.0 / * precision;
1152         run = 0;
1153         i = 0;
1154         lip = ip;
1155         while ( i < lsize ) {
1156             thiscoord = (int *)(lip) + i * 3;
1157
1158             if (bitsize == 0) {
1159                 thiscoord[0] = receivebits(buf, bitsizeint[0]);
1160                 thiscoord[1] = receivebits(buf, bitsizeint[1]);
1161                 thiscoord[2] = receivebits(buf, bitsizeint[2]);
1162             } else {
1163                 receiveints(buf, 3, bitsize, sizeint, thiscoord);
1164             }
1165
1166             i++;
1167             thiscoord[0] += minint[0];
1168             thiscoord[1] += minint[1];
1169             thiscoord[2] += minint[2];
1170
1171             prevcoord[0] = thiscoord[0];
1172             prevcoord[1] = thiscoord[1];
1173             prevcoord[2] = thiscoord[2];
1174
1175
1176             flag = receivebits(buf, 1);
1177             is_smaller = 0;
1178             if (flag == 1) {
1179                 run = receivebits(buf, 5);
1180                 is_smaller = run % 3;
1181                 run -= is_smaller;
1182                 is_smaller--;
1183             }
1184             if (run > 0) {
1185                 thiscoord += 3;
1186                 for (k = 0; k < run; k+=3) {
1187                     receiveints(buf, 3, smallidx, sizesmall, thiscoord);
1188                     i++;
1189                     thiscoord[0] += prevcoord[0] - small;
1190                     thiscoord[1] += prevcoord[1] - small;
1191                     thiscoord[2] += prevcoord[2] - small;
1192                     if (k == 0) {
1193                         /* interchange first with second atom for better
1194                          * compression of water molecules
1195                          */
1196                         tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1197                                 prevcoord[0] = tmp;
1198                         tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1199                                 prevcoord[1] = tmp;
1200                         tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1201                                 prevcoord[2] = tmp;
1202                         *lfp++ = prevcoord[0] * inv_precision;
1203                         *lfp++ = prevcoord[1] * inv_precision;
1204                         *lfp++ = prevcoord[2] * inv_precision;
1205                     } else {
1206                         prevcoord[0] = thiscoord[0];
1207                         prevcoord[1] = thiscoord[1];
1208                         prevcoord[2] = thiscoord[2];
1209                     }
1210                     *lfp++ = thiscoord[0] * inv_precision;
1211                     *lfp++ = thiscoord[1] * inv_precision;
1212                     *lfp++ = thiscoord[2] * inv_precision;
1213                 }
1214             } else {
1215                 *lfp++ = thiscoord[0] * inv_precision;
1216                 *lfp++ = thiscoord[1] * inv_precision;
1217                 *lfp++ = thiscoord[2] * inv_precision;
1218             }
1219             smallidx += is_smaller;
1220             if (is_smaller < 0) {
1221                 small = smaller;
1222                 if (smallidx > FIRSTIDX) {
1223                     smaller = magicints[smallidx - 1] /2;
1224                 } else {
1225                     smaller = 0;
1226                 }
1227             } else if (is_smaller > 0) {
1228                 smaller = small;
1229                 small = magicints[smallidx] / 2;
1230             }
1231             sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1232         }
1233     }
1234     return 1;
1235 }
1236
1237
1238