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