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