/* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*-
 *
 * $Id$
 *
 * Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, this
 * list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/* Get HAVE_RPC_XDR_H, F77_FUNC from config.h if available */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <limits.h>

#define _FILE_OFFSET_BITS 64

/* get fixed-width types if we are using ANSI C99 */
#ifdef HAVE_STDINT_H
#include <stdint.h>
#elif (defined HAVE_INTTYPES_H)
#include <inttypes.h>
#endif

#ifdef HAVE_RPC_XDR_H
#include <rpc/rpc.h>
#include <rpc/xdr.h>
#endif

#include "xdrfile.h"

/* Default FORTRAN name mangling is: lower case name, append underscore */
#ifndef F77_FUNC
#define F77_FUNC(name, NAME) name##_
#endif

char *exdr_message[exdrNR] = {
	"OK",
	"Header",
	"String",
	"Double",
	"Integer",
	"Float",
	"Unsigned integer",
	"Compressed 3D coordinate",
	"Closing file",
	"Magic number",
	"Not enough memory",
	"End of file",
	"File not found"};

/*
 * Declare our own XDR routines statically if no libraries are present.
 * Actual implementation is at the end of this file.
 *
 * We don't want the low-level XDR implementation as part of the Gromacs
 * documentation, so skip it for doxygen too...
 */
#if (!defined HAVE_RPC_XDR_H && !defined DOXYGEN)

enum xdr_op
{
	XDR_ENCODE = 0,
	XDR_DECODE = 1,
	XDR_FREE = 2
};

/* We need integer types that are guaranteed to be 4 bytes wide.
 * If ANSI C99 headers were included they are already defined
 * as int32_t and uint32_t. Check, and if not define them ourselves.
 * Since it is just our workaround for missing ANSI C99 types, avoid adding
 * it to the doxygen documentation.
 */
#if !(defined INT32_MAX || defined DOXYGEN)
#if (INT_MAX == 2147483647)
#define int32_t int
#define uint32_t unsigned int
#define INT32_MAX 2147483647
#elif (LONG_MAX == 2147483647)
#define int32_t long
#define uint32_t unsigned long
#define INT32_MAX 2147483647L
#else
#error ERROR: No 32 bit wide integer type found!
#error Use system XDR libraries instead, or update xdrfile.c
#endif
#endif

typedef struct XDR XDR;

struct XDR
{
	enum xdr_op x_op;
	struct xdr_ops
	{
		int (*x_getlong)(XDR *__xdrs, int32_t *__lp);
		int (*x_putlong)(XDR *__xdrs, int32_t *__lp);
		int (*x_getbytes)(XDR *__xdrs, char *__addr, unsigned int __len);
		int (*x_putbytes)(XDR *__xdrs, char *__addr, unsigned int __len);
		/* two next routines are not 64-bit IO safe - don't use! */
		unsigned int (*x_getpostn)(XDR *__xdrs);
		int (*x_setpostn)(XDR *__xdrs, unsigned int __pos);
		void (*x_destroy)(XDR *__xdrs);
	}
		*x_ops;
	char *x_private;
};

static int xdr_char(XDR *xdrs, char *ip);
static int xdr_u_char(XDR *xdrs, unsigned char *ip);
static int xdr_short(XDR *xdrs, short *ip);
static int xdr_u_short(XDR *xdrs, unsigned short *ip);
static int xdr_int(XDR *xdrs, int *ip);
static int xdr_u_int(XDR *xdrs, unsigned int *ip);
static int xdr_float(XDR *xdrs, float *ip);
static int xdr_double(XDR *xdrs, double *ip);
static int xdr_string(XDR *xdrs, char **ip, unsigned int maxsize);
static int xdr_opaque(XDR *xdrs, char *cp, unsigned int cnt);
static void xdrstdio_create(XDR *xdrs, FILE *fp, enum xdr_op xop);

#define xdr_getpos(xdrs) \
	(*(xdrs)->x_ops->x_getpostn)(xdrs)
#define xdr_setpos(xdrs, pos) \
	(*(xdrs)->x_ops->x_setpostn)(xdrs, pos)
#define xdr_destroy(xdrs)                      \
	do                                         \
	{                                          \
		if ((xdrs)->x_ops->x_destroy)          \
			(*(xdrs)->x_ops->x_destroy)(xdrs); \
	} while (0)
#endif /* end of our own XDR declarations */

/** Contents of the abstract XDRFILE data structure.
 *
 *  @internal
 *
 *  This structure is used to provide an XDR file interface that is
 *  virtual identical to the standard UNIX fopen/fread/fwrite/fclose.
 */
struct XDRFILE
{
	FILE *fp;	  /**< pointer to standard C library file handle */
	XDR *xdr;	  /**< pointer to corresponding XDR handle       */
	char mode;	  /**< r=read, w=write, a=append                 */
	int *buf1;	  /**< Buffer for internal use                   */
	int buf1size; /**< Current allocated length of buf1          */
	int *buf2;	  /**< Buffer for internal use                   */
	int buf2size; /**< Current allocated length of buf2          */
};

/*************************************************************
 * Implementation of higher-level routines to read/write     *
 * portable data based on the XDR standard. These should be  *
 * called from C - see further down for Fortran77 wrappers.  *
 *************************************************************/

XDRFILE *
xdrfile_open(const char *path, const char *mode)
{
	char newmode[5];
	enum xdr_op xdrmode;
	XDRFILE *xfp;

	/* make sure XDR files are opened in binary mode... */
	if (*mode == 'w' || *mode == 'W')
	{
		sprintf(newmode, "wb+");
		xdrmode = XDR_ENCODE;
	}
	else if (*mode == 'a' || *mode == 'A')
	{
		sprintf(newmode, "ab+");
		xdrmode = XDR_ENCODE;
	}
	else if (*mode == 'r' || *mode == 'R')
	{
		sprintf(newmode, "rb");
		xdrmode = XDR_DECODE;
	}
	else /* cannot determine mode */
		return NULL;

	if ((xfp = (XDRFILE *)malloc(sizeof(XDRFILE))) == NULL)
		return NULL;
	if ((xfp->fp = fopen(path, newmode)) == NULL)
	{
		free(xfp);
		return NULL;
	}
	if ((xfp->xdr = (XDR *)malloc(sizeof(XDR))) == NULL)
	{
		fclose(xfp->fp);
		free(xfp);
		return NULL;
	}
	xfp->mode = *mode;
	xdrstdio_create((XDR *)(xfp->xdr), xfp->fp, xdrmode);
	xfp->buf1 = xfp->buf2 = NULL;
	xfp->buf1size = xfp->buf2size = 0;
	return xfp;
}

int xdrfile_close(XDRFILE *xfp)
{
	int ret = exdrCLOSE;
	if (xfp)
	{
		/* flush and destroy XDR stream */
		if (xfp->xdr)
			xdr_destroy((XDR *)(xfp->xdr));
		free(xfp->xdr);
		/* close the file */
		ret = fclose(xfp->fp);
		if (xfp->buf1size)
			free(xfp->buf1);
		if (xfp->buf2size)
			free(xfp->buf2);
		free(xfp);
	}
	return ret; /* return 0 if ok */
}

int xdrfile_read_int(int *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_int((XDR *)(xfp->xdr), ptr + i))
		i++;

	return i;
}

int xdrfile_write_int(int *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_int((XDR *)(xfp->xdr), ptr + i))
		i++;

	return i;
}

int xdrfile_read_uint(unsigned int *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_u_int((XDR *)(xfp->xdr), ptr + i))
		i++;

	return i;
}

int xdrfile_write_uint(unsigned int *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_u_int((XDR *)(xfp->xdr), ptr + i))
		i++;
	return i;
}

int xdrfile_read_char(char *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_char((XDR *)(xfp->xdr), ptr + i))
		i++;

	return i;
}

int xdrfile_write_char(char *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_char((XDR *)(xfp->xdr), ptr + i))
		i++;
	return i;
}

int xdrfile_read_uchar(unsigned char *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_u_char((XDR *)(xfp->xdr), ptr + i))
		i++;

	return i;
}

int xdrfile_write_uchar(unsigned char *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_u_char((XDR *)(xfp->xdr), ptr + i))
		i++;
	return i;
}

int xdrfile_read_short(short *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_short((XDR *)(xfp->xdr), ptr + i))
		i++;

	return i;
}

int xdrfile_write_short(short *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_short((XDR *)(xfp->xdr), ptr + i))
		i++;
	return i;
}

int xdrfile_read_ushort(unsigned short *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_u_short((XDR *)(xfp->xdr), ptr + i))
		i++;

	return i;
}

int xdrfile_write_ushort(unsigned short *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;

	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_u_short((XDR *)(xfp->xdr), ptr + i))
		i++;
	return i;
}

int xdrfile_read_float(float *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;
	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_float((XDR *)(xfp->xdr), ptr + i))
		i++;
	return i;
}

int xdrfile_write_float(float *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;
	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_float((XDR *)(xfp->xdr), ptr + i))
		i++;
	return i;
}

int xdrfile_read_double(double *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;
	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_double((XDR *)(xfp->xdr), ptr + i))
		i++;
	return i;
}

int xdrfile_write_double(double *ptr, int ndata, XDRFILE *xfp)
{
	int i = 0;
	/* read write is encoded in the XDR struct */
	while (i < ndata && xdr_double((XDR *)(xfp->xdr), ptr + i))
		i++;
	return i;
}

int xdrfile_read_string(char *ptr, int maxlen, XDRFILE *xfp)
{
	int i;
	if (xdr_string((XDR *)(xfp->xdr), &ptr, maxlen))
	{
		i = 0;
		while (i < maxlen && ptr[i] != 0)
			i++;
		if (i == maxlen)
			return maxlen;
		else
			return i + 1;
	}
	else
		return 0;
}

