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