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