int xdrfile_write_string(char *ptr, XDRFILE *xfp)
{
	int len = strlen(ptr) + 1;

	if (xdr_string((XDR *)(xfp->xdr), &ptr, len))
		return len;
	else
		return 0;
}

int xdrfile_read_opaque(char *ptr, int cnt, XDRFILE *xfp)
{
	if (xdr_opaque((XDR *)(xfp->xdr), ptr, cnt))
		return cnt;
	else
		return 0;
}

int xdrfile_write_opaque(char *ptr, int cnt, XDRFILE *xfp)
{
	if (xdr_opaque((XDR *)(xfp->xdr), ptr, cnt))
		return cnt;
	else
		return 0;
}

/* Internal support routines for reading/writing compressed coordinates
 * sizeofint - calculate smallest number of bits necessary
 * to represent a certain integer.
 */
static int
sizeofint(int size)
{
	unsigned int num = 1;
	int num_of_bits = 0;

	while (size >= num && num_of_bits < 32)
	{
		num_of_bits++;
		num <<= 1;
	}
	return num_of_bits;
}

/*
 * sizeofints - calculate 'bitsize' of compressed ints
 *
 * given a number of small unsigned integers and the maximum value
 * return the number of bits needed to read or write them with the
 * routines encodeints/decodeints. You need this parameter when
 * calling those routines.
 * (However, in some cases we can just use the variable 'smallidx'
 * which is the exact number of bits, and them we dont need to call
 * this routine).
 */
static int
sizeofints(int num_of_ints, unsigned int sizes[])
{
	int i, num;
	unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
	num_of_bytes = 1;
	bytes[0] = 1;
	num_of_bits = 0;
	for (i = 0; i < num_of_ints; i++)
	{
		tmp = 0;
		for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
		{
			tmp = bytes[bytecnt] * sizes[i] + tmp;
			bytes[bytecnt] = tmp & 0xff;
			tmp >>= 8;
		}
		while (tmp != 0)
		{
			bytes[bytecnt++] = tmp & 0xff;
			tmp >>= 8;
		}
		num_of_bytes = bytecnt;
	}
	num = 1;
	num_of_bytes--;
	while (bytes[num_of_bytes] >= num)
	{
		num_of_bits++;
		num *= 2;
	}
	return num_of_bits + num_of_bytes * 8;
}

/*
 * encodebits - encode num into buf using the specified number of bits
 *
 * This routines appends the value of num to the bits already present in
 * the array buf. You need to give it the number of bits to use and you had
 * better make sure that this number of bits is enough to hold the value.
 * Num must also be positive.
 */
static void
encodebits(int buf[], int num_of_bits, int num)
{

	unsigned int cnt, lastbyte;
	int lastbits;
	unsigned char *cbuf;

	cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
	cnt = (unsigned int)buf[0];
	lastbits = buf[1];
	lastbyte = (unsigned int)buf[2];
	while (num_of_bits >= 8)
	{
		lastbyte = (lastbyte << 8) | ((num >> (num_of_bits - 8)) /* & 0xff*/);
		cbuf[cnt++] = lastbyte >> lastbits;
		num_of_bits -= 8;
	}
	if (num_of_bits > 0)
	{
		lastbyte = (lastbyte << num_of_bits) | num;
		lastbits += num_of_bits;
		if (lastbits >= 8)
		{
			lastbits -= 8;
			cbuf[cnt++] = lastbyte >> lastbits;
		}
	}
	buf[0] = cnt;
	buf[1] = lastbits;
	buf[2] = lastbyte;
	if (lastbits > 0)
	{
		cbuf[cnt] = lastbyte << (8 - lastbits);
	}
}

/*
 * encodeints - encode a small set of small integers in compressed format
 *
 * this routine is used internally by xdr3dfcoord, to encode a set of
 * small integers to the buffer for writing to a file.
 * Multiplication with fixed (specified maximum) sizes is used to get
 * to one big, multibyte integer. Allthough the routine could be
 * modified to handle sizes bigger than 16777216, or more than just
 * a few integers, this is not done because the gain in compression
 * isn't worth the effort. Note that overflowing the multiplication
 * or the byte buffer (32 bytes) is unchecked and whould cause bad results.
 * THese things are checked in the calling routines, so make sure not
 * to remove those checks...
 */

static void
encodeints(int buf[], int num_of_ints, int num_of_bits,
		   unsigned int sizes[], unsigned int nums[])
{

	int i;
	unsigned int bytes[32], num_of_bytes, bytecnt, tmp;

	tmp = nums[0];
	num_of_bytes = 0;
	do
	{
		bytes[num_of_bytes++] = tmp & 0xff;
		tmp >>= 8;
	} while (tmp != 0);

	for (i = 1; i < num_of_ints; i++)
	{
		if (nums[i] >= sizes[i])
		{
			fprintf(stderr, "major breakdown in encodeints - num %u doesn't "
							"match size %u\n",
					nums[i], sizes[i]);
			abort();
		}
		/* use one step multiply */
		tmp = nums[i];
		for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
		{
			tmp = bytes[bytecnt] * sizes[i] + tmp;
			bytes[bytecnt] = tmp & 0xff;
			tmp >>= 8;
		}
		while (tmp != 0)
		{
			bytes[bytecnt++] = tmp & 0xff;
			tmp >>= 8;
		}
		num_of_bytes = bytecnt;
	}
	if (num_of_bits >= num_of_bytes * 8)
	{
		for (i = 0; i < num_of_bytes; i++)
		{
			encodebits(buf, 8, bytes[i]);
		}
		encodebits(buf, num_of_bits - num_of_bytes * 8, 0);
	}
	else
	{
		for (i = 0; i < num_of_bytes - 1; i++)
		{
			encodebits(buf, 8, bytes[i]);
		}
		encodebits(buf, num_of_bits - (num_of_bytes - 1) * 8, bytes[i]);
	}
}

/*
 * decodebits - decode number from buf using specified number of bits
 *
 * extract the number of bits from the array buf and construct an integer
 * from it. Return that value.
 *
 */

static int
decodebits(int buf[], int num_of_bits)
{

	int cnt, num;
	unsigned int lastbits, lastbyte;
	unsigned char *cbuf;
	int mask = (1 << num_of_bits) - 1;

	cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
	cnt = buf[0];
	lastbits = (unsigned int)buf[1];
	lastbyte = (unsigned int)buf[2];

	num = 0;
	while (num_of_bits >= 8)
	{
		lastbyte = (lastbyte << 8) | cbuf[cnt++];
		num |= (lastbyte >> lastbits) << (num_of_bits - 8);
		num_of_bits -= 8;
	}
	if (num_of_bits > 0)
	{
		if (lastbits < num_of_bits)
		{
			lastbits += 8;
			lastbyte = (lastbyte << 8) | cbuf[cnt++];
		}
		lastbits -= num_of_bits;
		num |= (lastbyte >> lastbits) & ((1 << num_of_bits) - 1);
	}
	num &= mask;
	buf[0] = cnt;
	buf[1] = lastbits;
	buf[2] = lastbyte;
	return num;
}

/*
 * decodeints - decode 'small' integers from the buf array
 *
 * this routine is the inverse from encodeints() and decodes the small integers
 * written to buf by calculating the remainder and doing divisions with
 * the given sizes[]. You need to specify the total number of bits to be
 * used from buf in num_of_bits.
 *
 */

static void
decodeints(int buf[], int num_of_ints, int num_of_bits,
		   unsigned int sizes[], int nums[])
{

	int bytes[32];
	int i, j, num_of_bytes, p, num;

	bytes[1] = bytes[2] = bytes[3] = 0;
	num_of_bytes = 0;
	while (num_of_bits > 8)
	{
		bytes[num_of_bytes++] = decodebits(buf, 8);
		num_of_bits -= 8;
	}
	if (num_of_bits > 0)
	{
		bytes[num_of_bytes++] = decodebits(buf, num_of_bits);
	}
	for (i = num_of_ints - 1; i > 0; i--)
	{
		num = 0;
		for (j = num_of_bytes - 1; j >= 0; j--)
		{
			num = (num << 8) | bytes[j];
			p = num / sizes[i];
			bytes[j] = p;
			num = num - p * sizes[i];
		}
		nums[i] = num;
	}
	nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
}

static const int magicints[] =
	{
		0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
		80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 812, 1024, 1290,
		1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192, 10321, 13003,
		16384, 20642, 26007, 32768, 41285, 52015, 65536, 82570, 104031,
		131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
		832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021,
		4194304, 5284491, 6658042, 8388607, 10568983, 13316085, 16777216};

#define FIRSTIDX 9
/* note that magicints[FIRSTIDX-1] == 0 */
#define LASTIDX (sizeof(magicints) / sizeof(*magicints))

/* Compressed coordinate routines - modified from the original
 * implementation by Frans v. Hoesel to make them threadsafe.
 */
