source/wham/src-NEWSC/xdrf/libxdrf.m4.org

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