source/wham/src-M/xdrf.org/libxdrf.c

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