Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel
3 : All rights reserved.
4 :
5 : Redistribution and use in source and binary forms, with or without
6 : modification, are permitted provided that the following conditions are met:
7 :
8 : 1. Redistributions of source code must retain the above copyright notice, this
9 : list of conditions and the following disclaimer.
10 :
11 : 2. Redistributions in binary form must reproduce the above copyright notice,
12 : this list of conditions and the following disclaimer in the documentation
13 : and/or other materials provided with the distribution.
14 :
15 : THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
16 : AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17 : IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 : DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
19 : FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20 : DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
21 : SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
22 : CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
23 : OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
24 : OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
26 : /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*-
27 : *
28 : * $Id$
29 : *
30 : /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*-
31 : *
32 : * $Id$
33 : *
34 : * Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel
35 : * All rights reserved.
36 : *
37 : * Redistribution and use in source and binary forms, with or without
38 : * modification, are permitted provided that the following conditions are met:
39 : *
40 : * 1. Redistributions of source code must retain the above copyright notice, this
41 : * list of conditions and the following disclaimer.
42 : *
43 : * 2. Redistributions in binary form must reproduce the above copyright notice,
44 : * this list of conditions and the following disclaimer in the documentation
45 : * and/or other materials provided with the distribution.
46 : *
47 : * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
48 : * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
49 : * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
50 : * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
51 : * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
52 : * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
53 : * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
54 : * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
55 : * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
56 : * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
57 : */
58 :
59 : /* Get HAVE_RPC_XDR_H, F77_FUNC from config.h if available */
60 : #ifdef HAVE_CONFIG_H
61 : #include <config.h>
62 : #endif
63 :
64 : #include <stdio.h>
65 : #include <stdlib.h>
66 : #include <string.h>
67 : #include <math.h>
68 : #include <limits.h>
69 :
70 : #define _FILE_OFFSET_BITS 64
71 :
72 : /* get fixed-width types if we are using ANSI C99 */
73 : #ifdef HAVE_STDINT_H
74 : # include <stdint.h>
75 : #elif (defined HAVE_INTTYPES_H)
76 : # include <inttypes.h>
77 : #endif
78 :
79 :
80 : #ifdef HAVE_RPC_XDR_H
81 : # include <rpc/rpc.h>
82 : # include <rpc/xdr.h>
83 : #endif
84 :
85 : #include "xdrfile.h"
86 :
87 : /* Default FORTRAN name mangling is: lower case name, append underscore */
88 : #ifndef F77_FUNC
89 : #define F77_FUNC(name,NAME) name ## _
90 : #endif
91 :
92 : namespace PLMD
93 : {
94 : namespace xdrfile
95 : {
96 :
97 : const char *exdr_message[exdrNR] = {
98 : "OK",
99 : "Header",
100 : "String",
101 : "Double",
102 : "Integer",
103 : "Float",
104 : "Unsigned integer",
105 : "Compressed 3D coordinate",
106 : "Closing file",
107 : "Magic number",
108 : "Not enough memory",
109 : "End of file",
110 : "File not found"
111 : };
112 :
113 : /*
114 : * Declare our own XDR routines statically if no libraries are present.
115 : * Actual implementation is at the end of this file.
116 : *
117 : * We don't want the low-level XDR implementation as part of the Gromacs
118 : * documentation, so skip it for doxygen too...
119 : */
120 : #if (!defined HAVE_RPC_XDR_H && !defined DOXYGEN)
121 :
122 : enum xdr_op
123 : {
124 : XDR_ENCODE = 0,
125 : XDR_DECODE = 1,
126 : XDR_FREE = 2
127 : };
128 :
129 :
130 : /* We need integer types that are guaranteed to be 4 bytes wide.
131 : * If ANSI C99 headers were included they are already defined
132 : * as int32_t and uint32_t. Check, and if not define them ourselves.
133 : * Since it is just our workaround for missing ANSI C99 types, avoid adding
134 : * it to the doxygen documentation.
135 : */
136 : #if !(defined INT32_MAX || defined DOXYGEN)
137 : # if (INT_MAX == 2147483647)
138 : # define int32_t int
139 : # define uint32_t unsigned int
140 : # define INT32_MAX 2147483647
141 : # elif (LONG_MAX == 2147483647)
142 : # define int32_t long
143 : # define uint32_t unsigned long
144 : # define INT32_MAX 2147483647L
145 : # else
146 : # error ERROR: No 32 bit wide integer type found!
147 : # error Use system XDR libraries instead, or update xdrfile.c
148 : # endif
149 : #endif
150 :
151 : typedef struct XDR XDR;
152 :
153 : struct XDR
154 : {
155 : enum xdr_op x_op;
156 : struct xdr_ops
157 : {
158 : int (*x_getlong) (XDR *__xdrs, int32_t *__lp);
159 : int (*x_putlong) (XDR *__xdrs, int32_t *__lp);
160 : int (*x_getbytes) (XDR *__xdrs, char *__addr, unsigned int __len);
161 : int (*x_putbytes) (XDR *__xdrs, char *__addr, unsigned int __len);
162 : /* two next routines are not 64-bit IO safe - don't use! */
163 : unsigned int (*x_getpostn) (XDR *__xdrs);
164 : int (*x_setpostn) (XDR *__xdrs, unsigned int __pos);
165 : void (*x_destroy) (XDR *__xdrs);
166 : }
167 : *x_ops;
168 : char *x_private;
169 : };
170 :
171 : static int xdr_char (XDR *xdrs, char *ip);
172 : static int xdr_u_char (XDR *xdrs, unsigned char *ip);
173 : static int xdr_short (XDR *xdrs, short *ip);
174 : static int xdr_u_short (XDR *xdrs, unsigned short *ip);
175 : static int xdr_int (XDR *xdrs, int *ip);
176 : static int xdr_u_int (XDR *xdrs, unsigned int *ip);
177 : static int xdr_float (XDR *xdrs, float *ip);
178 : static int xdr_double (XDR *xdrs, double *ip);
179 : static int xdr_string (XDR *xdrs, char **ip, unsigned int maxsize);
180 : static int xdr_opaque (XDR *xdrs, char *cp, unsigned int cnt);
181 : static void xdrstdio_create (XDR *xdrs, FILE *fp, enum xdr_op xop);
182 :
183 : #define xdr_getpos(xdrs) \
184 : (*(xdrs)->x_ops->x_getpostn)(xdrs)
185 : #define xdr_setpos(xdrs, pos) \
186 : (*(xdrs)->x_ops->x_setpostn)(xdrs, pos)
187 : #define xdr_destroy(xdrs) \
188 : do { \
189 : if ((xdrs)->x_ops->x_destroy) \
190 : (*(xdrs)->x_ops->x_destroy)(xdrs); \
191 : } while (0)
192 : #endif /* end of our own XDR declarations */
193 :
194 :
195 :
196 :
197 :
198 : /** Contents of the abstract XDRFILE data structure.
199 : *
200 : * @internal
201 : *
202 : * This structure is used to provide an XDR file interface that is
203 : * virtual identical to the standard UNIX fopen/fread/fwrite/fclose.
204 : */
205 : struct XDRFILE
206 : {
207 : FILE * fp; /**< pointer to standard C library file handle */
208 : XDR * xdr; /**< pointer to corresponding XDR handle */
209 : char mode; /**< r=read, w=write, a=append */
210 : int * buf1; /**< Buffer for internal use */
211 : int buf1size; /**< Current allocated length of buf1 */
212 : int * buf2; /**< Buffer for internal use */
213 : int buf2size; /**< Current allocated length of buf2 */
214 : };
215 :
216 :
217 :
218 :
219 : /*************************************************************
220 : * Implementation of higher-level routines to read/write *
221 : * portable data based on the XDR standard. These should be *
222 : * called from C - see further down for Fortran77 wrappers. *
223 : *************************************************************/
224 :
225 : XDRFILE *
226 32 : xdrfile_open(const char *path, const char *mode)
227 : {
228 : char newmode[5];
229 : enum xdr_op xdrmode;
230 : XDRFILE *xfp;
231 :
232 : /* make sure XDR files are opened in binary mode... */
233 32 : if(*mode=='w' || *mode=='W')
234 : {
235 : snprintf(newmode,5,"wb+");
236 : xdrmode=XDR_ENCODE;
237 22 : } else if(*mode == 'a' || *mode == 'A')
238 : {
239 : snprintf(newmode,5,"ab+");
240 : xdrmode = XDR_ENCODE;
241 22 : } else if(*mode == 'r' || *mode == 'R')
242 : {
243 : snprintf(newmode,5,"rb");
244 : xdrmode = XDR_DECODE;
245 : } else /* cannot determine mode */
246 : return NULL;
247 :
248 32 : if((xfp=(XDRFILE *)malloc(sizeof(XDRFILE)))==NULL)
249 : return NULL;
250 32 : if((xfp->fp=fopen(path,newmode))==NULL)
251 : {
252 0 : free(xfp);
253 0 : return NULL;
254 : }
255 32 : if((xfp->xdr=(XDR *)malloc(sizeof(XDR)))==NULL)
256 : {
257 0 : fclose(xfp->fp);
258 0 : free(xfp);
259 0 : return NULL;
260 : }
261 32 : xfp->mode=*mode;
262 32 : xdrstdio_create((XDR *)(xfp->xdr),xfp->fp,xdrmode);
263 32 : xfp->buf1 = xfp->buf2 = NULL;
264 32 : xfp->buf1size = xfp->buf2size = 0;
265 32 : return xfp;
266 : }
267 :
268 : int
269 32 : xdrfile_close(XDRFILE *xfp)
270 : {
271 : int ret=exdrCLOSE;
272 32 : if(xfp)
273 : {
274 : /* flush and destroy XDR stream */
275 32 : if(xfp->xdr)
276 32 : xdr_destroy((XDR *)(xfp->xdr));
277 32 : free(xfp->xdr);
278 : /* close the file */
279 32 : ret=fclose(xfp->fp);
280 32 : if(xfp->buf1size)
281 7 : free(xfp->buf1);
282 32 : if(xfp->buf2size)
283 7 : free(xfp->buf2);
284 32 : free(xfp);
285 : }
286 32 : return ret; /* return 0 if ok */
287 : }
288 :
289 :
290 :
291 : int
292 1994 : xdrfile_read_int(int *ptr, int ndata, XDRFILE* xfp)
293 : {
294 : int i=0;
295 :
296 : /* read write is encoded in the XDR struct */
297 4091 : while(i<ndata && xdr_int((XDR *)(xfp->xdr),ptr+i))
298 2097 : i++;
299 :
300 1994 : return i;
301 : }
302 :
303 : int
304 620 : xdrfile_write_int(int *ptr, int ndata, XDRFILE* xfp)
305 : {
306 : int i=0;
307 :
308 : /* read write is encoded in the XDR struct */
309 1478 : while(i<ndata && xdr_int((XDR *)(xfp->xdr),ptr+i))
310 858 : i++;
311 620 : return i;
312 : }
313 :
314 :
315 : int
316 0 : xdrfile_read_uint(unsigned int *ptr, int ndata, XDRFILE* xfp)
317 : {
318 : int i=0;
319 :
320 : /* read write is encoded in the XDR struct */
321 0 : while(i<ndata && xdr_u_int((XDR *)(xfp->xdr),ptr+i))
322 0 : i++;
323 :
324 0 : return i;
325 : }
326 :
327 : int
328 0 : xdrfile_write_uint(unsigned int *ptr, int ndata, XDRFILE* xfp)
329 : {
330 : int i=0;
331 :
332 : /* read write is encoded in the XDR struct */
333 0 : while(i<ndata && xdr_u_int((XDR *)(xfp->xdr),ptr+i))
334 0 : i++;
335 0 : return i;
336 : }
337 :
338 : int
339 0 : xdrfile_read_char(char *ptr, int ndata, XDRFILE* xfp)
340 : {
341 : int i=0;
342 :
343 : /* read write is encoded in the XDR struct */
344 0 : while(i<ndata && xdr_char((XDR *)(xfp->xdr),ptr+i))
345 0 : i++;
346 :
347 0 : return i;
348 : }
349 :
350 : int
351 0 : xdrfile_write_char(char *ptr, int ndata, XDRFILE* xfp)
352 : {
353 : int i=0;
354 :
355 : /* read write is encoded in the XDR struct */
356 0 : while(i<ndata && xdr_char((XDR *)(xfp->xdr),ptr+i))
357 0 : i++;
358 0 : return i;
359 : }
360 :
361 :
362 : int
363 0 : xdrfile_read_uchar(unsigned char *ptr, int ndata, XDRFILE* xfp)
364 : {
365 : int i=0;
366 :
367 : /* read write is encoded in the XDR struct */
368 0 : while(i<ndata && xdr_u_char((XDR *)(xfp->xdr),ptr+i))
369 0 : i++;
370 :
371 0 : return i;
372 : }
373 :
374 : int
375 0 : xdrfile_write_uchar(unsigned char *ptr, int ndata, XDRFILE* xfp)
376 : {
377 : int i=0;
378 :
379 : /* read write is encoded in the XDR struct */
380 0 : while(i<ndata && xdr_u_char((XDR *)(xfp->xdr),ptr+i))
381 0 : i++;
382 0 : return i;
383 : }
384 :
385 : int
386 0 : xdrfile_read_short(short *ptr, int ndata, XDRFILE* xfp)
387 : {
388 : int i=0;
389 :
390 : /* read write is encoded in the XDR struct */
391 0 : while(i<ndata && xdr_short((XDR *)(xfp->xdr),ptr+i))
392 0 : i++;
393 :
394 0 : return i;
395 : }
396 :
397 : int
398 0 : xdrfile_write_short(short *ptr, int ndata, XDRFILE* xfp)
399 : {
400 : int i=0;
401 :
402 : /* read write is encoded in the XDR struct */
403 0 : while(i<ndata && xdr_short((XDR *)(xfp->xdr),ptr+i))
404 0 : i++;
405 0 : return i;
406 : }
407 :
408 :
409 : int
410 0 : xdrfile_read_ushort(unsigned short *ptr, int ndata, XDRFILE* xfp)
411 : {
412 : int i=0;
413 :
414 : /* read write is encoded in the XDR struct */
415 0 : while(i<ndata && xdr_u_short((XDR *)(xfp->xdr),ptr+i))
416 0 : i++;
417 :
418 0 : return i;
419 : }
420 :
421 : int
422 0 : xdrfile_write_ushort(unsigned short *ptr, int ndata, XDRFILE* xfp)
423 : {
424 : int i=0;
425 :
426 : /* read write is encoded in the XDR struct */
427 0 : while(i<ndata && xdr_u_short((XDR *)(xfp->xdr),ptr+i))
428 0 : i++;
429 0 : return i;
430 : }
431 :
432 : int
433 602 : xdrfile_read_float(float *ptr, int ndata, XDRFILE* xfp)
434 : {
435 : int i=0;
436 : /* read write is encoded in the XDR struct */
437 23916 : while(i<ndata && xdr_float((XDR *)(xfp->xdr),ptr+i))
438 23314 : i++;
439 602 : return i;
440 : }
441 :
442 : int
443 174 : xdrfile_write_float(float *ptr, int ndata, XDRFILE* xfp)
444 : {
445 : int i=0;
446 : /* read write is encoded in the XDR struct */
447 336 : while(i<ndata && xdr_float((XDR *)(xfp->xdr),ptr+i))
448 162 : i++;
449 174 : return i;
450 : }
451 :
452 : int
453 0 : xdrfile_read_double(double *ptr, int ndata, XDRFILE* xfp)
454 : {
455 : int i=0;
456 : /* read write is encoded in the XDR struct */
457 0 : while(i<ndata && xdr_double((XDR *)(xfp->xdr),ptr+i))
458 0 : i++;
459 0 : return i;
460 : }
461 :
462 : int
463 0 : xdrfile_write_double(double *ptr, int ndata, XDRFILE* xfp)
464 : {
465 : int i=0;
466 : /* read write is encoded in the XDR struct */
467 0 : while(i<ndata && xdr_double((XDR *)(xfp->xdr),ptr+i))
468 0 : i++;
469 0 : return i;
470 : }
471 :
472 : int
473 75 : xdrfile_read_string(char *ptr, int maxlen, XDRFILE* xfp)
474 : {
475 : int i;
476 75 : if(xdr_string((XDR *)(xfp->xdr),&ptr,maxlen)) {
477 : i=0;
478 975 : while(i<maxlen && ptr[i]!=0)
479 900 : i++;
480 75 : if(i==maxlen)
481 : return maxlen;
482 : else
483 75 : return i+1;
484 : } else
485 : return 0;
486 : }
487 :
488 : int
489 48 : xdrfile_write_string(char *ptr, XDRFILE* xfp)
490 : {
491 48 : int len=strlen(ptr)+1;
492 :
493 48 : if(xdr_string((XDR *)(xfp->xdr),&ptr,len))
494 : return len;
495 : else
496 0 : return 0;
497 : }
498 :
499 :
500 : int
501 28 : xdrfile_read_opaque(char *ptr, int cnt, XDRFILE* xfp)
502 : {
503 28 : if(xdr_opaque((XDR *)(xfp->xdr),ptr,cnt))
504 : return cnt;
505 : else
506 0 : return 0;
507 : }
508 :
509 :
510 : int
511 60 : xdrfile_write_opaque(char *ptr, int cnt, XDRFILE* xfp)
512 : {
513 60 : if(xdr_opaque((XDR *)(xfp->xdr),ptr,cnt))
514 : return cnt;
515 : else
516 0 : return 0;
517 : }
518 :
519 :
520 : /* Internal support routines for reading/writing compressed coordinates
521 : * sizeofint - calculate smallest number of bits necessary
522 : * to represent a certain integer.
523 : */
524 : static int
525 0 : sizeofint(int size) {
526 : unsigned int num = 1;
527 : int num_of_bits = 0;
528 :
529 0 : while (size >= num && num_of_bits < 32)
530 : {
531 0 : num_of_bits++;
532 0 : num <<= 1;
533 : }
534 0 : return num_of_bits;
535 : }
536 :
537 :
538 : /*
539 : * sizeofints - calculate 'bitsize' of compressed ints
540 : *
541 : * given a number of small unsigned integers and the maximum value
542 : * return the number of bits needed to read or write them with the
543 : * routines encodeints/decodeints. You need this parameter when
544 : * calling those routines.
545 : * (However, in some cases we can just use the variable 'smallidx'
546 : * which is the exact number of bits, and them we dont need to call
547 : * this routine).
548 : */
549 : static int
550 88 : sizeofints(int num_of_ints, unsigned int sizes[])
551 : {
552 : int i, num;
553 : unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
554 : num_of_bytes = 1;
555 88 : bytes[0] = 1;
556 : num_of_bits = 0;
557 352 : for (i=0; i < num_of_ints; i++)
558 : {
559 : tmp = 0;
560 876 : for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
561 : {
562 612 : tmp = bytes[bytecnt] * sizes[i] + tmp;
563 612 : bytes[bytecnt] = tmp & 0xff;
564 612 : tmp >>= 8;
565 : }
566 612 : while (tmp != 0)
567 : {
568 348 : bytes[bytecnt++] = tmp & 0xff;
569 348 : tmp >>= 8;
570 : }
571 : num_of_bytes = bytecnt;
572 : }
573 : num = 1;
574 88 : num_of_bytes--;
575 497 : while (bytes[num_of_bytes] >= num)
576 : {
577 409 : num_of_bits++;
578 409 : num *= 2;
579 : }
580 88 : return num_of_bits + num_of_bytes * 8;
581 :
582 : }
583 :
584 :
585 : /*
586 : * encodebits - encode num into buf using the specified number of bits
587 : *
588 : * This routines appends the value of num to the bits already present in
589 : * the array buf. You need to give it the number of bits to use and you had
590 : * better make sure that this number of bits is enough to hold the value.
591 : * Num must also be positive.
592 : */
593 : static void
594 21991 : encodebits(int buf[], int num_of_bits, int num)
595 : {
596 :
597 : unsigned int cnt, lastbyte;
598 : int lastbits;
599 : unsigned char * cbuf;
600 :
601 21991 : cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
602 21991 : cnt = (unsigned int) buf[0];
603 21991 : lastbits = buf[1];
604 21991 : lastbyte =(unsigned int) buf[2];
605 38553 : while (num_of_bits >= 8)
606 : {
607 16562 : lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
608 16562 : cbuf[cnt++] = lastbyte >> lastbits;
609 : num_of_bits -= 8;
610 : }
611 21991 : if (num_of_bits > 0)
612 : {
613 4729 : lastbyte = (lastbyte << num_of_bits) | num;
614 4729 : lastbits += num_of_bits;
615 4729 : if (lastbits >= 8)
616 : {
617 1916 : lastbits -= 8;
618 1916 : cbuf[cnt++] = lastbyte >> lastbits;
619 : }
620 : }
621 21991 : buf[0] = cnt;
622 21991 : buf[1] = lastbits;
623 21991 : buf[2] = lastbyte;
624 21991 : if (lastbits>0)
625 : {
626 19395 : cbuf[cnt] = lastbyte << (8 - lastbits);
627 : }
628 21991 : }
629 :
630 : /*
631 : * encodeints - encode a small set of small integers in compressed format
632 : *
633 : * this routine is used internally by xdr3dfcoord, to encode a set of
634 : * small integers to the buffer for writing to a file.
635 : * Multiplication with fixed (specified maximum) sizes is used to get
636 : * to one big, multibyte integer. Allthough the routine could be
637 : * modified to handle sizes bigger than 16777216, or more than just
638 : * a few integers, this is not done because the gain in compression
639 : * isn't worth the effort. Note that overflowing the multiplication
640 : * or the byte buffer (32 bytes) is unchecked and whould cause bad results.
641 : * THese things are checked in the calling routines, so make sure not
642 : * to remove those checks...
643 : */
644 :
645 : static void
646 3960 : encodeints(int buf[], int num_of_ints, int num_of_bits,
647 : unsigned int sizes[], unsigned int nums[])
648 : {
649 :
650 : int i;
651 : unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
652 :
653 3960 : tmp = nums[0];
654 : num_of_bytes = 0;
655 : do
656 : {
657 7242 : bytes[num_of_bytes++] = tmp & 0xff;
658 7242 : tmp >>= 8;
659 7242 : } while (tmp != 0);
660 :
661 11880 : for (i = 1; i < num_of_ints; i++)
662 : {
663 7920 : if (nums[i] >= sizes[i])
664 : {
665 0 : fprintf(stderr,"major breakdown in encodeints - num %u doesn't "
666 : "match size %u\n", nums[i], sizes[i]);
667 0 : abort();
668 : }
669 : /* use one step multiply */
670 : tmp = nums[i];
671 28600 : for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
672 : {
673 20680 : tmp = bytes[bytecnt] * sizes[i] + tmp;
674 20680 : bytes[bytecnt] = tmp & 0xff;
675 20680 : tmp >>= 8;
676 : }
677 20013 : while (tmp != 0)
678 : {
679 12093 : bytes[bytecnt++] = tmp & 0xff;
680 12093 : tmp >>= 8;
681 : }
682 : num_of_bytes = bytecnt;
683 : }
684 3960 : if (num_of_bits >= num_of_bytes * 8)
685 : {
686 7225 : for (i = 0; i < num_of_bytes; i++)
687 : {
688 6040 : encodebits(buf, 8, bytes[i]);
689 : }
690 1185 : encodebits(buf, num_of_bits - num_of_bytes * 8, 0);
691 : }
692 : else
693 : {
694 13295 : for (i = 0; i < num_of_bytes-1; i++)
695 : {
696 10520 : encodebits(buf, 8, bytes[i]);
697 : }
698 2775 : encodebits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
699 : }
700 3960 : }
701 :
702 :
703 : /*
704 : * decodebits - decode number from buf using specified number of bits
705 : *
706 : * extract the number of bits from the array buf and construct an integer
707 : * from it. Return that value.
708 : *
709 : */
710 :
711 : static int
712 14579 : decodebits(int buf[], int num_of_bits)
713 : {
714 :
715 : int cnt, num;
716 : unsigned int lastbits, lastbyte;
717 : unsigned char * cbuf;
718 14579 : int mask = (1 << num_of_bits) -1;
719 :
720 14579 : cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
721 14579 : cnt = buf[0];
722 14579 : lastbits = (unsigned int) buf[1];
723 14579 : lastbyte = (unsigned int) buf[2];
724 :
725 : num = 0;
726 26636 : while (num_of_bits >= 8)
727 : {
728 12057 : lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
729 12057 : num |= (lastbyte >> lastbits) << (num_of_bits - 8);
730 : num_of_bits -=8;
731 : }
732 14579 : if (num_of_bits > 0)
733 : {
734 2522 : if (lastbits < num_of_bits)
735 : {
736 998 : lastbits += 8;
737 998 : lastbyte = (lastbyte << 8) | cbuf[cnt++];
738 : }
739 2522 : lastbits -= num_of_bits;
740 2522 : num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
741 : }
742 14579 : num &= mask;
743 14579 : buf[0] = cnt;
744 14579 : buf[1] = lastbits;
745 14579 : buf[2] = lastbyte;
746 14579 : return num;
747 : }
748 :
749 : /*
750 : * decodeints - decode 'small' integers from the buf array
751 : *
752 : * this routine is the inverse from encodeints() and decodes the small integers
753 : * written to buf by calculating the remainder and doing divisions with
754 : * the given sizes[]. You need to specify the total number of bits to be
755 : * used from buf in num_of_bits.
756 : *
757 : */
758 :
759 : static void
760 1936 : decodeints(int buf[], int num_of_ints, int num_of_bits,
761 : unsigned int sizes[], int nums[])
762 : {
763 :
764 : int bytes[32];
765 : int i, j, num_of_bytes, p, num;
766 :
767 1936 : bytes[1] = bytes[2] = bytes[3] = 0;
768 : num_of_bytes = 0;
769 13916 : while (num_of_bits > 8)
770 : {
771 11980 : bytes[num_of_bytes++] = decodebits(buf, 8);
772 11980 : num_of_bits -= 8;
773 : }
774 1936 : if (num_of_bits > 0)
775 : {
776 1936 : bytes[num_of_bytes++] = decodebits(buf, num_of_bits);
777 : }
778 5808 : for (i = num_of_ints-1; i > 0; i--)
779 : {
780 : num = 0;
781 31704 : for (j = num_of_bytes-1; j >=0; j--)
782 : {
783 27832 : num = (num << 8) | bytes[j];
784 27832 : p = num / sizes[i];
785 27832 : bytes[j] = p;
786 27832 : num = num - p * sizes[i];
787 : }
788 3872 : nums[i] = num;
789 : }
790 1936 : nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
791 1936 : }
792 :
793 :
794 : static const int magicints[] =
795 : {
796 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
797 : 80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 812, 1024, 1290,
798 : 1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192, 10321, 13003,
799 : 16384, 20642, 26007, 32768, 41285, 52015, 65536,82570, 104031,
800 : 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
801 : 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021,
802 : 4194304, 5284491, 6658042, 8388607, 10568983, 13316085, 16777216
803 : };
804 :
805 : #define FIRSTIDX 9
806 : /* note that magicints[FIRSTIDX-1] == 0 */
807 : #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
808 :
809 : /* Compressed coordinate routines - modified from the original
810 : * implementation by Frans v. Hoesel to make them threadsafe.
811 : */
812 : int
813 28 : xdrfile_decompress_coord_float(float *ptr,
814 : int *size,
815 : float *precision,
816 : XDRFILE* xfp)
817 : {
818 : int minint[3], maxint[3], *lip;
819 : int smallidx, minidx, maxidx;
820 : unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3;
821 : int k, *buf1, *buf2, lsize, flag;
822 : int smallnum, smaller, larger, i, is_smaller, run;
823 : float *lfp, inv_precision;
824 : int tmp, *thiscoord, prevcoord[3];
825 : unsigned int bitsize;
826 :
827 : bitsizeint[0] = 0;
828 : bitsizeint[1] = 0;
829 : bitsizeint[2] = 0;
830 :
831 28 : if(xfp==NULL || ptr==NULL)
832 : return -1;
833 28 : tmp=xdrfile_read_int(&lsize,1,xfp);
834 28 : if(tmp==0)
835 : return -1; /* return if we could not read size */
836 28 : if (*size < lsize)
837 : {
838 0 : fprintf(stderr, "Requested to decompress %d coords, file contains %d\n",
839 : *size, lsize);
840 0 : return -1;
841 : }
842 28 : *size = lsize;
843 28 : size3 = *size * 3;
844 28 : if(size3>xfp->buf1size)
845 : {
846 2 : if((xfp->buf1=(int *)malloc(sizeof(int)*size3))==NULL)
847 : {
848 0 : fprintf(stderr,"Cannot allocate memory for decompressing coordinates.\n");
849 0 : return -1;
850 : }
851 2 : xfp->buf1size=size3;
852 2 : xfp->buf2size=size3*1.2;
853 2 : if((xfp->buf2=(int *)malloc(sizeof(int)*xfp->buf2size))==NULL)
854 : {
855 0 : fprintf(stderr,"Cannot allocate memory for decompressing coordinates.\n");
856 0 : return -1;
857 : }
858 : }
859 : /* Dont bother with compression for three atoms or less */
860 28 : if(*size<=9)
861 : {
862 0 : return xdrfile_read_float(ptr,size3,xfp)/3;
863 : /* return number of coords, not floats */
864 : }
865 : /* Compression-time if we got here. Read precision first */
866 28 : xdrfile_read_float(precision,1,xfp);
867 :
868 : /* avoid repeated pointer dereferencing. */
869 28 : buf1=xfp->buf1;
870 28 : buf2=xfp->buf2;
871 : /* buf2[0-2] are special and do not contain actual data */
872 28 : buf2[0] = buf2[1] = buf2[2] = 0;
873 28 : xdrfile_read_int(minint,3,xfp);
874 28 : xdrfile_read_int(maxint,3,xfp);
875 :
876 28 : sizeint[0] = maxint[0] - minint[0]+1;
877 28 : sizeint[1] = maxint[1] - minint[1]+1;
878 28 : sizeint[2] = maxint[2] - minint[2]+1;
879 :
880 : /* check if one of the sizes is to big to be multiplied */
881 28 : if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
882 : {
883 0 : bitsizeint[0] = sizeofint(sizeint[0]);
884 0 : bitsizeint[1] = sizeofint(sizeint[1]);
885 0 : bitsizeint[2] = sizeofint(sizeint[2]);
886 : bitsize = 0; /* flag the use of large sizes */
887 : }
888 : else
889 : {
890 28 : bitsize = sizeofints(3, sizeint);
891 : }
892 :
893 28 : if (xdrfile_read_int(&smallidx,1,xfp) == 0)
894 : return 0; /* not sure what has happened here or why we return... */
895 28 : tmp=smallidx+8;
896 : maxidx = (LASTIDX<tmp) ? LASTIDX : tmp;
897 : minidx = maxidx - 8; /* often this equal smallidx */
898 28 : tmp = smallidx-1;
899 28 : tmp = (FIRSTIDX>tmp) ? FIRSTIDX : tmp;
900 28 : smaller = magicints[tmp] / 2;
901 28 : smallnum = magicints[smallidx] / 2;
902 28 : sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
903 : larger = magicints[maxidx];
904 :
905 : /* buf2[0] holds the length in bytes */
906 :
907 28 : if (xdrfile_read_int(buf2,1,xfp) == 0)
908 : return 0;
909 28 : if (xdrfile_read_opaque((char *)&(buf2[3]),(unsigned int)buf2[0],xfp) == 0)
910 : return 0;
911 28 : buf2[0] = buf2[1] = buf2[2] = 0;
912 :
913 : lfp = ptr;
914 28 : inv_precision = 1.0 / * precision;
915 : run = 0;
916 : i = 0;
917 : lip = buf1;
918 387 : while ( i < lsize )
919 : {
920 359 : thiscoord = (int *)(lip) + i * 3;
921 :
922 359 : if (bitsize == 0)
923 : {
924 0 : thiscoord[0] = decodebits(buf2, bitsizeint[0]);
925 0 : thiscoord[1] = decodebits(buf2, bitsizeint[1]);
926 0 : thiscoord[2] = decodebits(buf2, bitsizeint[2]);
927 : }
928 : else
929 : {
930 359 : decodeints(buf2, 3, bitsize, sizeint, thiscoord);
931 : }
932 :
933 359 : i++;
934 359 : thiscoord[0] += minint[0];
935 359 : thiscoord[1] += minint[1];
936 359 : thiscoord[2] += minint[2];
937 :
938 : prevcoord[0] = thiscoord[0];
939 : prevcoord[1] = thiscoord[1];
940 : prevcoord[2] = thiscoord[2];
941 :
942 359 : flag = decodebits(buf2, 1);
943 : is_smaller = 0;
944 359 : if (flag == 1)
945 : {
946 304 : run = decodebits(buf2, 5);
947 304 : is_smaller = run % 3;
948 304 : run -= is_smaller;
949 304 : is_smaller--;
950 : }
951 359 : if (run > 0)
952 : {
953 241 : thiscoord += 3;
954 1818 : for (k = 0; k < run; k+=3)
955 : {
956 1577 : decodeints(buf2, 3, smallidx, sizesmall, thiscoord);
957 1577 : i++;
958 1577 : thiscoord[0] += prevcoord[0] - smallnum;
959 1577 : thiscoord[1] += prevcoord[1] - smallnum;
960 1577 : thiscoord[2] += prevcoord[2] - smallnum;
961 1577 : if (k == 0) {
962 : /* interchange first with second atom for better
963 : * compression of water molecules
964 : */
965 241 : tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
966 : prevcoord[0] = tmp;
967 241 : tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
968 : prevcoord[1] = tmp;
969 241 : tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
970 : prevcoord[2] = tmp;
971 241 : *lfp++ = prevcoord[0] * inv_precision;
972 241 : *lfp++ = prevcoord[1] * inv_precision;
973 241 : *lfp++ = prevcoord[2] * inv_precision;
974 : } else {
975 : prevcoord[0] = thiscoord[0];
976 : prevcoord[1] = thiscoord[1];
977 : prevcoord[2] = thiscoord[2];
978 : }
979 1577 : *lfp++ = thiscoord[0] * inv_precision;
980 1577 : *lfp++ = thiscoord[1] * inv_precision;
981 1577 : *lfp++ = thiscoord[2] * inv_precision;
982 : }
983 : }
984 : else
985 : {
986 118 : *lfp++ = thiscoord[0] * inv_precision;
987 118 : *lfp++ = thiscoord[1] * inv_precision;
988 118 : *lfp++ = thiscoord[2] * inv_precision;
989 : }
990 359 : smallidx += is_smaller;
991 359 : if (is_smaller < 0)
992 : {
993 : smallnum = smaller;
994 :
995 38 : if (smallidx > FIRSTIDX)
996 : {
997 38 : smaller = magicints[smallidx - 1] /2;
998 : }
999 : else
1000 : {
1001 : smaller = 0;
1002 : }
1003 : }
1004 321 : else if (is_smaller > 0)
1005 : {
1006 : smaller = smallnum;
1007 230 : smallnum = magicints[smallidx] / 2;
1008 : }
1009 359 : sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1010 : }
1011 28 : return *size;
1012 : }
1013 :
1014 : int
1015 72 : xdrfile_compress_coord_float(float *ptr,
1016 : int size,
1017 : float precision,
1018 : XDRFILE* xfp)
1019 : {
1020 : int minint[3], maxint[3], mindiff, *lip, diff;
1021 : int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
1022 : int minidx, maxidx;
1023 : unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
1024 : int k, *buf1, *buf2;
1025 : int smallnum, smaller, larger, i, j, is_small, is_smaller, run, prevrun;
1026 : float *lfp, lf;
1027 : int tmp, tmpsum, *thiscoord, prevcoord[3];
1028 : unsigned int tmpcoord[30];
1029 : int errval=1;
1030 : unsigned int bitsize;
1031 :
1032 72 : if(xfp==NULL)
1033 : return -1;
1034 72 : size3=3*size;
1035 :
1036 : bitsizeint[0] = 0;
1037 : bitsizeint[1] = 0;
1038 : bitsizeint[2] = 0;
1039 :
1040 72 : if(size3>xfp->buf1size)
1041 : {
1042 5 : if((xfp->buf1=(int *)malloc(sizeof(int)*size3))==NULL)
1043 : {
1044 0 : fprintf(stderr,"Cannot allocate memory for compressing coordinates.\n");
1045 0 : return -1;
1046 : }
1047 5 : xfp->buf1size=size3;
1048 5 : xfp->buf2size=size3*1.2;
1049 5 : if((xfp->buf2=(int *)malloc(sizeof(int)*xfp->buf2size))==NULL)
1050 : {
1051 0 : fprintf(stderr,"Cannot allocate memory for compressing coordinates.\n");
1052 0 : return -1;
1053 : }
1054 : }
1055 72 : if(xdrfile_write_int(&size,1,xfp)==0)
1056 : return -1; /* return if we could not write size */
1057 : /* Dont bother with compression for three atoms or less */
1058 72 : if(size<=9)
1059 : {
1060 12 : return xdrfile_write_float(ptr,size3,xfp)/3;
1061 : /* return number of coords, not floats */
1062 : }
1063 : /* Compression-time if we got here. Write precision first */
1064 60 : if (precision <= 0)
1065 0 : precision = 1000;
1066 60 : xdrfile_write_float(&precision,1,xfp);
1067 : /* avoid repeated pointer dereferencing. */
1068 60 : buf1=xfp->buf1;
1069 60 : buf2=xfp->buf2;
1070 : /* buf2[0-2] are special and do not contain actual data */
1071 60 : buf2[0] = buf2[1] = buf2[2] = 0;
1072 60 : minint[0] = minint[1] = minint[2] = INT_MAX;
1073 60 : maxint[0] = maxint[1] = maxint[2] = INT_MIN;
1074 : prevrun = -1;
1075 : lfp = ptr;
1076 : lip = buf1;
1077 : mindiff = INT_MAX;
1078 : oldlint1 = oldlint2 = oldlint3 = 0;
1079 4020 : while(lfp < ptr + size3 )
1080 : {
1081 : /* find nearest integer */
1082 3960 : if (*lfp >= 0.0)
1083 35 : lf = *lfp * precision + 0.5;
1084 : else
1085 3925 : lf = *lfp * precision - 0.5;
1086 3960 : if (fabs(lf) > INT_MAX-2)
1087 : {
1088 : /* scaling would cause overflow */
1089 0 : fprintf(stderr,"Internal overflow compressing coordinates.\n");
1090 : errval=0;
1091 : }
1092 3960 : lint1 = lf;
1093 3960 : if (lint1 < minint[0]) minint[0] = lint1;
1094 3960 : if (lint1 > maxint[0]) maxint[0] = lint1;
1095 3960 : *lip++ = lint1;
1096 : lfp++;
1097 3960 : if (*lfp >= 0.0)
1098 4 : lf = *lfp * precision + 0.5;
1099 : else
1100 3956 : lf = *lfp * precision - 0.5;
1101 3960 : if (fabs(lf) > INT_MAX-2)
1102 : {
1103 : /* scaling would cause overflow */
1104 0 : fprintf(stderr,"Internal overflow compressing coordinates.\n");
1105 : errval=0;
1106 : }
1107 3960 : lint2 = lf;
1108 3960 : if (lint2 < minint[1]) minint[1] = lint2;
1109 3960 : if (lint2 > maxint[1]) maxint[1] = lint2;
1110 3960 : *lip++ = lint2;
1111 : lfp++;
1112 3960 : if (*lfp >= 0.0)
1113 3960 : lf = *lfp * precision + 0.5;
1114 : else
1115 0 : lf = *lfp * precision - 0.5;
1116 : if (fabs(lf) > INT_MAX-2)
1117 : {
1118 : errval=0;
1119 : }
1120 3960 : lint3 = lf;
1121 3960 : if (lint3 < minint[2]) minint[2] = lint3;
1122 3960 : if (lint3 > maxint[2]) maxint[2] = lint3;
1123 3960 : *lip++ = lint3;
1124 3960 : lfp++;
1125 3960 : diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
1126 3960 : if (diff < mindiff && lfp > ptr + 3)
1127 : mindiff = diff;
1128 : oldlint1 = lint1;
1129 : oldlint2 = lint2;
1130 : oldlint3 = lint3;
1131 : }
1132 60 : xdrfile_write_int(minint,3,xfp);
1133 60 : xdrfile_write_int(maxint,3,xfp);
1134 :
1135 60 : if ((float)maxint[0] - (float)minint[0] >= INT_MAX-2 ||
1136 60 : (float)maxint[1] - (float)minint[1] >= INT_MAX-2 ||
1137 60 : (float)maxint[2] - (float)minint[2] >= INT_MAX-2) {
1138 : /* turning value in unsigned by subtracting minint
1139 : * would cause overflow
1140 : */
1141 0 : fprintf(stderr,"Internal overflow compressing coordinates.\n");
1142 : errval=0;
1143 : }
1144 60 : sizeint[0] = maxint[0] - minint[0]+1;
1145 60 : sizeint[1] = maxint[1] - minint[1]+1;
1146 60 : sizeint[2] = maxint[2] - minint[2]+1;
1147 :
1148 : /* check if one of the sizes is to big to be multiplied */
1149 60 : if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
1150 : {
1151 0 : bitsizeint[0] = sizeofint(sizeint[0]);
1152 0 : bitsizeint[1] = sizeofint(sizeint[1]);
1153 0 : bitsizeint[2] = sizeofint(sizeint[2]);
1154 : bitsize = 0; /* flag the use of large sizes */
1155 : }
1156 : else
1157 : {
1158 60 : bitsize = sizeofints(3, sizeint);
1159 : }
1160 : lip = buf1;
1161 : luip = (unsigned int *) buf1;
1162 60 : smallidx = FIRSTIDX;
1163 1046 : while (smallidx < LASTIDX && magicints[smallidx] < mindiff)
1164 : {
1165 986 : smallidx++;
1166 : }
1167 60 : xdrfile_write_int(&smallidx,1,xfp);
1168 60 : tmp=smallidx+8;
1169 60 : maxidx = (LASTIDX<tmp) ? LASTIDX : tmp;
1170 60 : minidx = maxidx - 8; /* often this equal smallidx */
1171 60 : tmp=smallidx-1;
1172 60 : tmp= (FIRSTIDX>tmp) ? FIRSTIDX : tmp;
1173 60 : smaller = magicints[tmp] / 2;
1174 60 : smallnum = magicints[smallidx] / 2;
1175 60 : sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1176 60 : larger = magicints[maxidx] / 2;
1177 : i = 0;
1178 814 : while (i < size)
1179 : {
1180 : is_small = 0;
1181 754 : thiscoord = (int *)(luip) + i * 3;
1182 754 : if (smallidx < maxidx && i >= 1 &&
1183 531 : abs(thiscoord[0] - prevcoord[0]) < larger &&
1184 531 : abs(thiscoord[1] - prevcoord[1]) < larger &&
1185 531 : abs(thiscoord[2] - prevcoord[2]) < larger) {
1186 : is_smaller = 1;
1187 : }
1188 223 : else if (smallidx > minidx)
1189 : {
1190 : is_smaller = -1;
1191 : }
1192 : else
1193 : {
1194 : is_smaller = 0;
1195 : }
1196 754 : if (i + 1 < size)
1197 : {
1198 738 : if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
1199 651 : abs(thiscoord[1] - thiscoord[4]) < smallnum &&
1200 580 : abs(thiscoord[2] - thiscoord[5]) < smallnum)
1201 : {
1202 : /* interchange first with second atom for better
1203 : * compression of water molecules
1204 : */
1205 532 : tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
1206 532 : thiscoord[3] = tmp;
1207 532 : tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
1208 532 : thiscoord[4] = tmp;
1209 532 : tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
1210 532 : thiscoord[5] = tmp;
1211 : is_small = 1;
1212 : }
1213 : }
1214 754 : tmpcoord[0] = thiscoord[0] - minint[0];
1215 754 : tmpcoord[1] = thiscoord[1] - minint[1];
1216 754 : tmpcoord[2] = thiscoord[2] - minint[2];
1217 754 : if (bitsize == 0)
1218 : {
1219 0 : encodebits(buf2, bitsizeint[0], tmpcoord[0]);
1220 0 : encodebits(buf2, bitsizeint[1], tmpcoord[1]);
1221 0 : encodebits(buf2, bitsizeint[2], tmpcoord[2]);
1222 : }
1223 : else
1224 : {
1225 754 : encodeints(buf2, 3, bitsize, sizeint, tmpcoord);
1226 : }
1227 754 : prevcoord[0] = thiscoord[0];
1228 754 : prevcoord[1] = thiscoord[1];
1229 754 : prevcoord[2] = thiscoord[2];
1230 754 : thiscoord = thiscoord + 3;
1231 : i++;
1232 :
1233 : run = 0;
1234 754 : if (is_small == 0 && is_smaller == -1)
1235 : is_smaller = 0;
1236 3960 : while (is_small && run < 8*3)
1237 : {
1238 : tmpsum=0;
1239 12824 : for(j=0;j<3;j++)
1240 : {
1241 9618 : tmp=thiscoord[j] - prevcoord[j];
1242 9618 : tmpsum+=tmp*tmp;
1243 : }
1244 3206 : if (is_smaller == -1 && tmpsum >= smaller * smaller)
1245 : {
1246 : is_smaller = 0;
1247 : }
1248 :
1249 3206 : tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
1250 3206 : tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
1251 3206 : tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
1252 :
1253 3206 : prevcoord[0] = thiscoord[0];
1254 3206 : prevcoord[1] = thiscoord[1];
1255 3206 : prevcoord[2] = thiscoord[2];
1256 :
1257 3206 : i++;
1258 3206 : thiscoord = thiscoord + 3;
1259 : is_small = 0;
1260 3206 : if (i < size &&
1261 3162 : abs(thiscoord[0] - prevcoord[0]) < smallnum &&
1262 3083 : abs(thiscoord[1] - prevcoord[1]) < smallnum &&
1263 3027 : abs(thiscoord[2] - prevcoord[2]) < smallnum)
1264 : {
1265 : is_small = 1;
1266 : }
1267 : }
1268 754 : if (run != prevrun || is_smaller != 0)
1269 : {
1270 : prevrun = run;
1271 717 : encodebits(buf2, 1, 1); /* flag the change in run-length */
1272 717 : encodebits(buf2, 5, run+is_smaller+1);
1273 : }
1274 : else
1275 : {
1276 37 : encodebits(buf2, 1, 0); /* flag the fact that runlength did not change */
1277 : }
1278 3960 : for (k=0; k < run; k+=3)
1279 : {
1280 3206 : encodeints(buf2, 3, smallidx, sizesmall, &tmpcoord[k]);
1281 : }
1282 754 : if (is_smaller != 0)
1283 : {
1284 652 : smallidx += is_smaller;
1285 652 : if (is_smaller < 0)
1286 : {
1287 : smallnum = smaller;
1288 121 : smaller = magicints[smallidx-1] / 2;
1289 : }
1290 : else
1291 : {
1292 : smaller = smallnum;
1293 531 : smallnum = magicints[smallidx] / 2;
1294 : }
1295 652 : sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1296 : }
1297 : }
1298 60 : if (buf2[1] != 0) buf2[0]++;
1299 60 : xdrfile_write_int(buf2,1,xfp); /* buf2[0] holds the length in bytes */
1300 60 : tmp=xdrfile_write_opaque((char *)&(buf2[3]),(unsigned int)buf2[0],xfp);
1301 60 : if(tmp==(unsigned int)buf2[0])
1302 60 : return size;
1303 : else
1304 : return -1;
1305 : }
1306 :
1307 :
1308 : int
1309 0 : xdrfile_decompress_coord_double(double *ptr,
1310 : int *size,
1311 : double *precision,
1312 : XDRFILE* xfp)
1313 : {
1314 : int minint[3], maxint[3], *lip;
1315 : int smallidx, minidx, maxidx;
1316 : unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3;
1317 : int k, *buf1, *buf2, lsize, flag;
1318 : int smallnum, smaller, larger, i, is_smaller, run;
1319 : double *lfp, inv_precision;
1320 : float float_prec, tmpdata[30];
1321 : int tmp, *thiscoord, prevcoord[3];
1322 : unsigned int bitsize;
1323 :
1324 : bitsizeint[0] = 0;
1325 : bitsizeint[1] = 0;
1326 : bitsizeint[2] = 0;
1327 :
1328 0 : if(xfp==NULL || ptr==NULL)
1329 : return -1;
1330 0 : tmp=xdrfile_read_int(&lsize,1,xfp);
1331 0 : if(tmp==0)
1332 : return -1; /* return if we could not read size */
1333 0 : if (*size < lsize)
1334 : {
1335 0 : fprintf(stderr, "Requested to decompress %d coords, file contains %d\n",
1336 : *size, lsize);
1337 0 : return -1;
1338 : }
1339 0 : *size = lsize;
1340 0 : size3 = *size * 3;
1341 0 : if(size3>xfp->buf1size)
1342 : {
1343 0 : if((xfp->buf1=(int *)malloc(sizeof(int)*size3))==NULL)
1344 : {
1345 0 : fprintf(stderr,"Cannot allocate memory for decompression coordinates.\n");
1346 0 : return -1;
1347 : }
1348 0 : xfp->buf1size=size3;
1349 0 : xfp->buf2size=size3*1.2;
1350 0 : if((xfp->buf2=(int *)malloc(sizeof(int)*xfp->buf2size))==NULL)
1351 : {
1352 0 : fprintf(stderr,"Cannot allocate memory for decompressing coordinates.\n");
1353 0 : return -1;
1354 : }
1355 : }
1356 : /* Dont bother with compression for three atoms or less */
1357 0 : if(*size<=9)
1358 : {
1359 0 : tmp=xdrfile_read_float(tmpdata,size3,xfp);
1360 0 : for(i=0;i<9*3;i++)
1361 0 : ptr[i]=tmpdata[i];
1362 0 : return tmp/3;
1363 : /* return number of coords, not floats */
1364 : }
1365 : /* Compression-time if we got here. Read precision first */
1366 0 : xdrfile_read_float(&float_prec,1,xfp);
1367 0 : *precision=float_prec;
1368 : /* avoid repeated pointer dereferencing. */
1369 0 : buf1=xfp->buf1;
1370 0 : buf2=xfp->buf2;
1371 : /* buf2[0-2] are special and do not contain actual data */
1372 0 : buf2[0] = buf2[1] = buf2[2] = 0;
1373 0 : xdrfile_read_int(minint,3,xfp);
1374 0 : xdrfile_read_int(maxint,3,xfp);
1375 :
1376 0 : sizeint[0] = maxint[0] - minint[0]+1;
1377 0 : sizeint[1] = maxint[1] - minint[1]+1;
1378 0 : sizeint[2] = maxint[2] - minint[2]+1;
1379 :
1380 : /* check if one of the sizes is to big to be multiplied */
1381 0 : if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
1382 : {
1383 0 : bitsizeint[0] = sizeofint(sizeint[0]);
1384 0 : bitsizeint[1] = sizeofint(sizeint[1]);
1385 0 : bitsizeint[2] = sizeofint(sizeint[2]);
1386 : bitsize = 0; /* flag the use of large sizes */
1387 : }
1388 : else
1389 : {
1390 0 : bitsize = sizeofints(3, sizeint);
1391 : }
1392 :
1393 0 : if (xdrfile_read_int(&smallidx,1,xfp) == 0)
1394 : return 0;
1395 0 : tmp=smallidx+8;
1396 : maxidx = (LASTIDX<tmp) ? LASTIDX : tmp;
1397 : minidx = maxidx - 8; /* often this equal smallidx */
1398 0 : tmp = smallidx-1;
1399 0 : tmp = (FIRSTIDX>tmp) ? FIRSTIDX : tmp;
1400 0 : smaller = magicints[tmp] / 2;
1401 0 : smallnum = magicints[smallidx] / 2;
1402 0 : sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1403 : larger = magicints[maxidx];
1404 :
1405 : /* buf2[0] holds the length in bytes */
1406 :
1407 0 : if (xdrfile_read_int(buf2,1,xfp) == 0)
1408 : return 0;
1409 0 : if (xdrfile_read_opaque((char *)&(buf2[3]),(unsigned int)buf2[0],xfp) == 0)
1410 : return 0;
1411 0 : buf2[0] = buf2[1] = buf2[2] = 0;
1412 :
1413 : lfp = ptr;
1414 0 : inv_precision = 1.0 / * precision;
1415 : run = 0;
1416 : i = 0;
1417 : lip = buf1;
1418 0 : while ( i < lsize )
1419 : {
1420 0 : thiscoord = (int *)(lip) + i * 3;
1421 :
1422 0 : if (bitsize == 0)
1423 : {
1424 0 : thiscoord[0] = decodebits(buf2, bitsizeint[0]);
1425 0 : thiscoord[1] = decodebits(buf2, bitsizeint[1]);
1426 0 : thiscoord[2] = decodebits(buf2, bitsizeint[2]);
1427 : } else {
1428 0 : decodeints(buf2, 3, bitsize, sizeint, thiscoord);
1429 : }
1430 :
1431 0 : i++;
1432 0 : thiscoord[0] += minint[0];
1433 0 : thiscoord[1] += minint[1];
1434 0 : thiscoord[2] += minint[2];
1435 :
1436 : prevcoord[0] = thiscoord[0];
1437 : prevcoord[1] = thiscoord[1];
1438 : prevcoord[2] = thiscoord[2];
1439 :
1440 0 : flag = decodebits(buf2, 1);
1441 : is_smaller = 0;
1442 0 : if (flag == 1)
1443 : {
1444 0 : run = decodebits(buf2, 5);
1445 0 : is_smaller = run % 3;
1446 0 : run -= is_smaller;
1447 0 : is_smaller--;
1448 : }
1449 0 : if (run > 0)
1450 : {
1451 0 : thiscoord += 3;
1452 0 : for (k = 0; k < run; k+=3)
1453 : {
1454 0 : decodeints(buf2, 3, smallidx, sizesmall, thiscoord);
1455 0 : i++;
1456 0 : thiscoord[0] += prevcoord[0] - smallnum;
1457 0 : thiscoord[1] += prevcoord[1] - smallnum;
1458 0 : thiscoord[2] += prevcoord[2] - smallnum;
1459 0 : if (k == 0)
1460 : {
1461 : /* interchange first with second atom for better
1462 : * compression of water molecules
1463 : */
1464 0 : tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1465 : prevcoord[0] = tmp;
1466 0 : tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1467 : prevcoord[1] = tmp;
1468 0 : tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1469 : prevcoord[2] = tmp;
1470 0 : *lfp++ = prevcoord[0] * inv_precision;
1471 0 : *lfp++ = prevcoord[1] * inv_precision;
1472 0 : *lfp++ = prevcoord[2] * inv_precision;
1473 : }
1474 : else
1475 : {
1476 : prevcoord[0] = thiscoord[0];
1477 : prevcoord[1] = thiscoord[1];
1478 : prevcoord[2] = thiscoord[2];
1479 : }
1480 0 : *lfp++ = thiscoord[0] * inv_precision;
1481 0 : *lfp++ = thiscoord[1] * inv_precision;
1482 0 : *lfp++ = thiscoord[2] * inv_precision;
1483 : }
1484 : } else {
1485 0 : *lfp++ = thiscoord[0] * inv_precision;
1486 0 : *lfp++ = thiscoord[1] * inv_precision;
1487 0 : *lfp++ = thiscoord[2] * inv_precision;
1488 : }
1489 0 : smallidx += is_smaller;
1490 0 : if (is_smaller < 0) {
1491 : smallnum = smaller;
1492 0 : if (smallidx > FIRSTIDX) {
1493 0 : smaller = magicints[smallidx - 1] /2;
1494 : } else {
1495 : smaller = 0;
1496 : }
1497 0 : } else if (is_smaller > 0) {
1498 : smaller = smallnum;
1499 0 : smallnum = magicints[smallidx] / 2;
1500 : }
1501 0 : sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1502 : }
1503 0 : return *size;
1504 : }
1505 :
1506 : int
1507 0 : xdrfile_compress_coord_double(double *ptr,
1508 : int size,
1509 : double precision,
1510 : XDRFILE* xfp)
1511 : {
1512 : int minint[3], maxint[3], mindiff, *lip, diff;
1513 : int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
1514 : int minidx, maxidx;
1515 : unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
1516 : int k, *buf1, *buf2;
1517 : int smallnum, smaller, larger, i, j, is_small, is_smaller, run, prevrun;
1518 : double *lfp;
1519 : float float_prec, lf,tmpdata[30];
1520 : int tmp, tmpsum, *thiscoord, prevcoord[3];
1521 : unsigned int tmpcoord[30];
1522 : int errval=1;
1523 : unsigned int bitsize;
1524 :
1525 : bitsizeint[0] = 0;
1526 : bitsizeint[1] = 0;
1527 : bitsizeint[2] = 0;
1528 :
1529 0 : if(xfp==NULL)
1530 : return -1;
1531 0 : size3=3*size;
1532 0 : if(size3>xfp->buf1size) {
1533 0 : if((xfp->buf1=(int *)malloc(sizeof(int)*size3))==NULL) {
1534 0 : fprintf(stderr,"Cannot allocate memory for compressing coordinates.\n");
1535 0 : return -1;
1536 : }
1537 0 : xfp->buf1size=size3;
1538 0 : xfp->buf2size=size3*1.2;
1539 0 : if((xfp->buf2=(int *)malloc(sizeof(int)*xfp->buf2size))==NULL) {
1540 0 : fprintf(stderr,"Cannot allocate memory for compressing coordinates.\n");
1541 0 : return -1;
1542 : }
1543 : }
1544 0 : if(xdrfile_write_int(&size,1,xfp)==0)
1545 : return -1; /* return if we could not write size */
1546 : /* Dont bother with compression for three atoms or less */
1547 0 : if(size<=9) {
1548 0 : for(i=0;i<9*3;i++)
1549 0 : tmpdata[i]=ptr[i];
1550 0 : return xdrfile_write_float(tmpdata,size3,xfp)/3;
1551 : /* return number of coords, not floats */
1552 : }
1553 : /* Compression-time if we got here. Write precision first */
1554 0 : if (precision <= 0)
1555 : precision = 1000;
1556 0 : float_prec=precision;
1557 0 : xdrfile_write_float(&float_prec,1,xfp);
1558 : /* avoid repeated pointer dereferencing. */
1559 0 : buf1=xfp->buf1;
1560 0 : buf2=xfp->buf2;
1561 : /* buf2[0-2] are special and do not contain actual data */
1562 0 : buf2[0] = buf2[1] = buf2[2] = 0;
1563 0 : minint[0] = minint[1] = minint[2] = INT_MAX;
1564 0 : maxint[0] = maxint[1] = maxint[2] = INT_MIN;
1565 : prevrun = -1;
1566 : lfp = ptr;
1567 : lip = buf1;
1568 : mindiff = INT_MAX;
1569 : oldlint1 = oldlint2 = oldlint3 = 0;
1570 0 : while(lfp < ptr + size3 ) {
1571 : /* find nearest integer */
1572 0 : if (*lfp >= 0.0)
1573 0 : lf = (float)*lfp * float_prec + 0.5;
1574 : else
1575 0 : lf = (float)*lfp * float_prec - 0.5;
1576 0 : if (fabs(lf) > INT_MAX-2) {
1577 : /* scaling would cause overflow */
1578 0 : fprintf(stderr,"Internal overflow compressing coordinates.\n");
1579 : errval=0;
1580 : }
1581 0 : lint1 = lf;
1582 0 : if (lint1 < minint[0]) minint[0] = lint1;
1583 0 : if (lint1 > maxint[0]) maxint[0] = lint1;
1584 0 : *lip++ = lint1;
1585 : lfp++;
1586 0 : if (*lfp >= 0.0)
1587 0 : lf = (float)*lfp * float_prec + 0.5;
1588 : else
1589 0 : lf = (float)*lfp * float_prec - 0.5;
1590 0 : if (fabs(lf) > INT_MAX-2) {
1591 : /* scaling would cause overflow */
1592 0 : fprintf(stderr,"Internal overflow compressing coordinates.\n");
1593 : errval=0;
1594 : }
1595 0 : lint2 = lf;
1596 0 : if (lint2 < minint[1]) minint[1] = lint2;
1597 0 : if (lint2 > maxint[1]) maxint[1] = lint2;
1598 0 : *lip++ = lint2;
1599 : lfp++;
1600 0 : if (*lfp >= 0.0)
1601 0 : lf = (float)*lfp * float_prec + 0.5;
1602 : else
1603 0 : lf = (float)*lfp * float_prec - 0.5;
1604 : if (fabs(lf) > INT_MAX-2) {
1605 : errval=0;
1606 : }
1607 0 : lint3 = lf;
1608 0 : if (lint3 < minint[2]) minint[2] = lint3;
1609 0 : if (lint3 > maxint[2]) maxint[2] = lint3;
1610 0 : *lip++ = lint3;
1611 0 : lfp++;
1612 0 : diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
1613 0 : if (diff < mindiff && lfp > ptr + 3)
1614 : mindiff = diff;
1615 : oldlint1 = lint1;
1616 : oldlint2 = lint2;
1617 : oldlint3 = lint3;
1618 : }
1619 0 : xdrfile_write_int(minint,3,xfp);
1620 0 : xdrfile_write_int(maxint,3,xfp);
1621 :
1622 0 : if ((float)maxint[0] - (float)minint[0] >= INT_MAX-2 ||
1623 0 : (float)maxint[1] - (float)minint[1] >= INT_MAX-2 ||
1624 0 : (float)maxint[2] - (float)minint[2] >= INT_MAX-2) {
1625 : /* turning value in unsigned by subtracting minint
1626 : * would cause overflow
1627 : */
1628 0 : fprintf(stderr,"Internal overflow compressing coordinates.\n");
1629 : errval=0;
1630 : }
1631 0 : sizeint[0] = maxint[0] - minint[0]+1;
1632 0 : sizeint[1] = maxint[1] - minint[1]+1;
1633 0 : sizeint[2] = maxint[2] - minint[2]+1;
1634 :
1635 : /* check if one of the sizes is to big to be multiplied */
1636 0 : if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
1637 0 : bitsizeint[0] = sizeofint(sizeint[0]);
1638 0 : bitsizeint[1] = sizeofint(sizeint[1]);
1639 0 : bitsizeint[2] = sizeofint(sizeint[2]);
1640 : bitsize = 0; /* flag the use of large sizes */
1641 : } else {
1642 0 : bitsize = sizeofints(3, sizeint);
1643 : }
1644 : lip = buf1;
1645 : luip = (unsigned int *) buf1;
1646 0 : smallidx = FIRSTIDX;
1647 0 : while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
1648 0 : smallidx++;
1649 : }
1650 0 : xdrfile_write_int(&smallidx,1,xfp);
1651 0 : tmp=smallidx+8;
1652 0 : maxidx = (LASTIDX<tmp) ? LASTIDX : tmp;
1653 0 : minidx = maxidx - 8; /* often this equal smallidx */
1654 0 : tmp=smallidx-1;
1655 0 : tmp= (FIRSTIDX>tmp) ? FIRSTIDX : tmp;
1656 0 : smaller = magicints[tmp] / 2;
1657 0 : smallnum = magicints[smallidx] / 2;
1658 0 : sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1659 0 : larger = magicints[maxidx] / 2;
1660 : i = 0;
1661 0 : while (i < size) {
1662 : is_small = 0;
1663 0 : thiscoord = (int *)(luip) + i * 3;
1664 0 : if (smallidx < maxidx && i >= 1 &&
1665 0 : abs(thiscoord[0] - prevcoord[0]) < larger &&
1666 0 : abs(thiscoord[1] - prevcoord[1]) < larger &&
1667 0 : abs(thiscoord[2] - prevcoord[2]) < larger) {
1668 : is_smaller = 1;
1669 0 : } else if (smallidx > minidx) {
1670 : is_smaller = -1;
1671 : } else {
1672 : is_smaller = 0;
1673 : }
1674 0 : if (i + 1 < size) {
1675 0 : if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
1676 0 : abs(thiscoord[1] - thiscoord[4]) < smallnum &&
1677 0 : abs(thiscoord[2] - thiscoord[5]) < smallnum) {
1678 : /* interchange first with second atom for better
1679 : * compression of water molecules
1680 : */
1681 0 : tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
1682 0 : thiscoord[3] = tmp;
1683 0 : tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
1684 0 : thiscoord[4] = tmp;
1685 0 : tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
1686 0 : thiscoord[5] = tmp;
1687 : is_small = 1;
1688 : }
1689 : }
1690 0 : tmpcoord[0] = thiscoord[0] - minint[0];
1691 0 : tmpcoord[1] = thiscoord[1] - minint[1];
1692 0 : tmpcoord[2] = thiscoord[2] - minint[2];
1693 0 : if (bitsize == 0) {
1694 0 : encodebits(buf2, bitsizeint[0], tmpcoord[0]);
1695 0 : encodebits(buf2, bitsizeint[1], tmpcoord[1]);
1696 0 : encodebits(buf2, bitsizeint[2], tmpcoord[2]);
1697 : } else {
1698 0 : encodeints(buf2, 3, bitsize, sizeint, tmpcoord);
1699 : }
1700 0 : prevcoord[0] = thiscoord[0];
1701 0 : prevcoord[1] = thiscoord[1];
1702 0 : prevcoord[2] = thiscoord[2];
1703 0 : thiscoord = thiscoord + 3;
1704 : i++;
1705 :
1706 : run = 0;
1707 0 : if (is_small == 0 && is_smaller == -1)
1708 : is_smaller = 0;
1709 0 : while (is_small && run < 8*3) {
1710 : tmpsum=0;
1711 0 : for(j=0;j<3;j++) {
1712 0 : tmp=thiscoord[j] - prevcoord[j];
1713 0 : tmpsum+=tmp*tmp;
1714 : }
1715 0 : if (is_smaller == -1 && tmpsum >= smaller * smaller) {
1716 : is_smaller = 0;
1717 : }
1718 :
1719 0 : tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
1720 0 : tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
1721 0 : tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
1722 :
1723 0 : prevcoord[0] = thiscoord[0];
1724 0 : prevcoord[1] = thiscoord[1];
1725 0 : prevcoord[2] = thiscoord[2];
1726 :
1727 0 : i++;
1728 0 : thiscoord = thiscoord + 3;
1729 : is_small = 0;
1730 0 : if (i < size &&
1731 0 : abs(thiscoord[0] - prevcoord[0]) < smallnum &&
1732 0 : abs(thiscoord[1] - prevcoord[1]) < smallnum &&
1733 0 : abs(thiscoord[2] - prevcoord[2]) < smallnum) {
1734 : is_small = 1;
1735 : }
1736 : }
1737 0 : if (run != prevrun || is_smaller != 0) {
1738 : prevrun = run;
1739 0 : encodebits(buf2, 1, 1); /* flag the change in run-length */
1740 0 : encodebits(buf2, 5, run+is_smaller+1);
1741 : } else {
1742 0 : encodebits(buf2, 1, 0); /* flag the fact that runlength did not change */
1743 : }
1744 0 : for (k=0; k < run; k+=3) {
1745 0 : encodeints(buf2, 3, smallidx, sizesmall, &tmpcoord[k]);
1746 : }
1747 0 : if (is_smaller != 0) {
1748 0 : smallidx += is_smaller;
1749 0 : if (is_smaller < 0) {
1750 : smallnum = smaller;
1751 0 : smaller = magicints[smallidx-1] / 2;
1752 : } else {
1753 : smaller = smallnum;
1754 0 : smallnum = magicints[smallidx] / 2;
1755 : }
1756 0 : sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1757 : }
1758 : }
1759 0 : if (buf2[1] != 0) buf2[0]++;
1760 0 : xdrfile_write_int(buf2,1,xfp); /* buf2[0] holds the length in bytes */
1761 0 : tmp=xdrfile_write_opaque((char *)&(buf2[3]),(unsigned int)buf2[0],xfp);
1762 0 : if(tmp==(unsigned int)buf2[0])
1763 0 : return size;
1764 : else
1765 : return -1;
1766 : }
1767 :
1768 :
1769 : /* Dont try do document Fortran interface, since
1770 : * Doxygen barfs at the F77_FUNC macro
1771 : */
1772 : #ifndef DOXYGEN
1773 :
1774 : /*************************************************************
1775 : * Fortran77 interface for reading/writing portable data *
1776 : * The routine are not threadsafe when called from Fortran *
1777 : * (as they are when called from C) unless you compile with *
1778 : * this file with posix thread support. *
1779 : * Note that these are not multithread-safe. *
1780 : *************************************************************/
1781 : #define MAX_FORTRAN_XDR 1024
1782 : static XDRFILE *f77xdr[MAX_FORTRAN_XDR]; /* array of file handles */
1783 : static int f77init = 1; /* zero array first time */
1784 :
1785 : /* internal to this file: C<-->Fortran string conversion */
1786 : static int ftocstr(char *dest, int dest_len, char *src, int src_len);
1787 : static int ctofstr(char *dest, int dest_len, char *src);
1788 :
1789 :
1790 : void
1791 0 : F77_FUNC(xdropen,XDROPEN)(int *fid, char *filename, char *mode,
1792 : int fn_len, int mode_len)
1793 : {
1794 : char cfilename[512];
1795 : char cmode[5];
1796 : int i;
1797 :
1798 : /* zero array at first invocation */
1799 0 : if(f77init) {
1800 0 : for(i=0;i<MAX_FORTRAN_XDR;i++)
1801 0 : f77xdr[i]=NULL;
1802 0 : f77init=0;
1803 : }
1804 : i=0;
1805 :
1806 : /* nf77xdr is always smaller or equal to MAX_FORTRAN_XDR */
1807 0 : while(i<MAX_FORTRAN_XDR && f77xdr[i]!=NULL)
1808 0 : i++;
1809 0 : if(i==MAX_FORTRAN_XDR) {
1810 0 : *fid = -1;
1811 0 : } else if (ftocstr(cfilename, sizeof(cfilename), filename, fn_len)) {
1812 0 : *fid = -1;
1813 0 : } else if (ftocstr(cmode, sizeof(cmode), mode,mode_len)) {
1814 0 : *fid = -1;
1815 : } else {
1816 0 : f77xdr[i]=xdrfile_open(cfilename,cmode);
1817 : /* return the index in the array as a fortran file handle */
1818 0 : *fid=i;
1819 : }
1820 0 : }
1821 :
1822 : void
1823 0 : F77_FUNC(xdrclose,XDRCLOSE)(int *fid)
1824 : {
1825 : /* first close it */
1826 0 : xdrfile_close(f77xdr[*fid]);
1827 : /* the remove it from file handle list */
1828 0 : f77xdr[*fid]=NULL;
1829 0 : }
1830 :
1831 :
1832 : void
1833 0 : F77_FUNC(xdrrint,XDRRINT)(int *fid, int *data, int *ndata, int *ret)
1834 : {
1835 0 : *ret = xdrfile_read_int(data,*ndata,f77xdr[*fid]);
1836 0 : }
1837 :
1838 : void
1839 0 : F77_FUNC(xdrwint,XDRWINT)(int *fid, int *data, int *ndata, int *ret)
1840 : {
1841 0 : *ret = xdrfile_write_int(data,*ndata,f77xdr[*fid]);
1842 0 : }
1843 :
1844 : void
1845 0 : F77_FUNC(xdrruint,XDRRUINT)(int *fid, unsigned int *data, int *ndata, int *ret)
1846 : {
1847 0 : *ret = xdrfile_read_uint(data,*ndata,f77xdr[*fid]);
1848 0 : }
1849 :
1850 : void
1851 0 : F77_FUNC(xdrwuint,XDRWUINT)(int *fid, unsigned int *data, int *ndata, int *ret)
1852 : {
1853 0 : *ret = xdrfile_write_uint(data,*ndata,f77xdr[*fid]);
1854 0 : }
1855 :
1856 : void
1857 0 : F77_FUNC(xdrrchar,XDRRCHAR)(int *fid, char *ip, int *ndata, int *ret)
1858 : {
1859 0 : *ret = xdrfile_read_char(ip,*ndata,f77xdr[*fid]);
1860 0 : }
1861 :
1862 : void
1863 0 : F77_FUNC(xdrwchar,XDRWCHAR)(int *fid, char *ip, int *ndata, int *ret)
1864 : {
1865 0 : *ret = xdrfile_write_char(ip,*ndata,f77xdr[*fid]);
1866 0 : }
1867 :
1868 : void
1869 0 : F77_FUNC(xdrruchar,XDRRUCHAR)(int *fid, unsigned char *ip, int *ndata, int *ret)
1870 : {
1871 0 : *ret = xdrfile_read_uchar(ip,*ndata,f77xdr[*fid]);
1872 0 : }
1873 :
1874 : void
1875 0 : F77_FUNC(xdrwuchar,XDRWUCHAR)(int *fid, unsigned char *ip, int *ndata, int *ret)
1876 : {
1877 0 : *ret = xdrfile_write_uchar(ip,*ndata,f77xdr[*fid]);
1878 0 : }
1879 :
1880 : void
1881 0 : F77_FUNC(xdrrshort,XDRRSHORT)(int *fid, short *ip, int *ndata, int *ret)
1882 : {
1883 0 : *ret = xdrfile_read_short(ip,*ndata,f77xdr[*fid]);
1884 0 : }
1885 :
1886 : void
1887 0 : F77_FUNC(xdrwshort,XDRWSHORT)(int *fid, short *ip, int *ndata, int *ret)
1888 : {
1889 0 : *ret = xdrfile_write_short(ip,*ndata,f77xdr[*fid]);
1890 0 : }
1891 :
1892 : void
1893 0 : F77_FUNC(xdrrushort,XDRRUSHORT)(int *fid, unsigned short *ip, int *ndata, int *ret)
1894 : {
1895 0 : *ret = xdrfile_read_ushort(ip,*ndata,f77xdr[*fid]);
1896 0 : }
1897 :
1898 : void
1899 0 : F77_FUNC(xdrwushort,XDRWUSHORT)(int *fid, unsigned short *ip, int *ndata, int *ret)
1900 : {
1901 0 : *ret = xdrfile_write_ushort(ip,*ndata,f77xdr[*fid]);
1902 0 : }
1903 :
1904 : void
1905 0 : F77_FUNC(xdrrsingle,XDRRSINGLE)(int *fid, float *data, int *ndata, int *ret)
1906 : {
1907 0 : *ret = xdrfile_read_float(data,*ndata,f77xdr[*fid]);
1908 0 : }
1909 :
1910 : void
1911 0 : F77_FUNC(xdrwsingle,XDRWSINGLE)(int *fid, float *data, int *ndata, int *ret)
1912 : {
1913 0 : *ret = xdrfile_write_float(data,*ndata,f77xdr[*fid]);
1914 0 : }
1915 :
1916 : void
1917 0 : F77_FUNC(xdrrdouble,XDRRDOUBLE)(int *fid, double *data, int *ndata, int *ret)
1918 : {
1919 0 : *ret = xdrfile_read_double(data,*ndata,f77xdr[*fid]);
1920 0 : }
1921 :
1922 : void
1923 0 : F77_FUNC(xdrwdouble,XDRWDOUBLE)(int *fid, double *data, int *ndata, int *ret)
1924 : {
1925 0 : *ret = xdrfile_write_double(data,*ndata,f77xdr[*fid]);
1926 0 : }
1927 :
1928 0 : static int ftocstr(char *dest, int destlen, char *src, int srclen)
1929 : {
1930 : char *p;
1931 :
1932 0 : p = src + srclen;
1933 0 : while ( --p >= src && *p == ' ' );
1934 0 : srclen = p - src + 1;
1935 : destlen--;
1936 0 : dest[0] = 0;
1937 0 : if (srclen > destlen)
1938 : return 1;
1939 0 : while (srclen--)
1940 0 : (*dest++ = *src++);
1941 0 : *dest = '\0';
1942 0 : return 0;
1943 : }
1944 :
1945 :
1946 0 : static int ctofstr(char *dest, int destlen, char *src)
1947 : {
1948 0 : while (destlen && *src) {
1949 0 : *dest++ = *src++;
1950 0 : destlen--;
1951 : }
1952 0 : while (destlen--)
1953 0 : *dest++ = ' ';
1954 0 : return 0;
1955 : }
1956 :
1957 :
1958 : void
1959 0 : F77_FUNC(xdrrstring,XDRRSTRING)(int *fid, char *str, int *ret, int len)
1960 : {
1961 : char *cstr;
1962 :
1963 0 : if((cstr=(char*)malloc((len+1)*sizeof(char)))==NULL) {
1964 0 : *ret = 0;
1965 0 : return;
1966 : }
1967 0 : if (ftocstr(cstr, len+1, str, len)) {
1968 0 : *ret = 0;
1969 0 : free(cstr);
1970 0 : return;
1971 : }
1972 :
1973 0 : *ret = xdrfile_read_string(cstr, len+1,f77xdr[*fid]);
1974 0 : ctofstr( str, len , cstr);
1975 0 : free(cstr);
1976 : }
1977 :
1978 : void
1979 0 : F77_FUNC(xdrwstring,XDRWSTRING)(int *fid, char *str, int *ret, int len)
1980 : {
1981 : char *cstr;
1982 :
1983 0 : if((cstr=(char*)malloc((len+1)*sizeof(char)))==NULL) {
1984 0 : *ret = 0;
1985 0 : return;
1986 : }
1987 0 : if (ftocstr(cstr, len+1, str, len)) {
1988 0 : *ret = 0;
1989 0 : free(cstr);
1990 0 : return;
1991 : }
1992 :
1993 0 : *ret = xdrfile_write_string(cstr, f77xdr[*fid]);
1994 0 : ctofstr( str, len , cstr);
1995 0 : free(cstr);
1996 : }
1997 :
1998 : void
1999 0 : F77_FUNC(xdrropaque,XDRROPAQUE)(int *fid, char *data, int *ndata, int *ret)
2000 : {
2001 0 : *ret = xdrfile_read_opaque(data,*ndata,f77xdr[*fid]);
2002 0 : }
2003 :
2004 : void
2005 0 : F77_FUNC(xdrwopaque,XDRWOPAQUE)(int *fid, char *data, int *ndata, int *ret)
2006 : {
2007 0 : *ret = xdrfile_write_opaque(data,*ndata,f77xdr[*fid]);
2008 0 : }
2009 :
2010 :
2011 : /* Write single-precision compressed 3d coordinates */
2012 : void
2013 0 : F77_FUNC(xdrccs,XDRCCS)(int *fid, float *data, int *ncoord,
2014 : float *precision, int *ret)
2015 : {
2016 0 : *ret = xdrfile_compress_coord_float(data,*ncoord,*precision,f77xdr[*fid]);
2017 0 : }
2018 :
2019 :
2020 : /* Read single-precision compressed 3d coordinates */
2021 : void
2022 0 : F77_FUNC(xdrdcs,XDRDCS)(int *fid, float *data, int *ncoord,
2023 : float *precision, int *ret)
2024 : {
2025 0 : *ret = xdrfile_decompress_coord_float(data,ncoord,precision,f77xdr[*fid]);
2026 0 : }
2027 :
2028 :
2029 : /* Write compressed 3d coordinates from double precision data */
2030 : void
2031 0 : F77_FUNC(xdrccd,XDRCCD)(int *fid, double *data, int *ncoord,
2032 : double *precision, int *ret)
2033 : {
2034 0 : *ret = xdrfile_compress_coord_double(data,*ncoord,*precision,f77xdr[*fid]);
2035 0 : }
2036 :
2037 : /* Read compressed 3d coordinates into double precision data */
2038 : void
2039 0 : F77_FUNC(xddcd,XDRDCD)(int *fid, double *data, int *ncoord,
2040 : double *precision, int *ret)
2041 : {
2042 0 : *ret = xdrfile_decompress_coord_double(data,ncoord,precision,f77xdr[*fid]);
2043 0 : }
2044 :
2045 :
2046 :
2047 :
2048 :
2049 :
2050 :
2051 : #endif /* DOXYGEN */
2052 :
2053 : /*************************************************************
2054 : * End of higher-level routines - dont change things below! *
2055 : *************************************************************/
2056 :
2057 :
2058 :
2059 :
2060 :
2061 :
2062 :
2063 :
2064 :
2065 :
2066 :
2067 :
2068 :
2069 :
2070 :
2071 :
2072 :
2073 :
2074 :
2075 :
2076 : /*************************************************************
2077 : * The rest of this file contains our own implementation *
2078 : * of the XDR calls in case you are compiling without them. *
2079 : * You do NOT want to change things here since it would make *
2080 : * things incompatible with the standard RPC/XDR routines. *
2081 : *************************************************************/
2082 : #ifndef HAVE_RPC_XDR_H
2083 :
2084 : /*
2085 : * What follows is a modified version of the Sun XDR code. For reference
2086 : * we include their copyright and license:
2087 : *
2088 : * Sun RPC is a product of Sun Microsystems, Inc. and is provided for
2089 : * unrestricted use provided that this legend is included on all tape
2090 : * media and as a part of the software program in whole or part. Users
2091 : * may copy or modify Sun RPC without charge, but are not authorized
2092 : * to license or distribute it to anyone else except as part of a product or
2093 : * program developed by the user.
2094 : *
2095 : * SUN RPC IS PROVIDED AS IS WITH NO WARRANTIES OF ANY KIND INCLUDING THE
2096 : * WARRANTIES OF DESIGN, MERCHANTIBILITY AND FITNESS FOR A PARTICULAR
2097 : * PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE PRACTICE.
2098 : *
2099 : * Sun RPC is provided with no support and without any obligation on the
2100 : * part of Sun Microsystems, Inc. to assist in its use, correction,
2101 : * modification or enhancement.
2102 : *
2103 : * SUN MICROSYSTEMS, INC. SHALL HAVE NO LIABILITY WITH RESPECT TO THE
2104 : * INFRINGEMENT OF COPYRIGHTS, TRADE SECRETS OR ANY PATENTS BY SUN RPC
2105 : * OR ANY PART THEREOF.
2106 : *
2107 : * In no event will Sun Microsystems, Inc. be liable for any lost revenue
2108 : * or profits or other special, indirect and consequential damages, even if
2109 : * Sun has been advised of the possibility of such damages.
2110 : *
2111 : * Sun Microsystems, Inc.
2112 : * 2550 Garcia Avenue
2113 : * Mountain View, California 94043
2114 : */
2115 :
2116 : /* INT_MAX is defined in limits.h according to ANSI C */
2117 : #if (INT_MAX > 2147483647)
2118 : # error Error: Cannot use builtin XDR support when size of int
2119 : # error is larger than 4 bytes. Use your system XDR libraries
2120 : # error instead, or modify the source code in xdrfile.c
2121 : #endif /* Check for 4 byte int type */
2122 :
2123 :
2124 :
2125 :
2126 :
2127 : typedef int (*xdrproc_t) (XDR *, void *,...);
2128 :
2129 : #define xdr_getlong(xdrs, longp) \
2130 : (*(xdrs)->x_ops->x_getlong)(xdrs, longp)
2131 : #define xdr_putlong(xdrs, longp) \
2132 : (*(xdrs)->x_ops->x_putlong)(xdrs, longp)
2133 : #define xdr_getbytes(xdrs, addr, len) \
2134 : (*(xdrs)->x_ops->x_getbytes)(xdrs, addr, len)
2135 : #define xdr_putbytes(xdrs, addr, len) \
2136 : (*(xdrs)->x_ops->x_putbytes)(xdrs, addr, len)
2137 :
2138 : #define BYTES_PER_XDR_UNIT 4
2139 : static char xdr_zero[BYTES_PER_XDR_UNIT] = {0, 0, 0, 0};
2140 :
2141 : static int32_t
2142 26554 : xdr_swapbytes(int32_t x)
2143 : {
2144 : int32_t y,i;
2145 : char *px=(char *)&x;
2146 : char *py=(char *)&y;
2147 :
2148 132770 : for(i=0;i<4;i++)
2149 106216 : py[i]=px[3-i];
2150 :
2151 26554 : return y;
2152 : }
2153 :
2154 : static int32_t
2155 9972 : xdr_htonl(int32_t x)
2156 : {
2157 : int s=0x1234;
2158 : if( *((char *)&s)==(char)0x34) {
2159 : /* smallendian,swap bytes */
2160 9972 : return xdr_swapbytes(x);
2161 : } else {
2162 : /* bigendian, do nothing */
2163 : return x;
2164 : }
2165 : }
2166 :
2167 : static int32_t
2168 16582 : xdr_ntohl(int x)
2169 : {
2170 : int s=0x1234;
2171 : if( *((char *)&s)==(char)0x34) {
2172 : /* smallendian, swap bytes */
2173 16582 : return xdr_swapbytes(x);
2174 : } else {
2175 : /* bigendian, do nothing */
2176 : return x;
2177 : }
2178 : }
2179 :
2180 : static int
2181 2966 : xdr_int (XDR *xdrs, int *ip)
2182 : {
2183 : int32_t i32;
2184 :
2185 2966 : switch (xdrs->x_op)
2186 : {
2187 1488 : case XDR_ENCODE:
2188 1488 : i32 = (int32_t) *ip;
2189 1488 : return xdr_putlong (xdrs, &i32);
2190 :
2191 1478 : case XDR_DECODE:
2192 1478 : if (!xdr_getlong (xdrs, &i32))
2193 : {
2194 : return 0;
2195 : }
2196 1467 : *ip = (int) i32;
2197 : case XDR_FREE:
2198 : return 1;
2199 : }
2200 : return 0;
2201 : }
2202 :
2203 : static int
2204 123 : xdr_u_int (XDR *xdrs, unsigned int *up)
2205 : {
2206 : uint32_t ui32;
2207 :
2208 123 : switch (xdrs->x_op)
2209 : {
2210 48 : case XDR_ENCODE:
2211 48 : ui32 = (uint32_t) * up;
2212 48 : return xdr_putlong (xdrs, (int32_t *)&ui32);
2213 :
2214 75 : case XDR_DECODE:
2215 75 : if (!xdr_getlong (xdrs, (int32_t *)&ui32))
2216 : {
2217 : return 0;
2218 : }
2219 75 : *up = (uint32_t) ui32;
2220 : case XDR_FREE:
2221 : return 1;
2222 : }
2223 : return 0;
2224 : }
2225 :
2226 : static int
2227 0 : xdr_short (XDR *xdrs, short *sp)
2228 : {
2229 : int32_t i32;
2230 :
2231 0 : switch (xdrs->x_op)
2232 : {
2233 0 : case XDR_ENCODE:
2234 0 : i32 = (int32_t) *sp;
2235 0 : return xdr_putlong (xdrs, &i32);
2236 :
2237 0 : case XDR_DECODE:
2238 0 : if (!xdr_getlong (xdrs, &i32))
2239 : {
2240 : return 0;
2241 : }
2242 0 : *sp = (short) i32;
2243 0 : return 1;
2244 :
2245 : case XDR_FREE:
2246 : return 1;
2247 : }
2248 0 : return 0;
2249 : }
2250 :
2251 : static int
2252 0 : xdr_u_short (XDR *xdrs, unsigned short *sp)
2253 : {
2254 : uint32_t ui32;
2255 :
2256 0 : switch (xdrs->x_op)
2257 : {
2258 0 : case XDR_ENCODE:
2259 0 : ui32 = (uint32_t) *sp;
2260 0 : return xdr_putlong (xdrs, (int32_t *)&ui32);
2261 :
2262 0 : case XDR_DECODE:
2263 0 : if (!xdr_getlong (xdrs, (int32_t *)&ui32))
2264 : {
2265 : return 0;
2266 : }
2267 0 : *sp = (unsigned short) ui32;
2268 0 : return 1;
2269 :
2270 : case XDR_FREE:
2271 : return 1;
2272 : }
2273 0 : return 0;
2274 : }
2275 :
2276 : static int
2277 0 : xdr_char (XDR *xdrs, char *cp)
2278 : {
2279 : int i;
2280 :
2281 0 : i = (*cp);
2282 0 : if (!xdr_int (xdrs, &i))
2283 : {
2284 : return 0;
2285 : }
2286 0 : *cp = i;
2287 0 : return 1;
2288 : }
2289 :
2290 : static int
2291 0 : xdr_u_char (XDR *xdrs, unsigned char *cp)
2292 : {
2293 : unsigned int u;
2294 :
2295 0 : u = (*cp);
2296 0 : if (!xdr_u_int (xdrs, &u))
2297 : {
2298 : return 0;
2299 : }
2300 0 : *cp = u;
2301 0 : return 1;
2302 : }
2303 :
2304 : /*
2305 : * XDR opaque data
2306 : * Allows the specification of a fixed size sequence of opaque bytes.
2307 : * cp points to the opaque object and cnt gives the byte length.
2308 : */
2309 : static int
2310 211 : xdr_opaque (XDR *xdrs, char *cp, unsigned int cnt)
2311 : {
2312 : unsigned int rndup;
2313 : static char crud[BYTES_PER_XDR_UNIT];
2314 :
2315 : /*
2316 : * if no data we are done
2317 : */
2318 211 : if (cnt == 0)
2319 : return 1;
2320 :
2321 : /*
2322 : * round byte count to full xdr units
2323 : */
2324 211 : rndup = cnt % BYTES_PER_XDR_UNIT;
2325 211 : if (rndup > 0)
2326 61 : rndup = BYTES_PER_XDR_UNIT - rndup;
2327 :
2328 211 : switch (xdrs->x_op)
2329 : {
2330 103 : case XDR_DECODE:
2331 103 : if (!xdr_getbytes (xdrs, cp, cnt))
2332 : {
2333 : return 0;
2334 : }
2335 103 : if (rndup == 0)
2336 : return 1;
2337 22 : return xdr_getbytes (xdrs, (char *)crud, rndup);
2338 :
2339 108 : case XDR_ENCODE:
2340 108 : if (!xdr_putbytes (xdrs, cp, cnt))
2341 : {
2342 : return 0;
2343 : }
2344 108 : if (rndup == 0)
2345 : return 1;
2346 39 : return xdr_putbytes (xdrs, xdr_zero, rndup);
2347 :
2348 : case XDR_FREE:
2349 : return 1;
2350 : }
2351 : #undef BYTES_PER_XDR_UNIT
2352 0 : return 0;
2353 : }
2354 :
2355 :
2356 : /*
2357 : * XDR null terminated ASCII strings
2358 : */
2359 : static int
2360 123 : xdr_string (XDR *xdrs, char **cpp, unsigned int maxsize)
2361 : {
2362 123 : char *sp = *cpp; /* sp is the actual string pointer */
2363 : unsigned int size;
2364 : unsigned int nodesize;
2365 :
2366 : /*
2367 : * first deal with the length since xdr strings are counted-strings
2368 : */
2369 123 : switch (xdrs->x_op)
2370 : {
2371 0 : case XDR_FREE:
2372 0 : if (sp == NULL)
2373 : {
2374 : return 1; /* already free */
2375 : }
2376 : /* fall through... */
2377 : case XDR_ENCODE:
2378 48 : if (sp == NULL)
2379 : return 0;
2380 48 : size = strlen (sp);
2381 48 : break;
2382 : case XDR_DECODE:
2383 : break;
2384 : }
2385 123 : if (!xdr_u_int (xdrs, &size))
2386 : {
2387 : return 0;
2388 : }
2389 123 : if (size > maxsize)
2390 : {
2391 : return 0;
2392 : }
2393 123 : nodesize = size + 1;
2394 :
2395 : /*
2396 : * now deal with the actual bytes
2397 : */
2398 123 : switch (xdrs->x_op)
2399 : {
2400 75 : case XDR_DECODE:
2401 75 : if (nodesize == 0)
2402 : {
2403 : return 1;
2404 : }
2405 75 : if (sp == NULL)
2406 0 : *cpp = sp = (char *) malloc (nodesize);
2407 75 : if (sp == NULL)
2408 : {
2409 0 : (void) fputs ("xdr_string: out of memory\n", stderr);
2410 0 : return 0;
2411 : }
2412 75 : sp[size] = 0;
2413 : /* fall into ... */
2414 :
2415 123 : case XDR_ENCODE:
2416 123 : return xdr_opaque (xdrs, sp, size);
2417 :
2418 0 : case XDR_FREE:
2419 0 : free (sp);
2420 0 : *cpp = NULL;
2421 0 : return 1;
2422 : }
2423 : return 0;
2424 : }
2425 :
2426 :
2427 :
2428 : /* Floating-point stuff */
2429 :
2430 : static int
2431 23476 : xdr_float(XDR *xdrs, float *fp)
2432 : {
2433 23476 : switch (xdrs->x_op) {
2434 :
2435 : case XDR_ENCODE:
2436 : if (sizeof(float) == sizeof(int32_t))
2437 8436 : return (xdr_putlong(xdrs, (int32_t *) (void*) fp));
2438 : else if (sizeof(float) == sizeof(int)) {
2439 : int32_t tmp = *(int *) (void*) fp;
2440 : return (xdr_putlong(xdrs, &tmp));
2441 : }
2442 : break;
2443 :
2444 : case XDR_DECODE:
2445 : if (sizeof(float) == sizeof(int32_t))
2446 15040 : return (xdr_getlong(xdrs, (int32_t *) (void*) fp));
2447 : else if (sizeof(float) == sizeof(int)) {
2448 : int32_t tmp;
2449 : if (xdr_getlong(xdrs, &tmp)) {
2450 : *(int *) (void*)fp = tmp;
2451 : return (1);
2452 : }
2453 : }
2454 : break;
2455 :
2456 : case XDR_FREE:
2457 : return (1);
2458 : }
2459 0 : return (0);
2460 : }
2461 :
2462 :
2463 : static int
2464 0 : xdr_double(XDR *xdrs, double *dp)
2465 : {
2466 : /* Gromacs detects floating-point stuff at compile time, which is faster */
2467 : #ifdef GROMACS
2468 : # ifndef FLOAT_FORMAT_IEEE754
2469 : # error non-IEEE floating point system, or you defined GROMACS yourself...
2470 : # endif
2471 : int LSW;
2472 : # ifdef IEEE754_BIG_ENDIAN_WORD_ORDER
2473 : int LSW=1;
2474 : # else
2475 : int LSW=0;
2476 : # endif /* Big endian word order */
2477 : #else
2478 : /* Outside Gromacs we rely on dynamic detection of FP order. */
2479 : int LSW; /* Least significant fp word */
2480 :
2481 : double x=0.987654321; /* Just a number */
2482 : unsigned char ix = *((char *)&x);
2483 :
2484 : /* Possible representations in IEEE double precision:
2485 : * (S=small endian, B=big endian)
2486 : *
2487 : * Byte order, Word order, Hex
2488 : * S S b8 56 0e 3c dd 9a ef 3f
2489 : * B S 3c 0e 56 b8 3f ef 9a dd
2490 : * S B dd 9a ef 3f b8 56 0e 3c
2491 : * B B 3f ef 9a dd 3c 0e 56 b8
2492 : */
2493 : if(ix==0xdd || ix==0x3f)
2494 : LSW=1; /* Big endian word order */
2495 : else if(ix==0xb8 || ix==0x3c)
2496 : LSW=0; /* Small endian word order */
2497 : else { /* Catch strange errors */
2498 : fprintf(stderr,"Cannot detect floating-point word order.\n"
2499 : "Do you have a non-IEEE system?\n"
2500 : "Use system XDR libraries or fix xdr_double().\n");
2501 : abort();
2502 : }
2503 : #endif /* end of dynamic detection of fp word order */
2504 :
2505 0 : switch (xdrs->x_op) {
2506 :
2507 : case XDR_ENCODE:
2508 : if (2*sizeof(int32_t) == sizeof(double)) {
2509 : int32_t *lp = (int32_t *) (void*) dp;
2510 0 : return (xdr_putlong(xdrs, lp+!LSW) &&
2511 0 : xdr_putlong(xdrs, lp+LSW));
2512 : } else if (2*sizeof(int) == sizeof(double)) {
2513 : int *ip = (int *) (void*) dp;
2514 : int32_t tmp[2];
2515 : tmp[0] = ip[!LSW];
2516 : tmp[1] = ip[LSW];
2517 : return (xdr_putlong(xdrs, tmp) &&
2518 : xdr_putlong(xdrs, tmp+1));
2519 : }
2520 : break;
2521 :
2522 : case XDR_DECODE:
2523 : if (2*sizeof(int32_t) == sizeof(double)) {
2524 : int32_t *lp = (int32_t *) (void*) dp;
2525 0 : return (xdr_getlong(xdrs, lp+!LSW) &&
2526 0 : xdr_getlong(xdrs, lp+LSW));
2527 : } else if (2*sizeof(int) == sizeof(double)) {
2528 : int *ip = (int *) (void*) dp;
2529 : int32_t tmp[2];
2530 : if (xdr_getlong(xdrs, tmp+!LSW) &&
2531 : xdr_getlong(xdrs, tmp+LSW)) {
2532 : ip[0] = tmp[0];
2533 : ip[1] = tmp[1];
2534 : return (1);
2535 : }
2536 : }
2537 : break;
2538 :
2539 : case XDR_FREE:
2540 : return (1);
2541 : }
2542 0 : return (0);
2543 : }
2544 :
2545 :
2546 : static int xdrstdio_getlong (XDR *, int32_t *);
2547 : static int xdrstdio_putlong (XDR *, int32_t *);
2548 : static int xdrstdio_getbytes (XDR *, char *, unsigned int);
2549 : static int xdrstdio_putbytes (XDR *, char *, unsigned int);
2550 : static unsigned int xdrstdio_getpos (XDR *);
2551 : static int xdrstdio_setpos (XDR *, unsigned int);
2552 : static void xdrstdio_destroy (XDR *);
2553 :
2554 : /*
2555 : * Ops vector for stdio type XDR
2556 : */
2557 : static const XDR::xdr_ops xdrstdio_ops =
2558 : {
2559 : xdrstdio_getlong, /* deserialize a long int */
2560 : xdrstdio_putlong, /* serialize a long int */
2561 : xdrstdio_getbytes, /* deserialize counted bytes */
2562 : xdrstdio_putbytes, /* serialize counted bytes */
2563 : xdrstdio_getpos, /* get offset in the stream */
2564 : xdrstdio_setpos, /* set offset in the stream */
2565 : xdrstdio_destroy, /* destroy stream */
2566 : };
2567 :
2568 : /*
2569 : * Initialize a stdio xdr stream.
2570 : * Sets the xdr stream handle xdrs for use on the stream file.
2571 : * Operation flag is set to op.
2572 : */
2573 : static void
2574 32 : xdrstdio_create (XDR *xdrs, FILE *file, enum xdr_op op)
2575 : {
2576 32 : xdrs->x_op = op;
2577 :
2578 32 : xdrs->x_ops = (XDR::xdr_ops *) &xdrstdio_ops;
2579 32 : xdrs->x_private = (char *) file;
2580 32 : }
2581 :
2582 : /*
2583 : * Destroy a stdio xdr stream.
2584 : * Cleans up the xdr stream handle xdrs previously set up by xdrstdio_create.
2585 : */
2586 : static void
2587 32 : xdrstdio_destroy (XDR *xdrs)
2588 : {
2589 32 : (void) fflush ((FILE *) xdrs->x_private);
2590 : /* xx should we close the file ?? */
2591 32 : }
2592 :
2593 : static int
2594 16593 : xdrstdio_getlong (XDR *xdrs, int32_t *lp)
2595 : {
2596 : int32_t mycopy;
2597 :
2598 16593 : if (fread ((char *) & mycopy, 4, 1, (FILE *) xdrs->x_private) != 1)
2599 : return 0;
2600 16582 : *lp = (int32_t) xdr_ntohl (mycopy);
2601 16582 : return 1;
2602 : }
2603 :
2604 : static int
2605 9972 : xdrstdio_putlong (XDR *xdrs, int32_t *lp)
2606 : {
2607 9972 : int32_t mycopy = xdr_htonl (*lp);
2608 : lp = &mycopy;
2609 9972 : if (fwrite ((char *) lp, 4, 1, (FILE *) xdrs->x_private) != 1)
2610 0 : return 0;
2611 : return 1;
2612 : }
2613 :
2614 : static int
2615 125 : xdrstdio_getbytes (XDR *xdrs, char *addr, unsigned int len)
2616 : {
2617 250 : if ((len != 0) && (fread (addr, (int) len, 1,
2618 125 : (FILE *) xdrs->x_private) != 1))
2619 0 : return 0;
2620 : return 1;
2621 : }
2622 :
2623 : static int
2624 147 : xdrstdio_putbytes (XDR *xdrs, char *addr, unsigned int len)
2625 : {
2626 294 : if ((len != 0) && (fwrite (addr, (int) len, 1,
2627 147 : (FILE *) xdrs->x_private) != 1))
2628 0 : return 0;
2629 : return 1;
2630 : }
2631 :
2632 : /* 32 bit fileseek operations */
2633 : static unsigned int
2634 0 : xdrstdio_getpos (XDR *xdrs)
2635 : {
2636 0 : return (unsigned int) ftell ((FILE *) xdrs->x_private);
2637 : }
2638 :
2639 : static int
2640 0 : xdrstdio_setpos (XDR *xdrs, unsigned int pos)
2641 : {
2642 0 : return fseek ((FILE *) xdrs->x_private, pos, 0) < 0 ? 0 : 1;
2643 : }
2644 :
2645 : #endif /* HAVE_RPC_XDR_H not defined */
2646 :
2647 : }
2648 : }
|