int xdrfile_decompress_coord_float(float *ptr,
								   int *size,
								   float *precision,
								   XDRFILE *xfp)
{
	int minint[3], maxint[3], *lip;
	int smallidx = 0;
	unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3;
	int k, *buf1, *buf2, flag;
	int lsize = 0;
	int smallnum, smaller, i, is_smaller, run;
	float *lfp, inv_precision;
	int tmp, *thiscoord, prevcoord[3];
	unsigned int bitsize;

	bitsizeint[0] = 0;
	bitsizeint[1] = 0;
	bitsizeint[2] = 0;

	if (xfp == NULL || ptr == NULL)
		return -1;
	tmp = xdrfile_read_int(&lsize, 1, xfp);
	if (tmp == 0)
		return -1; /* return if we could not read size */
	if (*size < lsize)
	{
		fprintf(stderr, "Requested to decompress %d coords, file contains %d\n",
				*size, lsize);
		return -1;
	}
	*size = lsize;
	size3 = *size * 3;
	if (size3 > xfp->buf1size)
	{
		if ((xfp->buf1 = (int *)malloc(sizeof(int) * size3)) == NULL)
		{
			fprintf(stderr, "Cannot allocate memory for decompressing coordinates.\n");
			return -1;
		}
		xfp->buf1size = size3;
		xfp->buf2size = size3 * 1.2;
		if ((xfp->buf2 = (int *)malloc(sizeof(int) * xfp->buf2size)) == NULL)
		{
			fprintf(stderr, "Cannot allocate memory for decompressing coordinates.\n");
			return -1;
		}
	}
	/* Dont bother with compression for three atoms or less */
	if (*size <= 9)
	{
		return xdrfile_read_float(ptr, size3, xfp) / 3;
		/* return number of coords, not floats */
	}
	/* Compression-time if we got here. Read precision first */
	xdrfile_read_float(precision, 1, xfp);

	/* avoid repeated pointer dereferencing. */
	buf1 = xfp->buf1;
	buf2 = xfp->buf2;
	/* buf2[0-2] are special and do not contain actual data */
	buf2[0] = buf2[1] = buf2[2] = 0;
	xdrfile_read_int(minint, 3, xfp);
	xdrfile_read_int(maxint, 3, xfp);

	sizeint[0] = maxint[0] - minint[0] + 1;
	sizeint[1] = maxint[1] - minint[1] + 1;
	sizeint[2] = maxint[2] - minint[2] + 1;

	/* check if one of the sizes is to big to be multiplied */
	if ((sizeint[0] | sizeint[1] | sizeint[2]) > 0xffffff)
	{
		bitsizeint[0] = sizeofint(sizeint[0]);
		bitsizeint[1] = sizeofint(sizeint[1]);
		bitsizeint[2] = sizeofint(sizeint[2]);
		bitsize = 0; /* flag the use of large sizes */
	}
	else
	{
		bitsize = sizeofints(3, sizeint);
	}

	if (xdrfile_read_int(&smallidx, 1, xfp) == 0)
		return 0; /* not sure what has happened here or why we return... */
	tmp = smallidx + 8;
	// maxidx = (LASTIDX<tmp) ? LASTIDX : tmp;
	// minidx = maxidx - 8; /* often this equal smallidx */
	tmp = smallidx - 1;
	tmp = (FIRSTIDX > tmp) ? FIRSTIDX : tmp;
	smaller = magicints[tmp] / 2;
	smallnum = magicints[smallidx] / 2;
	sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
	// larger = magicints[maxidx];

	/* buf2[0] holds the length in bytes */

	if (xdrfile_read_int(buf2, 1, xfp) == 0)
		return 0;
	if (xdrfile_read_opaque((char *)&(buf2[3]), (unsigned int)buf2[0], xfp) == 0)
		return 0;
	buf2[0] = buf2[1] = buf2[2] = 0;

	lfp = ptr;
	inv_precision = 1.0 / *precision;
	run = 0;
	i = 0;
	lip = buf1;
	while (i < lsize)
	{
		thiscoord = (int *)(lip) + i * 3;

		if (bitsize == 0)
		{
			thiscoord[0] = decodebits(buf2, bitsizeint[0]);
			thiscoord[1] = decodebits(buf2, bitsizeint[1]);
			thiscoord[2] = decodebits(buf2, bitsizeint[2]);
		}
		else
		{
			decodeints(buf2, 3, bitsize, sizeint, thiscoord);
		}

		i++;
		thiscoord[0] += minint[0];
		thiscoord[1] += minint[1];
		thiscoord[2] += minint[2];

		prevcoord[0] = thiscoord[0];
		prevcoord[1] = thiscoord[1];
		prevcoord[2] = thiscoord[2];

		flag = decodebits(buf2, 1);
		is_smaller = 0;
		if (flag == 1)
		{
			run = decodebits(buf2, 5);
			is_smaller = run % 3;
			run -= is_smaller;
			is_smaller--;
		}
		if (run > 0)
		{
			thiscoord += 3;
			for (k = 0; k < run; k += 3)
			{
				decodeints(buf2, 3, smallidx, sizesmall, thiscoord);
				i++;
				thiscoord[0] += prevcoord[0] - smallnum;
				thiscoord[1] += prevcoord[1] - smallnum;
				thiscoord[2] += prevcoord[2] - smallnum;
				if (k == 0)
				{
					/* interchange first with second atom for better
					 * compression of water molecules
					 */
					tmp = thiscoord[0];
					thiscoord[0] = prevcoord[0];
					prevcoord[0] = tmp;
					tmp = thiscoord[1];
					thiscoord[1] = prevcoord[1];
					prevcoord[1] = tmp;
					tmp = thiscoord[2];
					thiscoord[2] = prevcoord[2];
					prevcoord[2] = tmp;
					*lfp++ = prevcoord[0] * inv_precision;
					*lfp++ = prevcoord[1] * inv_precision;
					*lfp++ = prevcoord[2] * inv_precision;
				}
				else
				{
					prevcoord[0] = thiscoord[0];
					prevcoord[1] = thiscoord[1];
					prevcoord[2] = thiscoord[2];
				}
				*lfp++ = thiscoord[0] * inv_precision;
				*lfp++ = thiscoord[1] * inv_precision;
				*lfp++ = thiscoord[2] * inv_precision;
			}
		}
		else
		{
			*lfp++ = thiscoord[0] * inv_precision;
			*lfp++ = thiscoord[1] * inv_precision;
			*lfp++ = thiscoord[2] * inv_precision;
		}
		smallidx += is_smaller;
		if (is_smaller < 0)
		{
			smallnum = smaller;

			if (smallidx > FIRSTIDX)
			{
				smaller = magicints[smallidx - 1] / 2;
			}
			else
			{
				smaller = 0;
			}
		}
		else if (is_smaller > 0)
		{
			smaller = smallnum;
			smallnum = magicints[smallidx] / 2;
		}
		sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
	}
	return *size;
}

