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