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