source/unres/src_MD/src/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 = "a+";
 460             lmode = XDR_ENCODE;
 461     } else {
 462             type = "r";
 463             type1 = "r";
 464             lmode = XDR_DECODE;
 465     }
 466     xdrfiles[xdrid] = fopen(filename, type1);
 467     if (xdrfiles[xdrid] == NULL) {
 468         xdrs = NULL;
 469         return 0;
 470     }
 471     xdrmodes[xdrid] = *type;
 472     /* next test isn't usefull in the case of C language
 473      * but is used for the Fortran interface
 474      * (C users are expected to pass the address of an already allocated
 475      * XDR staructure)
 476      */
 477     if (xdrs == NULL) {
 478         xdridptr[xdrid] = (XDR *) malloc(sizeof(XDR));
 479         xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode);
 480     } else {
 481         xdridptr[xdrid] = xdrs;
 482         xdrstdio_create(xdrs, xdrfiles[xdrid], lmode);
 483     }
 484     return xdrid;
 485 }
 486
 487 /*_________________________________________________________________________
 488  |
 489  | xdrclose - close a xdr file
 490  |
 491  | This will flush the xdr buffers, and destroy the xdr stream.
 492  | It also closes the associated file descriptor (this is *not*
 493  | done by xdr_destroy).
 494  |
 495 */
 496
 497 int xdrclose(XDR *xdrs) {
 498     int xdrid;
 499
 500     if (xdrs == NULL) {
 501         fprintf(stderr, "xdrclose: passed a NULL pointer\n");
 502         exit(1);
 503     }
 504     for (xdrid = 1; xdrid < MAXID; xdrid++) {
 505         if (xdridptr[xdrid] == xdrs) {
 506
 507             xdr_destroy(xdrs);
 508             fclose(xdrfiles[xdrid]);
 509             xdridptr[xdrid] = NULL;
 510             return 1;
 511         }
 512     }
 513     fprintf(stderr, "xdrclose: no such open xdr file\n");
 514     exit(1);
 515
 516 }
 517
 518 /*____________________________________________________________________________
 519  |
 520  | sendbits - encode num into buf using the specified number of bits
 521  |
 522  | This routines appends the value of num to the bits already present in
 523  | the array buf. You need to give it the number of bits to use and you
 524  | better make sure that this number of bits is enough to hold the value
 525  | Also num must be positive.
 526  |
 527 */
 528
 529 static void sendbits(int buf[], int num_of_bits, int num) {
 530
 531     unsigned int cnt, lastbyte;
 532     int lastbits;
 533     unsigned char * cbuf;
 534
 535     cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
 536     cnt = (unsigned int) buf[0];
 537     lastbits = buf[1];
 538     lastbyte =(unsigned int) buf[2];
 539     while (num_of_bits >= 8) {
 540         lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
 541         cbuf[cnt++] = lastbyte >> lastbits;
 542         num_of_bits -= 8;
 543     }
 544     if (num_of_bits > 0) {
 545         lastbyte = (lastbyte << num_of_bits) | num;
 546         lastbits += num_of_bits;
 547         if (lastbits >= 8) {
 548             lastbits -= 8;
 549             cbuf[cnt++] = lastbyte >> lastbits;
 550         }
 551     }
 552     buf[0] = cnt;
 553     buf[1] = lastbits;
 554     buf[2] = lastbyte;
 555     if (lastbits>0) {
 556         cbuf[cnt] = lastbyte << (8 - lastbits);
 557     }
 558 }
 559
 560 /*_________________________________________________________________________
 561  |
 562  | sizeofint - calculate bitsize of an integer
 563  |
 564  | return the number of bits needed to store an integer with given max size
 565  |
 566 */
 567
 568 static int sizeofint(const int size) {
 569     unsigned int num = 1;
 570     int num_of_bits = 0;
 571
 572     while (size >= num && num_of_bits < 32) {
 573         num_of_bits++;
 574         num <<= 1;
 575     }
 576     return num_of_bits;
 577 }
 578
 579 /*___________________________________________________________________________
 580  |
 581  | sizeofints - calculate 'bitsize' of compressed ints
 582  |
 583  | given the number of small unsigned integers and the maximum value
 584  | return the number of bits needed to read or write them with the
 585  | routines receiveints and sendints. You need this parameter when
 586  | calling these routines. Note that for many calls I can use
 587  | the variable 'smallidx' which is exactly the number of bits, and
 588  | So I don't need to call 'sizeofints for those calls.
 589 */
 590
 591 static int sizeofints( const int num_of_ints, unsigned int sizes[]) {
 592     int i, num;
 593     unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
 594     num_of_bytes = 1;
 595     bytes[0] = 1;
 596     num_of_bits = 0;
 597     for (i=0; i < num_of_ints; i++) {
 598         tmp = 0;
 599         for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
 600             tmp = bytes[bytecnt] * sizes[i] + tmp;
 601             bytes[bytecnt] = tmp & 0xff;
 602             tmp >>= 8;
 603         }
 604         while (tmp != 0) {
 605             bytes[bytecnt++] = tmp & 0xff;
 606             tmp >>= 8;
 607         }
 608         num_of_bytes = bytecnt;
 609     }
 610     num = 1;
 611     num_of_bytes--;
 612     while (bytes[num_of_bytes] >= num) {
 613         num_of_bits++;
 614         num *= 2;
 615     }
 616     return num_of_bits + num_of_bytes * 8;
 617
 618 }
 619
 620 /*____________________________________________________________________________
 621  |
 622  | sendints - send a small set of small integers in compressed format
 623  |
 624  | this routine is used internally by xdr3dfcoord, to send a set of
 625  | small integers to the buffer.
 626  | Multiplication with fixed (specified maximum ) sizes is used to get
 627  | to one big, multibyte integer. Allthough the routine could be
 628  | modified to handle sizes bigger than 16777216, or more than just
 629  | a few integers, this is not done, because the gain in compression
 630  | isn't worth the effort. Note that overflowing the multiplication
 631  | or the byte buffer (32 bytes) is unchecked and causes bad results.
 632  |
 633  */
 634
 635 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
 636         unsigned int sizes[], unsigned int nums[]) {
 637
 638     int i;
 639     unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
 640
 641     tmp = nums[0];
 642     num_of_bytes = 0;
 643     do {
 644         bytes[num_of_bytes++] = tmp & 0xff;
 645         tmp >>= 8;
 646     } while (tmp != 0);
 647
 648     for (i = 1; i < num_of_ints; i++) {
 649         if (nums[i] >= sizes[i]) {
 650             fprintf(stderr,"major breakdown in sendints num %d doesn't "
 651                     "match size %d\n", nums[i], sizes[i]);
 652             exit(1);
 653         }
 654         /* use one step multiply */
 655         tmp = nums[i];
 656         for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
 657             tmp = bytes[bytecnt] * sizes[i] + tmp;
 658             bytes[bytecnt] = tmp & 0xff;
 659             tmp >>= 8;
 660         }
 661         while (tmp != 0) {
 662             bytes[bytecnt++] = tmp & 0xff;
 663             tmp >>= 8;
 664         }
 665         num_of_bytes = bytecnt;
 666     }
 667     if (num_of_bits >= num_of_bytes * 8) {
 668         for (i = 0; i < num_of_bytes; i++) {
 669             sendbits(buf, 8, bytes[i]);
 670         }
 671         sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
 672     } else {
 673         for (i = 0; i < num_of_bytes-1; i++) {
 674             sendbits(buf, 8, bytes[i]);
 675         }
 676         sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
 677     }
 678 }
 679
 680
 681 /*___________________________________________________________________________
 682  |
 683  | receivebits - decode number from buf using specified number of bits
 684  |
 685  | extract the number of bits from the array buf and construct an integer
 686  | from it. Return that value.
 687  |
 688 */
 689
 690 static int receivebits(int buf[], int num_of_bits) {
 691
 692     int cnt, num;
 693     unsigned int lastbits, lastbyte;
 694     unsigned char * cbuf;
 695     int mask = (1 << num_of_bits) -1;
 696
 697     cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
 698     cnt = buf[0];
 699     lastbits = (unsigned int) buf[1];
 700     lastbyte = (unsigned int) buf[2];
 701
 702     num = 0;
 703     while (num_of_bits >= 8) {
 704         lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
 705         num |=  (lastbyte >> lastbits) << (num_of_bits - 8);
 706         num_of_bits -=8;
 707     }
 708     if (num_of_bits > 0) {
 709         if (lastbits < num_of_bits) {
 710             lastbits += 8;
 711             lastbyte = (lastbyte << 8) | cbuf[cnt++];
 712         }
 713         lastbits -= num_of_bits;
 714         num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
 715     }
 716     num &= mask;
 717     buf[0] = cnt;
 718     buf[1] = lastbits;
 719     buf[2] = lastbyte;
 720     return num;
 721 }
 722
 723 /*____________________________________________________________________________
 724  |
 725  | receiveints - decode 'small' integers from the buf array
 726  |
 727  | this routine is the inverse from sendints() and decodes the small integers
 728  | written to buf by calculating the remainder and doing divisions with
 729  | the given sizes[]. You need to specify the total number of bits to be
 730  | used from buf in num_of_bits.
 731  |
 732 */
 733
 734 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
 735         unsigned int sizes[], int nums[]) {
 736     int bytes[32];
 737     int i, j, num_of_bytes, p, num;
 738
 739     bytes[1] = bytes[2] = bytes[3] = 0;
 740     num_of_bytes = 0;
 741     while (num_of_bits > 8) {
 742         bytes[num_of_bytes++] = receivebits(buf, 8);
 743         num_of_bits -= 8;
 744     }
 745     if (num_of_bits > 0) {
 746         bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
 747     }
 748     for (i = num_of_ints-1; i > 0; i--) {
 749         num = 0;
 750         for (j = num_of_bytes-1; j >=0; j--) {
 751             num = (num << 8) | bytes[j];
 752             p = num / sizes[i];
 753             bytes[j] = p;
 754             num = num - p * sizes[i];
 755         }
 756         nums[i] = num;
 757     }
 758     nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
 759 }
 760
 761 /*____________________________________________________________________________
 762  |
 763  | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
 764  |
 765  | this routine reads or writes (depending on how you opened the file with
 766  | xdropen() ) a large number of 3d coordinates (stored in *fp).
 767  | The number of coordinates triplets to write is given by *size. On
 768  | read this number may be zero, in which case it reads as many as were written
 769  | or it may specify the number if triplets to read (which should match the
 770  | number written).
 771  | Compression is achieved by first converting all floating numbers to integer
 772  | using multiplication by *precision and rounding to the nearest integer.
 773  | Then the minimum and maximum value are calculated to determine the range.
 774  | The limited range of integers so found, is used to compress the coordinates.
 775  | In addition the differences between succesive coordinates is calculated.
 776  | If the difference happens to be 'small' then only the difference is saved,
 777  | compressing the data even more. The notion of 'small' is changed dynamically
 778  | and is enlarged or reduced whenever needed or possible.
 779  | Extra compression is achieved in the case of GROMOS and coordinates of
 780  | water molecules. GROMOS first writes out the Oxygen position, followed by
 781  | the two hydrogens. In order to make the differences smaller (and thereby
 782  | compression the data better) the order is changed into first one hydrogen
 783  | then the oxygen, followed by the other hydrogen. This is rather special, but
 784  | it shouldn't harm in the general case.
 785  |
 786  */
 787
 788 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
 789
 790
 791     static int *ip = NULL;
 792     static int oldsize;
 793     static int *buf;
 794
 795     int minint[3], maxint[3], mindiff, *lip, diff;
 796     int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
 797     int minidx, maxidx;
 798     unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
 799     int flag, k;
 800     int small, smaller, larger, i, is_small, is_smaller, run, prevrun;
 801     float *lfp, lf;
 802     int tmp, *thiscoord,  prevcoord[3];
 803     unsigned int tmpcoord[30];
 804
 805     int bufsize, xdrid, lsize;
 806     unsigned int bitsize;
 807     float inv_precision;
 808     int errval = 1;
 809
 810     /* find out if xdrs is opened for reading or for writing */
 811     xdrid = 0;
 812     while (xdridptr[xdrid] != xdrs) {
 813         xdrid++;
 814         if (xdrid >= MAXID) {
 815             fprintf(stderr, "xdr error. no open xdr stream\n");
 816             exit (1);
 817         }
 818     }
 819     if (xdrmodes[xdrid] == 'w') {
 820
 821         /* xdrs is open for writing */
 822
 823         if (xdr_int(xdrs, size) == 0)
 824             return 0;
 825         size3 = *size * 3;
 826         /* when the number of coordinates is small, don't try to compress; just
 827          * write them as floats using xdr_vector
 828          */
 829         if (*size <= 9 ) {
 830             return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
 831                 (xdrproc_t)xdr_float));
 832         }
 833
 834         xdr_float(xdrs, precision);
 835         if (ip == NULL) {
 836             ip = (int *)malloc(size3 * sizeof(*ip));
 837             if (ip == NULL) {
 838                 fprintf(stderr,"malloc failed\n");
 839                 exit(1);
 840             }
 841             bufsize = size3 * 1.2;
 842             buf = (int *)malloc(bufsize * sizeof(*buf));
 843             if (buf == NULL) {
 844                 fprintf(stderr,"malloc failed\n");
 845                 exit(1);
 846             }
 847             oldsize = *size;
 848         } else if (*size > oldsize) {
 849             ip = (int *)realloc(ip, size3 * sizeof(*ip));
 850             if (ip == NULL) {
 851                 fprintf(stderr,"malloc failed\n");
 852                 exit(1);
 853             }
 854             bufsize = size3 * 1.2;
 855             buf = (int *)realloc(buf, bufsize * sizeof(*buf));
 856             if (buf == NULL) {
 857                 fprintf(stderr,"malloc failed\n");
 858                 exit(1);
 859             }
 860             oldsize = *size;
 861         }
 862         /* buf[0-2] are special and do not contain actual data */
 863         buf[0] = buf[1] = buf[2] = 0;
 864         minint[0] = minint[1] = minint[2] = INT_MAX;
 865         maxint[0] = maxint[1] = maxint[2] = INT_MIN;
 866         prevrun = -1;
 867         lfp = fp;
 868         lip = ip;
 869         mindiff = INT_MAX;
 870         oldlint1 = oldlint2 = oldlint3 = 0;
 871         while(lfp < fp + size3 ) {
 872             /* find nearest integer */
 873             if (*lfp >= 0.0)
 874                 lf = *lfp * *precision + 0.5;
 875             else
 876                 lf = *lfp * *precision - 0.5;
 877             if (fabs(lf) > MAXABS) {
 878                 /* scaling would cause overflow */
 879                 errval = 0;
 880             }
 881             lint1 = lf;
 882             if (lint1 < minint[0]) minint[0] = lint1;
 883             if (lint1 > maxint[0]) maxint[0] = lint1;
 884             *lip++ = lint1;
 885             lfp++;
 886             if (*lfp >= 0.0)
 887                 lf = *lfp * *precision + 0.5;
 888             else
 889                 lf = *lfp * *precision - 0.5;
 890             if (fabs(lf) > MAXABS) {
 891                 /* scaling would cause overflow */
 892                 errval = 0;
 893             }
 894             lint2 = lf;
 895             if (lint2 < minint[1]) minint[1] = lint2;
 896             if (lint2 > maxint[1]) maxint[1] = lint2;
 897             *lip++ = lint2;
 898             lfp++;
 899             if (*lfp >= 0.0)
 900                 lf = *lfp * *precision + 0.5;
 901             else
 902                 lf = *lfp * *precision - 0.5;
 903             if (fabs(lf) > MAXABS) {
 904                 /* scaling would cause overflow */
 905                 errval = 0;
 906             }
 907             lint3 = lf;
 908             if (lint3 < minint[2]) minint[2] = lint3;
 909             if (lint3 > maxint[2]) maxint[2] = lint3;
 910             *lip++ = lint3;
 911             lfp++;
 912             diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
 913             if (diff < mindiff && lfp > fp + 3)
 914                 mindiff = diff;
 915             oldlint1 = lint1;
 916             oldlint2 = lint2;
 917             oldlint3 = lint3;
 918         }
 919         xdr_int(xdrs, &(minint[0]));
 920         xdr_int(xdrs, &(minint[1]));
 921         xdr_int(xdrs, &(minint[2]));
 922
 923         xdr_int(xdrs, &(maxint[0]));
 924         xdr_int(xdrs, &(maxint[1]));
 925         xdr_int(xdrs, &(maxint[2]));
 926
 927         if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
 928                 (float)maxint[1] - (float)minint[1] >= MAXABS ||
 929                 (float)maxint[2] - (float)minint[2] >= MAXABS) {
 930             /* turning value in unsigned by subtracting minint
 931              * would cause overflow
 932              */
 933             errval = 0;
 934         }
 935         sizeint[0] = maxint[0] - minint[0]+1;
 936         sizeint[1] = maxint[1] - minint[1]+1;
 937         sizeint[2] = maxint[2] - minint[2]+1;
 938
 939         /* check if one of the sizes is to big to be multiplied */
 940         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
 941             bitsizeint[0] = sizeofint(sizeint[0]);
 942             bitsizeint[1] = sizeofint(sizeint[1]);
 943             bitsizeint[2] = sizeofint(sizeint[2]);
 944             bitsize = 0; /* flag the use of large sizes */
 945         } else {
 946             bitsize = sizeofints(3, sizeint);
 947         }
 948         lip = ip;
 949         luip = (unsigned int *) ip;
 950         smallidx = FIRSTIDX;
 951         while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
 952             smallidx++;
 953         }
 954         xdr_int(xdrs, &smallidx);
 955         maxidx = MIN(LASTIDX, smallidx + 8) ;
 956         minidx = maxidx - 8; /* often this equal smallidx */
 957         smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
 958         small = magicints[smallidx] / 2;
 959         sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
 960         larger = magicints[maxidx] / 2;
 961         i = 0;
 962         while (i < *size) {
 963             is_small = 0;
 964             thiscoord = (int *)(luip) + i * 3;
 965             if (smallidx < maxidx && i >= 1 &&
 966                     abs(thiscoord[0] - prevcoord[0]) < larger &&
 967                     abs(thiscoord[1] - prevcoord[1]) < larger &&
 968                     abs(thiscoord[2] - prevcoord[2]) < larger) {
 969                 is_smaller = 1;
 970             } else if (smallidx > minidx) {
 971                 is_smaller = -1;
 972             } else {
 973                 is_smaller = 0;
 974             }
 975             if (i + 1 < *size) {
 976                 if (abs(thiscoord[0] - thiscoord[3]) < small &&
 977                         abs(thiscoord[1] - thiscoord[4]) < small &&
 978                         abs(thiscoord[2] - thiscoord[5]) < small) {
 979                     /* interchange first with second atom for better
 980                      * compression of water molecules
 981                      */
 982                     tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
 983                         thiscoord[3] = tmp;
 984                     tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
 985                         thiscoord[4] = tmp;
 986                     tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
 987                         thiscoord[5] = tmp;
 988                     is_small = 1;
 989                 }
 990
 991             }
 992             tmpcoord[0] = thiscoord[0] - minint[0];
 993             tmpcoord[1] = thiscoord[1] - minint[1];
 994             tmpcoord[2] = thiscoord[2] - minint[2];
 995             if (bitsize == 0) {
 996                 sendbits(buf, bitsizeint[0], tmpcoord[0]);
 997                 sendbits(buf, bitsizeint[1], tmpcoord[1]);
 998                 sendbits(buf, bitsizeint[2], tmpcoord[2]);
 999             } else {
1000                 sendints(buf, 3, bitsize, sizeint, tmpcoord);
1001             }
1002             prevcoord[0] = thiscoord[0];
1003             prevcoord[1] = thiscoord[1];
1004             prevcoord[2] = thiscoord[2];
1005             thiscoord = thiscoord + 3;
1006             i++;
1007
1008             run = 0;
1009             if (is_small == 0 && is_smaller == -1)
1010                 is_smaller = 0;
1011             while (is_small && run < 8*3) {
1012                 if (is_smaller == -1 && (
1013                         SQR(thiscoord[0] - prevcoord[0]) +
1014                         SQR(thiscoord[1] - prevcoord[1]) +
1015                         SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) {
1016                     is_smaller = 0;
1017                 }
1018
1019                 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small;
1020                 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small;
1021                 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small;
1022
1023                 prevcoord[0] = thiscoord[0];
1024                 prevcoord[1] = thiscoord[1];
1025                 prevcoord[2] = thiscoord[2];
1026
1027                 i++;
1028                 thiscoord = thiscoord + 3;
1029                 is_small = 0;
1030                 if (i < *size &&
1031                         abs(thiscoord[0] - prevcoord[0]) < small &&
1032                         abs(thiscoord[1] - prevcoord[1]) < small &&
1033                         abs(thiscoord[2] - prevcoord[2]) < small) {
1034                     is_small = 1;
1035                 }
1036             }
1037             if (run != prevrun || is_smaller != 0) {
1038                 prevrun = run;
1039                 sendbits(buf, 1, 1); /* flag the change in run-length */
1040                 sendbits(buf, 5, run+is_smaller+1);
1041             } else {
1042                 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
1043             }
1044             for (k=0; k < run; k+=3) {
1045                 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
1046             }
1047             if (is_smaller != 0) {
1048                 smallidx += is_smaller;
1049                 if (is_smaller < 0) {
1050                     small = smaller;
1051                     smaller = magicints[smallidx-1] / 2;
1052                 } else {
1053                     smaller = small;
1054                     small = magicints[smallidx] / 2;
1055                 }
1056                 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1057             }
1058         }
1059         if (buf[1] != 0) buf[0]++;;
1060         xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */
1061         return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]));
1062     } else {
1063
1064         /* xdrs is open for reading */
1065
1066         if (xdr_int(xdrs, &lsize) == 0)
1067             return 0;
1068         if (*size != 0 && lsize != *size) {
1069             fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; "
1070                     "%d arg vs %d in file", *size, lsize);
1071         }
1072         *size = lsize;
1073         size3 = *size * 3;
1074         if (*size <= 9) {
1075             return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
1076                 (xdrproc_t)xdr_float));
1077         }
1078         xdr_float(xdrs, precision);
1079         if (ip == NULL) {
1080             ip = (int *)malloc(size3 * sizeof(*ip));
1081             if (ip == NULL) {
1082                 fprintf(stderr,"malloc failed\n");
1083                 exit(1);
1084             }
1085             bufsize = size3 * 1.2;
1086             buf = (int *)malloc(bufsize * sizeof(*buf));
1087             if (buf == NULL) {
1088                 fprintf(stderr,"malloc failed\n");
1089                 exit(1);
1090             }
1091             oldsize = *size;
1092         } else if (*size > oldsize) {
1093             ip = (int *)realloc(ip, size3 * sizeof(*ip));
1094             if (ip == NULL) {
1095                 fprintf(stderr,"malloc failed\n");
1096                 exit(1);
1097             }
1098             bufsize = size3 * 1.2;
1099             buf = (int *)realloc(buf, bufsize * sizeof(*buf));
1100             if (buf == NULL) {
1101                 fprintf(stderr,"malloc failed\n");
1102                 exit(1);
1103             }
1104             oldsize = *size;
1105         }
1106         buf[0] = buf[1] = buf[2] = 0;
1107
1108         xdr_int(xdrs, &(minint[0]));
1109         xdr_int(xdrs, &(minint[1]));
1110         xdr_int(xdrs, &(minint[2]));
1111
1112         xdr_int(xdrs, &(maxint[0]));
1113         xdr_int(xdrs, &(maxint[1]));
1114         xdr_int(xdrs, &(maxint[2]));
1115
1116         sizeint[0] = maxint[0] - minint[0]+1;
1117         sizeint[1] = maxint[1] - minint[1]+1;
1118         sizeint[2] = maxint[2] - minint[2]+1;
1119
1120         /* check if one of the sizes is to big to be multiplied */
1121         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
1122             bitsizeint[0] = sizeofint(sizeint[0]);
1123             bitsizeint[1] = sizeofint(sizeint[1]);
1124             bitsizeint[2] = sizeofint(sizeint[2]);
1125             bitsize = 0; /* flag the use of large sizes */
1126         } else {
1127             bitsize = sizeofints(3, sizeint);
1128         }
1129
1130         xdr_int(xdrs, &smallidx);
1131         maxidx = MIN(LASTIDX, smallidx + 8) ;
1132         minidx = maxidx - 8; /* often this equal smallidx */
1133         smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
1134         small = magicints[smallidx] / 2;
1135         sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1136         larger = magicints[maxidx];
1137
1138         /* buf[0] holds the length in bytes */
1139
1140         if (xdr_int(xdrs, &(buf[0])) == 0)
1141             return 0;
1142         if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0)
1143             return 0;
1144         buf[0] = buf[1] = buf[2] = 0;
1145
1146         lfp = fp;
1147         inv_precision = 1.0 / * precision;
1148         run = 0;
1149         i = 0;
1150         lip = ip;
1151         while ( i < lsize ) {
1152             thiscoord = (int *)(lip) + i * 3;
1153
1154             if (bitsize == 0) {
1155                 thiscoord[0] = receivebits(buf, bitsizeint[0]);
1156                 thiscoord[1] = receivebits(buf, bitsizeint[1]);
1157                 thiscoord[2] = receivebits(buf, bitsizeint[2]);
1158             } else {
1159                 receiveints(buf, 3, bitsize, sizeint, thiscoord);
1160             }
1161
1162             i++;
1163             thiscoord[0] += minint[0];
1164             thiscoord[1] += minint[1];
1165             thiscoord[2] += minint[2];
1166
1167             prevcoord[0] = thiscoord[0];
1168             prevcoord[1] = thiscoord[1];
1169             prevcoord[2] = thiscoord[2];
1170
1171
1172             flag = receivebits(buf, 1);
1173             is_smaller = 0;
1174             if (flag == 1) {
1175                 run = receivebits(buf, 5);
1176                 is_smaller = run % 3;
1177                 run -= is_smaller;
1178                 is_smaller--;
1179             }
1180             if (run > 0) {
1181                 thiscoord += 3;
1182                 for (k = 0; k < run; k+=3) {
1183                     receiveints(buf, 3, smallidx, sizesmall, thiscoord);
1184                     i++;
1185                     thiscoord[0] += prevcoord[0] - small;
1186                     thiscoord[1] += prevcoord[1] - small;
1187                     thiscoord[2] += prevcoord[2] - small;
1188                     if (k == 0) {
1189                         /* interchange first with second atom for better
1190                          * compression of water molecules
1191                          */
1192                         tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1193                                 prevcoord[0] = tmp;
1194                         tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1195                                 prevcoord[1] = tmp;
1196                         tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1197                                 prevcoord[2] = tmp;
1198                         *lfp++ = prevcoord[0] * inv_precision;
1199                         *lfp++ = prevcoord[1] * inv_precision;
1200                         *lfp++ = prevcoord[2] * inv_precision;
1201                     } else {
1202                         prevcoord[0] = thiscoord[0];
1203                         prevcoord[1] = thiscoord[1];
1204                         prevcoord[2] = thiscoord[2];
1205                     }
1206                     *lfp++ = thiscoord[0] * inv_precision;
1207                     *lfp++ = thiscoord[1] * inv_precision;
1208                     *lfp++ = thiscoord[2] * inv_precision;
1209                 }
1210             } else {
1211                 *lfp++ = thiscoord[0] * inv_precision;
1212                 *lfp++ = thiscoord[1] * inv_precision;
1213                 *lfp++ = thiscoord[2] * inv_precision;
1214             }
1215             smallidx += is_smaller;
1216             if (is_smaller < 0) {
1217                 small = smaller;
1218                 if (smallidx > FIRSTIDX) {
1219                     smaller = magicints[smallidx - 1] /2;
1220                 } else {
1221                     smaller = 0;
1222                 }
1223             } else if (is_smaller > 0) {
1224                 smaller = small;
1225                 small = magicints[smallidx] / 2;
1226             }
1227             sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1228         }
1229     }
1230     return 1;
1231 }
1232
1233
1234