int xdrfile_compress_coord_float(float *ptr,
								 int size,
								 float precision,
								 XDRFILE *xfp)
{
	int minint[3], maxint[3], mindiff, *lip, diff;
	int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
	int minidx, maxidx;
	unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
	int k, *buf1, *buf2;
	int smallnum, smaller, larger, i, j, is_small, is_smaller, run, prevrun;
	float *lfp, lf;
	int tmp, tmpsum, *thiscoord;
	int prevcoord[3] = {
		0,
		0,
	};
	unsigned int tmpcoord[30];
	// int errval=1;
	unsigned int bitsize;

	if (xfp == NULL)
		return -1;
	size3 = 3 * size;

	bitsizeint[0] = 0;
	bitsizeint[1] = 0;
	bitsizeint[2] = 0;

	if (size3 > xfp->buf1size)
	{
		if ((xfp->buf1 = (int *)malloc(sizeof(int) * size3)) == NULL)
		{
			fprintf(stderr, "Cannot allocate memory for compressing coordinates.\n");
			return -1;
		}
		xfp->buf1size = size3;
		xfp->buf2size = size3 * 1.2;
		if ((xfp->buf2 = (int *)malloc(sizeof(int) * xfp->buf2size)) == NULL)
		{
			fprintf(stderr, "Cannot allocate memory for compressing coordinates.\n");
			return -1;
		}
	}
	if (xdrfile_write_int(&size, 1, xfp) == 0)
		return -1; /* return if we could not write size */
	/* Dont bother with compression for three atoms or less */
	if (size <= 9)
	{
		return xdrfile_write_float(ptr, size3, xfp) / 3;
		/* return number of coords, not floats */
	}
	/* Compression-time if we got here. Write precision first */
	if (precision <= 0)
		precision = 1000;
	xdrfile_write_float(&precision, 1, xfp);
	/* avoid repeated pointer dereferencing. */
	buf1 = xfp->buf1;
	buf2 = xfp->buf2;
	/* buf2[0-2] are special and do not contain actual data */
	buf2[0] = buf2[1] = buf2[2] = 0;
	minint[0] = minint[1] = minint[2] = INT_MAX;
	maxint[0] = maxint[1] = maxint[2] = INT_MIN;
	prevrun = -1;
	lfp = ptr;
	lip = buf1;
	mindiff = INT_MAX;
	oldlint1 = oldlint2 = oldlint3 = 0;
	while (lfp < ptr + size3)
	{
		/* find nearest integer */
		if (*lfp >= 0.0)
			lf = *lfp * precision + 0.5;
		else
			lf = *lfp * precision - 0.5;
		if (fabs(lf) > INT_MAX - 2)
		{
			/* scaling would cause overflow */
			fprintf(stderr, "Internal overflow compressing coordinates.\n");
			// errval=0;
		}
		lint1 = lf;
		if (lint1 < minint[0])
			minint[0] = lint1;
		if (lint1 > maxint[0])
			maxint[0] = lint1;
		*lip++ = lint1;
		lfp++;
		if (*lfp >= 0.0)
			lf = *lfp * precision + 0.5;
		else
			lf = *lfp * precision - 0.5;
		if (fabs(lf) > INT_MAX - 2)
		{
			/* scaling would cause overflow */
			fprintf(stderr, "Internal overflow compressing coordinates.\n");
			// errval=0;
		}
		lint2 = lf;
		if (lint2 < minint[1])
			minint[1] = lint2;
		if (lint2 > maxint[1])
			maxint[1] = lint2;
		*lip++ = lint2;
		lfp++;
		if (*lfp >= 0.0)
			lf = *lfp * precision + 0.5;
		else
			lf = *lfp * precision - 0.5;
		if (fabs(lf) > INT_MAX - 2)
		{
			// errval=0;
		}
		lint3 = lf;
		if (lint3 < minint[2])
			minint[2] = lint3;
		if (lint3 > maxint[2])
			maxint[2] = lint3;
		*lip++ = lint3;
		lfp++;
		diff = abs(oldlint1 - lint1) + abs(oldlint2 - lint2) + abs(oldlint3 - lint3);
		if (diff < mindiff && lfp > ptr + 3)
			mindiff = diff;
		oldlint1 = lint1;
		oldlint2 = lint2;
		oldlint3 = lint3;
	}
	xdrfile_write_int(minint, 3, xfp);
	xdrfile_write_int(maxint, 3, xfp);

	if ((float)maxint[0] - (float)minint[0] >= INT_MAX - 2 ||
		(float)maxint[1] - (float)minint[1] >= INT_MAX - 2 ||
		(float)maxint[2] - (float)minint[2] >= INT_MAX - 2)
	{
		/* turning value in unsigned by subtracting minint
		 * would cause overflow
		 */
		fprintf(stderr, "Internal overflow compressing coordinates.\n");
		// errval=0;
	}
	sizeint[0] = maxint[0] - minint[0] + 1;
	sizeint[1] = maxint[1] - minint[1] + 1;
	sizeint[2] = maxint[2] - minint[2] + 1;

	/* check if one of the sizes is to big to be multiplied */
	if ((sizeint[0] | sizeint[1] | sizeint[2]) > 0xffffff)
	{
		bitsizeint[0] = sizeofint(sizeint[0]);
		bitsizeint[1] = sizeofint(sizeint[1]);
		bitsizeint[2] = sizeofint(sizeint[2]);
		bitsize = 0; /* flag the use of large sizes */
	}
	else
	{
		bitsize = sizeofints(3, sizeint);
	}
	lip = buf1;
	luip = (unsigned int *)buf1;
	smallidx = FIRSTIDX;
	while (smallidx < LASTIDX && magicints[smallidx] < mindiff)
	{
		smallidx++;
	}
	xdrfile_write_int(&smallidx, 1, xfp);
	tmp = smallidx + 8;
	maxidx = (LASTIDX < tmp) ? LASTIDX : tmp;
	minidx = maxidx - 8; /* often this equal smallidx */
	tmp = smallidx - 1;
	tmp = (FIRSTIDX > tmp) ? FIRSTIDX : tmp;
	smaller = magicints[tmp] / 2;
	smallnum = magicints[smallidx] / 2;
	sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
	larger = magicints[maxidx] / 2;
	i = 0;
	while (i < size)
	{
		is_small = 0;
		thiscoord = (int *)(luip) + i * 3;
		if (smallidx < maxidx && i >= 1 &&
			abs(thiscoord[0] - prevcoord[0]) < larger &&
			abs(thiscoord[1] - prevcoord[1]) < larger &&
			abs(thiscoord[2] - prevcoord[2]) < larger)
		{
			is_smaller = 1;
		}
		else if (smallidx > minidx)
		{
			is_smaller = -1;
		}
		else
		{
			is_smaller = 0;
		}
		if (i + 1 < size)
		{
			if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
				abs(thiscoord[1] - thiscoord[4]) < smallnum &&
				abs(thiscoord[2] - thiscoord[5]) < smallnum)
			{
				/* interchange first with second atom for better
				 * compression of water molecules
				 */
				tmp = thiscoord[0];
				thiscoord[0] = thiscoord[3];
				thiscoord[3] = tmp;
				tmp = thiscoord[1];
				thiscoord[1] = thiscoord[4];
				thiscoord[4] = tmp;
				tmp = thiscoord[2];
				thiscoord[2] = thiscoord[5];
				thiscoord[5] = tmp;
				is_small = 1;
			}
		}
		tmpcoord[0] = thiscoord[0] - minint[0];
		tmpcoord[1] = thiscoord[1] - minint[1];
		tmpcoord[2] = thiscoord[2] - minint[2];
		if (bitsize == 0)
		{
			encodebits(buf2, bitsizeint[0], tmpcoord[0]);
			encodebits(buf2, bitsizeint[1], tmpcoord[1]);
			encodebits(buf2, bitsizeint[2], tmpcoord[2]);
		}
		else
		{
			encodeints(buf2, 3, bitsize, sizeint, tmpcoord);
		}
		prevcoord[0] = thiscoord[0];
		prevcoord[1] = thiscoord[1];
		prevcoord[2] = thiscoord[2];
		thiscoord = thiscoord + 3;
		i++;

		run = 0;
		if (is_small == 0 && is_smaller == -1)
			is_smaller = 0;
		while (is_small && run < 8 * 3)
		{
			tmpsum = 0;
			for (j = 0; j < 3; j++)
			{
				tmp = thiscoord[j] - prevcoord[j];
				tmpsum += tmp * tmp;
			}
			if (is_smaller == -1 && tmpsum >= smaller * smaller)
			{
				is_smaller = 0;
			}

			tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
			tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
			tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;

			prevcoord[0] = thiscoord[0];
			prevcoord[1] = thiscoord[1];
			prevcoord[2] = thiscoord[2];

			i++;
			thiscoord = thiscoord + 3;
			is_small = 0;
			if (i < size &&
				abs(thiscoord[0] - prevcoord[0]) < smallnum &&
				abs(thiscoord[1] - prevcoord[1]) < smallnum &&
				abs(thiscoord[2] - prevcoord[2]) < smallnum)
			{
				is_small = 1;
			}
		}
		if (run != prevrun || is_smaller != 0)
		{
			prevrun = run;
			encodebits(buf2, 1, 1); /* flag the change in run-length */
			encodebits(buf2, 5, run + is_smaller + 1);
		}
		else
		{
			encodebits(buf2, 1, 0); /* flag the fact that runlength did not change */
		}
		for (k = 0; k < run; k += 3)
		{
			encodeints(buf2, 3, smallidx, sizesmall, &tmpcoord[k]);
		}
		if (is_smaller != 0)
		{
			smallidx += is_smaller;
			if (is_smaller < 0)
			{
				smallnum = smaller;
				smaller = magicints[smallidx - 1] / 2;
			}
			else
			{
				smaller = smallnum;
				smallnum = magicints[smallidx] / 2;
			}
			sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
		}
	}
	if (buf2[1] != 0)
		buf2[0]++;
	xdrfile_write_int(buf2, 1, xfp); /* buf2[0] holds the length in bytes */
	tmp = xdrfile_write_opaque((char *)&(buf2[3]), (unsigned int)buf2[0], xfp);
	if (tmp == (unsigned int)buf2[0])
		return size;
	else
		return -1;
}

