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