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