int xdrfile_decompress_coord_double(double *ptr,
									int *size,
									double *precision,
									XDRFILE *xfp)
{
	int minint[3], maxint[3], *lip;
	int smallidx = 0;
	unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3;
	int k, *buf1, *buf2, flag;
	int lsize = 0;
	int smallnum, smaller, i, is_smaller, run;
	double *lfp, inv_precision;
	float float_prec, tmpdata[30];
	int tmp, *thiscoord, prevcoord[3];
	unsigned int bitsize;

	bitsizeint[0] = 0;
	bitsizeint[1] = 0;
	bitsizeint[2] = 0;

	if (xfp == NULL || ptr == NULL)
		return -1;
	tmp = xdrfile_read_int(&lsize, 1, xfp);
	if (tmp == 0)
		return -1; /* return if we could not read size */
	if (*size < lsize)
	{
		fprintf(stderr, "Requested to decompress %d coords, file contains %d\n",
				*size, lsize);
		return -1;
	}
	*size = lsize;
	size3 = *size * 3;
	if (size3 > xfp->buf1size)
	{
		if ((xfp->buf1 = (int *)malloc(sizeof(int) * size3)) == NULL)
		{
			fprintf(stderr, "Cannot allocate memory for decompression coordinates.\n");
			return -1;
		}
		xfp->buf1size = size3;
		xfp->buf2size = size3 * 1.2;
		if ((xfp->buf2 = (int *)malloc(sizeof(int) * xfp->buf2size)) == NULL)
		{
			fprintf(stderr, "Cannot allocate memory for decompressing coordinates.\n");
			return -1;
		}
	}
	/* Dont bother with compression for three atoms or less */
	if (*size <= 9)
	{
		tmp = xdrfile_read_float(tmpdata, size3, xfp);
		for (i = 0; i < 9 * 3; i++)
			ptr[i] = tmpdata[i];
		return tmp / 3;
		/* return number of coords, not floats */
	}
	/* Compression-time if we got here. Read precision first */
	xdrfile_read_float(&float_prec, 1, xfp);
	*precision = float_prec;
	/* avoid repeated pointer dereferencing. */
	buf1 = xfp->buf1;
	buf2 = xfp->buf2;
	/* buf2[0-2] are special and do not contain actual data */
	buf2[0] = buf2[1] = buf2[2] = 0;
	xdrfile_read_int(minint, 3, xfp);
	xdrfile_read_int(maxint, 3, xfp);

	sizeint[0] = maxint[0] - minint[0] + 1;
	sizeint[1] = maxint[1] - minint[1] + 1;
	sizeint[2] = maxint[2] - minint[2] + 1;

	/* check if one of the sizes is to big to be multiplied */
	if ((sizeint[0] | sizeint[1] | sizeint[2]) > 0xffffff)
	{
		bitsizeint[0] = sizeofint(sizeint[0]);
		bitsizeint[1] = sizeofint(sizeint[1]);
		bitsizeint[2] = sizeofint(sizeint[2]);
		bitsize = 0; /* flag the use of large sizes */
	}
	else
	{
		bitsize = sizeofints(3, sizeint);
	}

	if (xdrfile_read_int(&smallidx, 1, xfp) == 0)
		return 0;
	tmp = smallidx + 8;
	// maxidx = (LASTIDX<tmp) ? LASTIDX : tmp;
	// minidx = maxidx - 8; /* often this equal smallidx */
	tmp = smallidx - 1;
	tmp = (FIRSTIDX > tmp) ? FIRSTIDX : tmp;
	smaller = magicints[tmp] / 2;
	smallnum = magicints[smallidx] / 2;
	sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
	// larger = magicints[maxidx];

	/* buf2[0] holds the length in bytes */

	if (xdrfile_read_int(buf2, 1, xfp) == 0)
		return 0;
	if (xdrfile_read_opaque((char *)&(buf2[3]), (unsigned int)buf2[0], xfp) == 0)
		return 0;
	buf2[0] = buf2[1] = buf2[2] = 0;

	lfp = ptr;
	inv_precision = 1.0 / *precision;
	run = 0;
	i = 0;
	lip = buf1;
	while (i < lsize)
	{
		thiscoord = (int *)(lip) + i * 3;

		if (bitsize == 0)
		{
			thiscoord[0] = decodebits(buf2, bitsizeint[0]);
			thiscoord[1] = decodebits(buf2, bitsizeint[1]);
			thiscoord[2] = decodebits(buf2, bitsizeint[2]);
		}
		else
		{
			decodeints(buf2, 3, bitsize, sizeint, thiscoord);
		}

		i++;
		thiscoord[0] += minint[0];
		thiscoord[1] += minint[1];
		thiscoord[2] += minint[2];

		prevcoord[0] = thiscoord[0];
		prevcoord[1] = thiscoord[1];
		prevcoord[2] = thiscoord[2];

		flag = decodebits(buf2, 1);
		is_smaller = 0;
		if (flag == 1)
		{
			run = decodebits(buf2, 5);
			is_smaller = run % 3;
			run -= is_smaller;
			is_smaller--;
		}
		if (run > 0)
		{
			thiscoord += 3;
			for (k = 0; k < run; k += 3)
			{
				decodeints(buf2, 3, smallidx, sizesmall, thiscoord);
				i++;
				thiscoord[0] += prevcoord[0] - smallnum;
				thiscoord[1] += prevcoord[1] - smallnum;
				thiscoord[2] += prevcoord[2] - smallnum;
				if (k == 0)
				{
					/* interchange first with second atom for better
					 * compression of water molecules
					 */
					tmp = thiscoord[0];
					thiscoord[0] = prevcoord[0];
					prevcoord[0] = tmp;
					tmp = thiscoord[1];
					thiscoord[1] = prevcoord[1];
					prevcoord[1] = tmp;
					tmp = thiscoord[2];
					thiscoord[2] = prevcoord[2];
					prevcoord[2] = tmp;
					*lfp++ = prevcoord[0] * inv_precision;
					*lfp++ = prevcoord[1] * inv_precision;
					*lfp++ = prevcoord[2] * inv_precision;
				}
				else
				{
					prevcoord[0] = thiscoord[0];
					prevcoord[1] = thiscoord[1];
					prevcoord[2] = thiscoord[2];
				}
				*lfp++ = thiscoord[0] * inv_precision;
				*lfp++ = thiscoord[1] * inv_precision;
				*lfp++ = thiscoord[2] * inv_precision;
			}
		}
		else
		{
			*lfp++ = thiscoord[0] * inv_precision;
			*lfp++ = thiscoord[1] * inv_precision;
			*lfp++ = thiscoord[2] * inv_precision;
		}
		smallidx += is_smaller;
		if (is_smaller < 0)
		{
			smallnum = smaller;
			if (smallidx > FIRSTIDX)
			{
				smaller = magicints[smallidx - 1] / 2;
			}
			else
			{
				smaller = 0;
			}
		}
		else if (is_smaller > 0)
		{
			smaller = smallnum;
			smallnum = magicints[smallidx] / 2;
		}
		sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
	}
	return *size;
}

int xdrfile_compress_coord_double(double *ptr,
								  int size,
								  double precision,
								  XDRFILE *xfp)
{
	int minint[3], maxint[3], mindiff, *lip, diff;
	int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
	int minidx, maxidx;
	unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
	int k, *buf1, *buf2;
	int smallnum, smaller, larger, i, j, is_small, is_smaller, run, prevrun;
	double *lfp;
	float float_prec, lf, tmpdata[30];
	int tmp, tmpsum, *thiscoord;
	int prevcoord[3] = {0, 0, 0};
	unsigned int tmpcoord[30];
	// int errval=1;
	unsigned int bitsize;

	bitsizeint[0] = 0;
	bitsizeint[1] = 0;
	bitsizeint[2] = 0;

	if (xfp == NULL)
		return -1;
	size3 = 3 * size;
	if (size3 > xfp->buf1size)
	{
		if ((xfp->buf1 = (int *)malloc(sizeof(int) * size3)) == NULL)
		{
			fprintf(stderr, "Cannot allocate memory for compressing coordinates.\n");
			return -1;
		}
		xfp->buf1size = size3;
		xfp->buf2size = size3 * 1.2;
		if ((xfp->buf2 = (int *)malloc(sizeof(int) * xfp->buf2size)) == NULL)
		{
			fprintf(stderr, "Cannot allocate memory for compressing coordinates.\n");
			return -1;
		}
	}
	if (xdrfile_write_int(&size, 1, xfp) == 0)
		return -1; /* return if we could not write size */
	/* Dont bother with compression for three atoms or less */
	if (size <= 9)
	{
		for (i = 0; i < 9 * 3; i++)
			tmpdata[i] = ptr[i];
		return xdrfile_write_float(tmpdata, size3, xfp) / 3;
		/* return number of coords, not floats */
	}
	/* Compression-time if we got here. Write precision first */
	if (precision <= 0)
		precision = 1000;
	float_prec = precision;
	xdrfile_write_float(&float_prec, 1, xfp);
	/* avoid repeated pointer dereferencing. */
	buf1 = xfp->buf1;
	buf2 = xfp->buf2;
	/* buf2[0-2] are special and do not contain actual data */
	buf2[0] = buf2[1] = buf2[2] = 0;
	minint[0] = minint[1] = minint[2] = INT_MAX;
	maxint[0] = maxint[1] = maxint[2] = INT_MIN;
	prevrun = -1;
	lfp = ptr;
	lip = buf1;
	mindiff = INT_MAX;
	oldlint1 = oldlint2 = oldlint3 = 0;
	while (lfp < ptr + size3)
	{
		/* find nearest integer */
		if (*lfp >= 0.0)
			lf = (float)*lfp * float_prec + 0.5;
		else
			lf = (float)*lfp * float_prec - 0.5;
		if (fabs(lf) > INT_MAX - 2)
		{
			/* scaling would cause overflow */
			fprintf(stderr, "Internal overflow compressing coordinates.\n");
			// errval=0;
		}
		lint1 = lf;
		if (lint1 < minint[0])
			minint[0] = lint1;
		if (lint1 > maxint[0])
			maxint[0] = lint1;
		*lip++ = lint1;
		lfp++;
		if (*lfp >= 0.0)
			lf = (float)*lfp * float_prec + 0.5;
		else
			lf = (float)*lfp * float_prec - 0.5;
		if (fabs(lf) > INT_MAX - 2)
		{
			/* scaling would cause overflow */
			fprintf(stderr, "Internal overflow compressing coordinates.\n");
			// errval=0;
		}
		lint2 = lf;
		if (lint2 < minint[1])
			minint[1] = lint2;
		if (lint2 > maxint[1])
			maxint[1] = lint2;
		*lip++ = lint2;
		lfp++;
		if (*lfp >= 0.0)
			lf = (float)*lfp * float_prec + 0.5;
		else
			lf = (float)*lfp * float_prec - 0.5;
		if (fabs(lf) > INT_MAX - 2)
		{
			// errval=0;
		}
		lint3 = lf;
		if (lint3 < minint[2])
			minint[2] = lint3;
		if (lint3 > maxint[2])
			maxint[2] = lint3;
		*lip++ = lint3;
		lfp++;
		diff = abs(oldlint1 - lint1) + abs(oldlint2 - lint2) + abs(oldlint3 - lint3);
		if (diff < mindiff && lfp > ptr + 3)
			mindiff = diff;
		oldlint1 = lint1;
		oldlint2 = lint2;
		oldlint3 = lint3;
	}
	xdrfile_write_int(minint, 3, xfp);
	xdrfile_write_int(maxint, 3, xfp);

	if ((float)maxint[0] - (float)minint[0] >= INT_MAX - 2 ||
		(float)maxint[1] - (float)minint[1] >= INT_MAX - 2 ||
		(float)maxint[2] - (float)minint[2] >= INT_MAX - 2)
	{
		/* turning value in unsigned by subtracting minint
		 * would cause overflow
		 */
		fprintf(stderr, "Internal overflow compressing coordinates.\n");
		// errval=0;
	}
	sizeint[0] = maxint[0] - minint[0] + 1;
	sizeint[1] = maxint[1] - minint[1] + 1;
	sizeint[2] = maxint[2] - minint[2] + 1;

	/* check if one of the sizes is to big to be multiplied */
	if ((sizeint[0] | sizeint[1] | sizeint[2]) > 0xffffff)
	{
		bitsizeint[0] = sizeofint(sizeint[0]);
		bitsizeint[1] = sizeofint(sizeint[1]);
		bitsizeint[2] = sizeofint(sizeint[2]);
		bitsize = 0; /* flag the use of large sizes */
	}
	else
	{
		bitsize = sizeofints(3, sizeint);
	}
	lip = buf1;
	luip = (unsigned int *)buf1;
	smallidx = FIRSTIDX;
	while (smallidx < LASTIDX && magicints[smallidx] < mindiff)
	{
		smallidx++;
	}
	xdrfile_write_int(&smallidx, 1, xfp);
	tmp = smallidx + 8;
	maxidx = (LASTIDX < tmp) ? LASTIDX : tmp;
	minidx = maxidx - 8; /* often this equal smallidx */
	tmp = smallidx - 1;
	tmp = (FIRSTIDX > tmp) ? FIRSTIDX : tmp;
	smaller = magicints[tmp] / 2;
	smallnum = magicints[smallidx] / 2;
	sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
	larger = magicints[maxidx] / 2;
	i = 0;
	while (i < size)
	{
		is_small = 0;
		thiscoord = (int *)(luip) + i * 3;
		if (smallidx < maxidx && i >= 1 &&
			abs(thiscoord[0] - prevcoord[0]) < larger &&
			abs(thiscoord[1] - prevcoord[1]) < larger &&
			abs(thiscoord[2] - prevcoord[2]) < larger)
		{
			is_smaller = 1;
		}
		else if (smallidx > minidx)
		{
			is_smaller = -1;
		}
		else
		{
			is_smaller = 0;
		}
		if (i + 1 < size)
		{
			if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
				abs(thiscoord[1] - thiscoord[4]) < smallnum &&
				abs(thiscoord[2] - thiscoord[5]) < smallnum)
			{
				/* interchange first with second atom for better
				 * compression of water molecules
				 */
				tmp = thiscoord[0];
				thiscoord[0] = thiscoord[3];
				thiscoord[3] = tmp;
				tmp = thiscoord[1];
				thiscoord[1] = thiscoord[4];
				thiscoord[4] = tmp;
				tmp = thiscoord[2];
				thiscoord[2] = thiscoord[5];
				thiscoord[5] = tmp;
				is_small = 1;
			}
		}
		tmpcoord[0] = thiscoord[0] - minint[0];
		tmpcoord[1] = thiscoord[1] - minint[1];
		tmpcoord[2] = thiscoord[2] - minint[2];
		if (bitsize == 0)
		{
			encodebits(buf2, bitsizeint[0], tmpcoord[0]);
			encodebits(buf2, bitsizeint[1], tmpcoord[1]);
			encodebits(buf2, bitsizeint[2], tmpcoord[2]);
		}
		else
		{
			encodeints(buf2, 3, bitsize, sizeint, tmpcoord);
		}
		prevcoord[0] = thiscoord[0];
		prevcoord[1] = thiscoord[1];
		prevcoord[2] = thiscoord[2];
		thiscoord = thiscoord + 3;
		i++;

		run = 0;
		if (is_small == 0 && is_smaller == -1)
			is_smaller = 0;
		while (is_small && run < 8 * 3)
		{
			tmpsum = 0;
			for (j = 0; j < 3; j++)
			{
				tmp = thiscoord[j] - prevcoord[j];
				tmpsum += tmp * tmp;
			}
			if (is_smaller == -1 && tmpsum >= smaller * smaller)
			{
				is_smaller = 0;
			}

			tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
			tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
			tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;

			prevcoord[0] = thiscoord[0];
			prevcoord[1] = thiscoord[1];
			prevcoord[2] = thiscoord[2];

			i++;
			thiscoord = thiscoord + 3;
			is_small = 0;
			if (i < size &&
				abs(thiscoord[0] - prevcoord[0]) < smallnum &&
				abs(thiscoord[1] - prevcoord[1]) < smallnum &&
				abs(thiscoord[2] - prevcoord[2]) < smallnum)
			{
				is_small = 1;
			}
		}
		if (run != prevrun || is_smaller != 0)
		{
			prevrun = run;
			encodebits(buf2, 1, 1); /* flag the change in run-length */
			encodebits(buf2, 5, run + is_smaller + 1);
		}
		else
		{
			encodebits(buf2, 1, 0); /* flag the fact that runlength did not change */
		}
		for (k = 0; k < run; k += 3)
		{
			encodeints(buf2, 3, smallidx, sizesmall, &tmpcoord[k]);
		}
		if (is_smaller != 0)
		{
			smallidx += is_smaller;
			if (is_smaller < 0)
			{
				smallnum = smaller;
				smaller = magicints[smallidx - 1] / 2;
			}
			else
			{
				smaller = smallnum;
				smallnum = magicints[smallidx] / 2;
			}
			sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
		}
	}
	if (buf2[1] != 0)
		buf2[0]++;
	xdrfile_write_int(buf2, 1, xfp); /* buf2[0] holds the length in bytes */
	tmp = xdrfile_write_opaque((char *)&(buf2[3]), (unsigned int)buf2[0], xfp);
	if (tmp == (unsigned int)buf2[0])
		return size;
	else
		return -1;
}

