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