int xdr_jump(XDRFILE *xdp, long offset)
{
	if (fseek(xdp->fp, offset, SEEK_CUR) < 0)
		return 1;

	return 0;
}

/* Dont try do document Fortran interface, since
 * Doxygen barfs at the F77_FUNC macro
 */
#ifndef DOXYGEN

/*************************************************************
 * Fortran77 interface for reading/writing portable data     *
 * The routine are not threadsafe when called from Fortran   *
 * (as they are when called from C) unless you compile with  *
 * this file with posix thread support.                      *
 * Note that these are not multithread-safe.                 *
 *************************************************************/
#define MAX_FORTRAN_XDR 1024
static XDRFILE *f77xdr[MAX_FORTRAN_XDR]; /* array of file handles */
static int f77init = 1;					 /* zero array first time */

/* internal to this file: C<-->Fortran string conversion */
static int ftocstr(char *dest, int dest_len, char *src, int src_len);
static int ctofstr(char *dest, int dest_len, char *src);

void F77_FUNC(xdropen, XDROPEN)(int *fid, char *filename, char *mode,
								int fn_len, int mode_len)
{
	char cfilename[512];
	char cmode[5];
	int i;

	/* zero array at first invocation */
	if (f77init)
	{
		for (i = 0; i < MAX_FORTRAN_XDR; i++)
			f77xdr[i] = NULL;
		f77init = 0;
	}
	i = 0;

	/* nf77xdr is always smaller or equal to MAX_FORTRAN_XDR */
	while (i < MAX_FORTRAN_XDR && f77xdr[i] != NULL)
		i++;
	if (i == MAX_FORTRAN_XDR)
	{
		*fid = -1;
	}
	else if (ftocstr(cfilename, sizeof(cfilename), filename, fn_len))
	{
		*fid = -1;
	}
	else if (ftocstr(cmode, sizeof(cmode), mode, mode_len))
	{
		*fid = -1;
	}
	else
	{
		f77xdr[i] = xdrfile_open(cfilename, cmode);
		/* return the index in the array as a fortran file handle */
		*fid = i;
	}
}

void F77_FUNC(xdrclose, XDRCLOSE)(int *fid)
{
	/* first close it */
	xdrfile_close(f77xdr[*fid]);
	/* the remove it from file handle list */
	f77xdr[*fid] = NULL;
}

void F77_FUNC(xdrrint, XDRRINT)(int *fid, int *data, int *ndata, int *ret)
{
	*ret = xdrfile_read_int(data, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrwint, XDRWINT)(int *fid, int *data, int *ndata, int *ret)
{
	*ret = xdrfile_write_int(data, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrruint, XDRRUINT)(int *fid, unsigned int *data, int *ndata, int *ret)
{
	*ret = xdrfile_read_uint(data, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrwuint, XDRWUINT)(int *fid, unsigned int *data, int *ndata, int *ret)
{
	*ret = xdrfile_write_uint(data, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrrchar, XDRRCHAR)(int *fid, char *ip, int *ndata, int *ret)
{
	*ret = xdrfile_read_char(ip, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrwchar, XDRWCHAR)(int *fid, char *ip, int *ndata, int *ret)
{
	*ret = xdrfile_write_char(ip, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrruchar, XDRRUCHAR)(int *fid, unsigned char *ip, int *ndata, int *ret)
{
	*ret = xdrfile_read_uchar(ip, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrwuchar, XDRWUCHAR)(int *fid, unsigned char *ip, int *ndata, int *ret)
{
	*ret = xdrfile_write_uchar(ip, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrrshort, XDRRSHORT)(int *fid, short *ip, int *ndata, int *ret)
{
	*ret = xdrfile_read_short(ip, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrwshort, XDRWSHORT)(int *fid, short *ip, int *ndata, int *ret)
{
	*ret = xdrfile_write_short(ip, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrrushort, XDRRUSHORT)(int *fid, unsigned short *ip, int *ndata, int *ret)
{
	*ret = xdrfile_read_ushort(ip, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrwushort, XDRWUSHORT)(int *fid, unsigned short *ip, int *ndata, int *ret)
{
	*ret = xdrfile_write_ushort(ip, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrrsingle, XDRRSINGLE)(int *fid, float *data, int *ndata, int *ret)
{
	*ret = xdrfile_read_float(data, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrwsingle, XDRWSINGLE)(int *fid, float *data, int *ndata, int *ret)
{
	*ret = xdrfile_write_float(data, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrrdouble, XDRRDOUBLE)(int *fid, double *data, int *ndata, int *ret)
{
	*ret = xdrfile_read_double(data, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrwdouble, XDRWDOUBLE)(int *fid, double *data, int *ndata, int *ret)
{
	*ret = xdrfile_write_double(data, *ndata, f77xdr[*fid]);
}

static int ftocstr(char *dest, int destlen, char *src, int srclen)
{
	char *p;

	p = src + srclen;
	while (--p >= src && *p == ' ')
		;
	srclen = p - src + 1;
	destlen--;
	dest[0] = 0;
	if (srclen > destlen)
		return 1;
	while (srclen--)
		(*dest++ = *src++);
	*dest = '\0';
	return 0;
}

static int ctofstr(char *dest, int destlen, char *src)
{
	while (destlen && *src)
	{
		*dest++ = *src++;
		destlen--;
	}
	while (destlen--)
		*dest++ = ' ';
	return 0;
}

void F77_FUNC(xdrrstring, XDRRSTRING)(int *fid, char *str, int *ret, int len)
{
	char *cstr;

	if ((cstr = (char *)malloc((len + 1) * sizeof(char))) == NULL)
	{
		*ret = 0;
		return;
	}
	if (ftocstr(cstr, len + 1, str, len))
	{
		*ret = 0;
		free(cstr);
		return;
	}

	*ret = xdrfile_read_string(cstr, len + 1, f77xdr[*fid]);
	ctofstr(str, len, cstr);
	free(cstr);
}

void F77_FUNC(xdrwstring, XDRWSTRING)(int *fid, char *str, int *ret, int len)
{
	char *cstr;

	if ((cstr = (char *)malloc((len + 1) * sizeof(char))) == NULL)
	{
		*ret = 0;
		return;
	}
	if (ftocstr(cstr, len + 1, str, len))
	{
		*ret = 0;
		free(cstr);
		return;
	}

	*ret = xdrfile_write_string(cstr, f77xdr[*fid]);
	ctofstr(str, len, cstr);
	free(cstr);
}

void F77_FUNC(xdrropaque, XDRROPAQUE)(int *fid, char *data, int *ndata, int *ret)
{
	*ret = xdrfile_read_opaque(data, *ndata, f77xdr[*fid]);
}

void F77_FUNC(xdrwopaque, XDRWOPAQUE)(int *fid, char *data, int *ndata, int *ret)
{
	*ret = xdrfile_write_opaque(data, *ndata, f77xdr[*fid]);
}

/* Write single-precision compressed 3d coordinates */
void F77_FUNC(xdrccs, XDRCCS)(int *fid, float *data, int *ncoord,
							  float *precision, int *ret)
{
	*ret = xdrfile_compress_coord_float(data, *ncoord, *precision, f77xdr[*fid]);
}

/* Read single-precision compressed 3d coordinates */
void F77_FUNC(xdrdcs, XDRDCS)(int *fid, float *data, int *ncoord,
							  float *precision, int *ret)
{
	*ret = xdrfile_decompress_coord_float(data, ncoord, precision, f77xdr[*fid]);
}

/* Write compressed 3d coordinates from double precision data */
void F77_FUNC(xdrccd, XDRCCD)(int *fid, double *data, int *ncoord,
							  double *precision, int *ret)
{
	*ret = xdrfile_compress_coord_double(data, *ncoord, *precision, f77xdr[*fid]);
}

/* Read compressed 3d coordinates into double precision data */
void F77_FUNC(xddcd, XDRDCD)(int *fid, double *data, int *ncoord,
							 double *precision, int *ret)
{
	*ret = xdrfile_decompress_coord_double(data, ncoord, precision, f77xdr[*fid]);
}

#endif /* DOXYGEN */

/*************************************************************
 * End of higher-level routines - dont change things below!  *
 *************************************************************/

/*************************************************************
 * The rest of this file contains our own implementation     *
 * of the XDR calls in case you are compiling without them.  *
 * You do NOT want to change things here since it would make *
 * things incompatible with the standard RPC/XDR routines.   *
 *************************************************************/
#ifndef HAVE_RPC_XDR_H

/*
 * What follows is a modified version of the Sun XDR code. For reference
 * we include their copyright and license:
 *
 * Sun RPC is a product of Sun Microsystems, Inc. and is provided for
 * unrestricted use provided that this legend is included on all tape
 * media and as a part of the software program in whole or part.  Users
 * may copy or modify Sun RPC without charge, but are not authorized
 * to license or distribute it to anyone else except as part of a product or
 * program developed by the user.
 *
 * SUN RPC IS PROVIDED AS IS WITH NO WARRANTIES OF ANY KIND INCLUDING THE
 * WARRANTIES OF DESIGN, MERCHANTIBILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE PRACTICE.
 *
 * Sun RPC is provided with no support and without any obligation on the
 * part of Sun Microsystems, Inc. to assist in its use, correction,
 * modification or enhancement.
 *
 * SUN MICROSYSTEMS, INC. SHALL HAVE NO LIABILITY WITH RESPECT TO THE
 * INFRINGEMENT OF COPYRIGHTS, TRADE SECRETS OR ANY PATENTS BY SUN RPC
 * OR ANY PART THEREOF.
 *
 * In no event will Sun Microsystems, Inc. be liable for any lost revenue
 * or profits or other special, indirect and consequential damages, even if
 * Sun has been advised of the possibility of such damages.
 *
 * Sun Microsystems, Inc.
 * 2550 Garcia Avenue
 * Mountain View, California  94043
 */

/* INT_MAX is defined in limits.h according to ANSI C */
#if (INT_MAX > 2147483647)
#error Error: Cannot use builtin XDR support when size of int
#error is larger than 4 bytes. Use your system XDR libraries
#error instead, or modify the source code in xdrfile.c
#endif /* Check for 4 byte int type */

typedef int (*xdrproc_t)(XDR *, void *, ...);

#define xdr_getlong(xdrs, longp) \
	(*(xdrs)->x_ops->x_getlong)(xdrs, longp)
#define xdr_putlong(xdrs, longp) \
	(*(xdrs)->x_ops->x_putlong)(xdrs, longp)
#define xdr_getbytes(xdrs, addr, len) \
	(*(xdrs)->x_ops->x_getbytes)(xdrs, addr, len)
#define xdr_putbytes(xdrs, addr, len) \
	(*(xdrs)->x_ops->x_putbytes)(xdrs, addr, len)

#define BYTES_PER_XDR_UNIT 4
static char xdr_zero[BYTES_PER_XDR_UNIT] = {0, 0, 0, 0};

static int32_t
xdr_swapbytes(int32_t x)
{
	int32_t y, i;
	char *px = (char *)&x;
	char *py = (char *)&y;

	for (i = 0; i < 4; i++)
		py[i] = px[3 - i];

	return y;
}

static int32_t
xdr_htonl(int32_t x)
{
	int s = 0x1234;
	if (*((char *)&s) == (char)0x34)
	{
		/* smallendian,swap bytes */
		return xdr_swapbytes(x);
	}
	else
	{
		/* bigendian, do nothing */
		return x;
	}
}

static int32_t
xdr_ntohl(int x)
{
	int s = 0x1234;
	if (*((char *)&s) == (char)0x34)
	{
		/* smallendian, swap bytes */
		return xdr_swapbytes(x);
	}
	else
	{
		/* bigendian, do nothing */
		return x;
	}
}

static int
xdr_int(XDR *xdrs, int *ip)
{
	int32_t i32;

	switch (xdrs->x_op)
	{
	case XDR_ENCODE:
		i32 = (int32_t)*ip;
		return xdr_putlong(xdrs, &i32);

	case XDR_DECODE:
		if (!xdr_getlong(xdrs, &i32))
		{
			return 0;
		}
		*ip = (int)i32;
	case XDR_FREE:
		return 1;
	}
	return 0;
}

static int
xdr_u_int(XDR *xdrs, unsigned int *up)
{
	uint32_t ui32;

	switch (xdrs->x_op)
	{
	case XDR_ENCODE:
		ui32 = (uint32_t)*up;
		return xdr_putlong(xdrs, (int32_t *)&ui32);

	case XDR_DECODE:
		if (!xdr_getlong(xdrs, (int32_t *)&ui32))
		{
			return 0;
		}
		*up = (uint32_t)ui32;
	case XDR_FREE:
		return 1;
	}
	return 0;
}

static int
xdr_short(XDR *xdrs, short *sp)
{
	int32_t i32;

	switch (xdrs->x_op)
	{
	case XDR_ENCODE:
		i32 = (int32_t)*sp;
		return xdr_putlong(xdrs, &i32);

	case XDR_DECODE:
		if (!xdr_getlong(xdrs, &i32))
		{
			return 0;
		}
		*sp = (short)i32;
		return 1;

	case XDR_FREE:
		return 1;
	}
	return 0;
}

static int
xdr_u_short(XDR *xdrs, unsigned short *sp)
{
	uint32_t ui32;

	switch (xdrs->x_op)
	{
	case XDR_ENCODE:
		ui32 = (uint32_t)*sp;
		return xdr_putlong(xdrs, (int32_t *)&ui32);

	case XDR_DECODE:
		if (!xdr_getlong(xdrs, (int32_t *)&ui32))
		{
			return 0;
		}
		*sp = (unsigned short)ui32;
		return 1;

	case XDR_FREE:
		return 1;
	}
	return 0;
}

static int
xdr_char(XDR *xdrs, char *cp)
{
	int i;

	i = (*cp);
	if (!xdr_int(xdrs, &i))
	{
		return 0;
	}
	*cp = i;
	return 1;
}

static int
xdr_u_char(XDR *xdrs, unsigned char *cp)
{
	unsigned int u;

	u = (*cp);
	if (!xdr_u_int(xdrs, &u))
	{
		return 0;
	}
	*cp = u;
	return 1;
}

/*
 * XDR opaque data
 * Allows the specification of a fixed size sequence of opaque bytes.
 * cp points to the opaque object and cnt gives the byte length.
 */
static int
xdr_opaque(XDR *xdrs, char *cp, unsigned int cnt)
{
	unsigned int rndup;
	static char crud[BYTES_PER_XDR_UNIT];

	/*
	 * if no data we are done
	 */
	if (cnt == 0)
		return 1;

	/*
	 * round byte count to full xdr units
	 */
	rndup = cnt % BYTES_PER_XDR_UNIT;
	if (rndup > 0)
		rndup = BYTES_PER_XDR_UNIT - rndup;

	switch (xdrs->x_op)
	{
	case XDR_DECODE:
		if (!xdr_getbytes(xdrs, cp, cnt))
		{
			return 0;
		}
		if (rndup == 0)
			return 1;
		return xdr_getbytes(xdrs, (char *)crud, rndup);

	case XDR_ENCODE:
		if (!xdr_putbytes(xdrs, cp, cnt))
		{
			return 0;
		}
		if (rndup == 0)
			return 1;
		return xdr_putbytes(xdrs, xdr_zero, rndup);

	case XDR_FREE:
		return 1;
	}
#undef BYTES_PER_XDR_UNIT
	return 0;
}

/*
 * XDR null terminated ASCII strings
 */
static int
xdr_string(XDR *xdrs, char **cpp, unsigned int maxsize)
{
	char *sp = *cpp; /* sp is the actual string pointer */
	unsigned int size = 0;
	unsigned int nodesize = 0;

	/*
	 * first deal with the length since xdr strings are counted-strings
	 */
	switch (xdrs->x_op)
	{
	case XDR_FREE:
		if (sp == NULL)
		{
			return 1; /* already free */
		}
		/* fall through... */
	case XDR_ENCODE:
		if (sp == NULL)
			return 0;
		size = strlen(sp);
		break;
	case XDR_DECODE:
		break;
	}
	if (!xdr_u_int(xdrs, &size))
	{
		return 0;
	}
	if (size > maxsize)
	{
		return 0;
	}
	nodesize = size + 1;

	/*
	 * now deal with the actual bytes
	 */
	switch (xdrs->x_op)
	{
	case XDR_DECODE:
		if (nodesize == 0)
		{
			return 1;
		}
		if (sp == NULL)
			*cpp = sp = (char *)malloc(nodesize);
		if (sp == NULL)
		{
			(void)fputs("xdr_string: out of memory\n", stderr);
			return 0;
		}
		sp[size] = 0;
		/* fall into ... */

	case XDR_ENCODE:
		return xdr_opaque(xdrs, sp, size);

	case XDR_FREE:
		free(sp);
		*cpp = NULL;
		return 1;
	}
	return 0;
}

/* Floating-point stuff */

static int
xdr_float(XDR *xdrs, float *fp)
{
	switch (xdrs->x_op)
	{

	case XDR_ENCODE:
		if (sizeof(float) == sizeof(int32_t))
			return (xdr_putlong(xdrs, (int32_t *)fp));
		else if (sizeof(float) == sizeof(int))
		{
			int32_t tmp = *(int *)fp;
			return (xdr_putlong(xdrs, &tmp));
		}
		break;

	case XDR_DECODE:
		if (sizeof(float) == sizeof(int32_t))
			return (xdr_getlong(xdrs, (int32_t *)fp));
		else if (sizeof(float) == sizeof(int))
		{
			int32_t tmp;
			if (xdr_getlong(xdrs, &tmp))
			{
				*(int *)fp = tmp;
				return (1);
			}
		}
		break;

	case XDR_FREE:
		return (1);
	}
	return (0);
}

static int
xdr_double(XDR *xdrs, double *dp)
{
	/* Gromacs detects floating-point stuff at compile time, which is faster */
#ifdef GROMACS
#ifndef FLOAT_FORMAT_IEEE754
#error non-IEEE floating point system, or you defined GROMACS yourself...
#endif
	int LSW;
#ifdef IEEE754_BIG_ENDIAN_WORD_ORDER
	int LSW = 1;
#else
	int LSW = 0;
#endif /* Big endian word order */
#else
	/* Outside Gromacs we rely on dynamic detection of FP order. */
	int LSW; /* Least significant fp word */

	double x = 0.987654321; /* Just a number */
	unsigned char ix = *((char *)&x);

	/* Possible representations in IEEE double precision:
	 * (S=small endian, B=big endian)
	 *
	 * Byte order, Word order, Hex
	 *     S           S       b8 56 0e 3c dd 9a ef 3f
	 *     B           S       3c 0e 56 b8 3f ef 9a dd
	 *     S           B       dd 9a ef 3f b8 56 0e 3c
	 *     B           B       3f ef 9a dd 3c 0e 56 b8
	 */
	if (ix == 0xdd || ix == 0x3f)
		LSW = 1; /* Big endian word order */
	else if (ix == 0xb8 || ix == 0x3c)
		LSW = 0; /* Small endian word order */
	else
	{ /* Catch strange errors */
		fprintf(stderr, "Cannot detect floating-point word order.\n"
						"Do you have a non-IEEE system?\n"
						"Use system XDR libraries or fix xdr_double().\n");
		abort();
	}
#endif /* end of dynamic detection of fp word order */

	switch (xdrs->x_op)
	{

	case XDR_ENCODE:
		if (2 * sizeof(int32_t) == sizeof(double))
		{
			int32_t *lp = (int32_t *)dp;
			return (xdr_putlong(xdrs, lp + !LSW) &&
					xdr_putlong(xdrs, lp + LSW));
		}
		else if (2 * sizeof(int) == sizeof(double))
		{
			int *ip = (int *)dp;
			int32_t tmp[2];
			tmp[0] = ip[!LSW];
			tmp[1] = ip[LSW];
			return (xdr_putlong(xdrs, tmp) &&
					xdr_putlong(xdrs, tmp + 1));
		}
		break;

	case XDR_DECODE:
		if (2 * sizeof(int32_t) == sizeof(double))
		{
			int32_t *lp = (int32_t *)dp;
			return (xdr_getlong(xdrs, lp + !LSW) &&
					xdr_getlong(xdrs, lp + LSW));
		}
		else if (2 * sizeof(int) == sizeof(double))
		{
			int *ip = (int *)dp;
			int32_t tmp[2];
			if (xdr_getlong(xdrs, tmp + !LSW) &&
				xdr_getlong(xdrs, tmp + LSW))
			{
				ip[0] = tmp[0];
				ip[1] = tmp[1];
				return (1);
			}
		}
		break;

	case XDR_FREE:
		return (1);
	}
	return (0);
}

static int xdrstdio_getlong(XDR *, int32_t *);
static int xdrstdio_putlong(XDR *, int32_t *);
static int xdrstdio_getbytes(XDR *, char *, unsigned int);
static int xdrstdio_putbytes(XDR *, char *, unsigned int);
static unsigned int xdrstdio_getpos(XDR *);
static int xdrstdio_setpos(XDR *, unsigned int);
static void xdrstdio_destroy(XDR *);

/*
 * Ops vector for stdio type XDR
 */
static const struct xdr_ops xdrstdio_ops =
	{
		xdrstdio_getlong,  /* deserialize a long int */
		xdrstdio_putlong,  /* serialize a long int */
		xdrstdio_getbytes, /* deserialize counted bytes */
		xdrstdio_putbytes, /* serialize counted bytes */
		xdrstdio_getpos,   /* get offset in the stream */
		xdrstdio_setpos,   /* set offset in the stream */
		xdrstdio_destroy,  /* destroy stream */
};

/*
 * Initialize a stdio xdr stream.
 * Sets the xdr stream handle xdrs for use on the stream file.
 * Operation flag is set to op.
 */
static void
xdrstdio_create(XDR *xdrs, FILE *file, enum xdr_op op)
{
	xdrs->x_op = op;

	xdrs->x_ops = (struct xdr_ops *)&xdrstdio_ops;
	xdrs->x_private = (char *)file;
}

/*
 * Destroy a stdio xdr stream.
 * Cleans up the xdr stream handle xdrs previously set up by xdrstdio_create.
 */
static void
xdrstdio_destroy(XDR *xdrs)
{
	(void)fflush((FILE *)xdrs->x_private);
	/* xx should we close the file ?? */
}

static int
xdrstdio_getlong(XDR *xdrs, int32_t *lp)
{
	int32_t mycopy;

	if (fread((char *)&mycopy, 4, 1, (FILE *)xdrs->x_private) != 1)
		return 0;
	*lp = (int32_t)xdr_ntohl(mycopy);
	return 1;
}

static int
xdrstdio_putlong(XDR *xdrs, int32_t *lp)
{
	int32_t mycopy = xdr_htonl(*lp);
	lp = &mycopy;
	if (fwrite((char *)lp, 4, 1, (FILE *)xdrs->x_private) != 1)
		return 0;
	return 1;
}

static int
xdrstdio_getbytes(XDR *xdrs, char *addr, unsigned int len)
{
	if ((len != 0) && (fread(addr, (int)len, 1,
							 (FILE *)xdrs->x_private) != 1))
		return 0;
	return 1;
}

static int
xdrstdio_putbytes(XDR *xdrs, char *addr, unsigned int len)
{
	if ((len != 0) && (fwrite(addr, (int)len, 1,
							  (FILE *)xdrs->x_private) != 1))
		return 0;
	return 1;
}

/* 32 bit fileseek operations */
static unsigned int
xdrstdio_getpos(XDR *xdrs)
{
	return (unsigned int)ftell((FILE *)xdrs->x_private);
}

static int
xdrstdio_setpos(XDR *xdrs, unsigned int pos)
{
	return fseek((FILE *)xdrs->x_private, pos, 0) < 0 ? 0 : 1;
}

#endif /* HAVE_RPC_XDR_H not defined */
