/*
 * $Id$
 * 
 *                This source code is part of
 * 
 *                 G   R   O   M   A   C   S
 * 
 *          GROningen MAchine for Chemical Simulations
 * 
 *                        VERSION 3.2.0
 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 * Copyright (c) 2001-2004, The GROMACS development team,
 * check out http://www.gromacs.org for more information.

 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 * 
 * If you want to redistribute modifications, please consider that
 * scientific software is very special. Version control is crucial -
 * bugs must be traceable. We will be happy to consider code for
 * inclusion in the official distribution, but derived work must not
 * be called official GROMACS. Details are found in the README & COPYING
 * files - if they are missing, get the official version at www.gromacs.org.
 * 
 * To help us fund GROMACS development, we humbly ask that you cite
 * the papers on the package - you can find them in the top README file.
 * 
 * For more info, check our website at http://www.gromacs.org
 * 
 * And Hey:
 * GROningen Mixture of Alchemy and Childrens' Stories
 */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "statutil.h"
#include "xdrf.h"




#ifdef HAVE_FSEEKO
#define gmx_fseek(A,B,C) fseeko(A,B,C)
#define gmx_ftell(A) ftello(A)
#define gmx_off_t off_t
#ifndef _LARGEFILE_SOURCE
#warning Declaring fseeko and ftello explicitly so as to bypass an autoconf 2.61 bug
#warning If you have an incompatible declaration error please comment out the fseeko and ftello declarations just after this warning
int fseeko(FILE *stream, off_t offset, int whence);
off_t ftello(FILE *stream);
#endif
#else
#define gmx_fseek(A,B,C) fseek(A,B,C)
#define gmx_ftell(A) ftell(A)
#define gmx_off_t int
#endif


#define MAXID 256
static FILE *xdrfiles[MAXID];
static XDR *xdridptr[MAXID];
static char xdrmodes[MAXID];
static unsigned int cnt;


#ifdef GMX_FORTRAN

typedef void (* F77_FUNC(xdrfproc,XDRFPROC))(int *, void *, int *);

int ftocstr(char *ds, int dl, char *ss, int sl)
    /* dst, src ptrs */
    /* dst max len */
    /* src len */
{
    char *p;

    p = ss + sl;
    while ( --p >= ss && *p == ' ' );
    sl = p - ss + 1;
    dl--;
    ds[0] = 0;
    if (sl > dl)
      return 1;
    while (sl--)
      (*ds++ = *ss++);
    *ds = '\0';
    return 0;
}


int ctofstr(char *ds, int dl, char *ss)
     /* dest space */
     /* max dest length */
     /* src string (0-term) */
{
    while (dl && *ss) {
    *ds++ = *ss++;
    dl--;
    }
    while (dl--)
    *ds++ = ' ';
    return 0;
}

void
F77_FUNC(xdrfbool,XDRFBOOL)(int *xdrid, int *pb, int *ret) 
{
    *ret = xdr_bool(xdridptr[*xdrid], pb);
    cnt += sizeof(int);
}

void
F77_FUNC(xdrfchar,XDRFCHAR)(int *xdrid, char *cp, int *ret)
{
    *ret = xdr_char(xdridptr[*xdrid], cp);
    cnt += sizeof(char);
}

void
F77_FUNC(xdrfdouble,XDRFDOUBLE)(int *xdrid, double *dp, int *ret)
{
    *ret = xdr_double(xdridptr[*xdrid], dp);
    cnt += sizeof(double);
}

void
F77_FUNC(xdrffloat,XDRFFLOAT)(int *xdrid, float *fp, int *ret)
{
    *ret = xdr_float(xdridptr[*xdrid], fp);
    cnt += sizeof(float);
}

void
F77_FUNC(xdrfint,XDRFINT)(int *xdrid, int *ip, int *ret)
{
    *ret = xdr_int(xdridptr[*xdrid], ip);
    cnt += sizeof(int);
}

F77_FUNC(xdrfshort,XDRFSHORT)(int *xdrid, short *sp, int *ret)
{
    *ret = xdr_short(xdridptr[*xdrid], sp);
      cnt += sizeof(sp);
}

void
F77_FUNC(xdrfuchar,XDRFUCHAR)(int *xdrid, unsigned char *ucp, int *ret)
{
    *ret = xdr_u_char(xdridptr[*xdrid], (u_char *)ucp);
    cnt += sizeof(char);
}


void
F77_FUNC(xdrfushort,XDRFUSHORT)(int *xdrid, unsigned short *usp, int *ret)
{
    *ret = xdr_u_short(xdridptr[*xdrid], (unsigned short *)usp);
    cnt += sizeof(unsigned short);
}

void 
F77_FUNC(xdrf3dfcoord,XDRF3DFCOORD)(int *xdrid, float *fp, int *size, float *precision, int *ret)
{
    *ret = xdr3dfcoord(xdridptr[*xdrid], fp, size, precision);
}

void
F77_FUNC(xdrfstring,XDRFSTRING)(int *xdrid, char * sp_ptr,
                int *maxsize, int *ret, int sp_len)
{
    char *tsp;

    tsp = (char*) malloc((size_t)(((sp_len) + 1) * sizeof(char)));
    if (tsp == NULL) {
        *ret = -1;
        return;
    }
    if (ftocstr(tsp, *maxsize+1, sp_ptr, sp_len)) {
        *ret = -1;
        free(tsp);
        return;
    }
        *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (unsigned int) *maxsize);
    ctofstr( sp_ptr, sp_len , tsp);
    cnt += *maxsize;
    free(tsp);
}

void
F77_FUNC(xdrfwrapstring,XDRFWRAPSTRING)(int *xdrid, char *sp_ptr,
                    int *ret, int sp_len)
{
    char *tsp;
    int maxsize;
    maxsize = (sp_len) + 1;
    tsp = (char*) malloc((size_t)(maxsize * sizeof(char)));
    if (tsp == NULL) {
        *ret = -1;
        return;
    }
    if (ftocstr(tsp, maxsize, sp_ptr, sp_len)) {
        *ret = -1;
        free(tsp);
        return;
    }
    *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int)maxsize);
    ctofstr( sp_ptr, sp_len, tsp);
    cnt += maxsize;
    free(tsp);
}

void
F77_FUNC(xdrfopaque,XDRFOPAQUE)(int *xdrid, caddr_t *cp, int *ccnt, int *ret)
{
    *ret = xdr_opaque(xdridptr[*xdrid], (caddr_t)*cp, (u_int)*ccnt);
    cnt += *ccnt;
}

void
F77_FUNC(xdrfsetpos,XDRFSETPOS)(int *xdrid, int *pos, int *ret)
{
    *ret = xdr_setpos(xdridptr[*xdrid], (u_int) *pos);
}


void
F77_FUNC(xdrf,XDRF)(int *xdrid, int *pos)
{
    *pos = xdr_getpos(xdridptr[*xdrid]);
}

void
F77_FUNC(xdrfvector,XDRFVECTOR)(int *xdrid, char *cp, int *size, F77_FUNC(xdrfproc,XDRFPROC) elproc, int *ret) 
{
    int lcnt;
    cnt = 0;
    for (lcnt = 0; lcnt < *size; lcnt++) {
        elproc(xdrid, (cp+cnt) , ret);
    }
}


void
F77_FUNC(xdrfclose,XDRFCLOSE)(int *xdrid, int *ret)
{
    *ret = xdrclose(xdridptr[*xdrid]);
    cnt = 0;
}

void
F77_FUNC(xdrfopen,XDRFOPEN)(int *xdrid, char *fp_ptr, char *mode_ptr,
                int *ret, int fp_len, int mode_len)
{
    char fname[512];
    char fmode[3];

    if (ftocstr(fname, sizeof(fname), fp_ptr, fp_len)) {
        *ret = 0;
    }
    if (ftocstr(fmode, sizeof(fmode), mode_ptr,
            mode_len)) {
        *ret = 0;
    }

    *xdrid = xdropen(NULL, fname, fmode);
    if (*xdrid == 0)
        *ret = 0;
    else 
        *ret = 1;    
}
#endif /* GMX_FORTRAN */

/*___________________________________________________________________________
 |
 | what follows are the C routines for opening, closing xdr streams
 | and the routine to read/write compressed coordinates together
 | with some routines to assist in this task (those are marked
 | static and cannot be called from user programs)
*/
#define MAXABS INT_MAX-2

#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x):(y))
#endif
#ifndef MAX
#define MAX(x,y) ((x) > (y) ? (x):(y))
#endif
#ifndef SQR
#define SQR(x) ((x)*(x))
#endif
static 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))


/*__________________________________________________________________________
 |
 | xdropen - open xdr file
 |
 | This versions differs from xdrstdio_create, because I need to know
 | the state of the file (read or write) so I can use xdr3dfcoord
 | in eigther read or write mode, and the file descriptor
 | so I can close the file (something xdr_destroy doesn't do).
 |
*/

int xdropen(XDR *xdrs, const char *filename, const char *type) {
    static int init_done = 0;
    enum xdr_op lmode;
    int xdrid;
    char newtype[5];

    if (init_done == 0) {
    for (xdrid = 1; xdrid < MAXID; xdrid++) {
        xdridptr[xdrid] = NULL;
    }
    init_done = 1;
    }
    xdrid = 1;
    while (xdrid < MAXID && xdridptr[xdrid] != NULL) {
    xdrid++;
    }
    if (xdrid == MAXID) {
    return 0;
    }
    if (*type == 'w' || *type == 'W') {
            strcpy(newtype,"wb+");
        lmode = XDR_ENCODE;
    } else if (*type == 'a' || *type == 'A') {
            strcpy(newtype,"ab+");
            lmode = XDR_ENCODE;
    } else {
            strcpy(newtype,"rb");
        lmode = XDR_DECODE;
    }
    xdrfiles[xdrid] = fopen(filename, newtype);
    if (xdrfiles[xdrid] == NULL) {
    xdrs = NULL;
    return 0;
    }
    xdrmodes[xdrid] = *type;
    /* next test isn't usefull in the case of C language
     * but is used for the Fortran interface
     * (C users are expected to pass the address of an already allocated
     * XDR staructure)
     */
    if (xdrs == NULL) {
    xdridptr[xdrid] = (XDR *) malloc((size_t)sizeof(XDR));
    xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode);
    } else {
    xdridptr[xdrid] = xdrs;
    xdrstdio_create(xdrs, xdrfiles[xdrid], lmode);
    }
    return xdrid;
}

/*_________________________________________________________________________
 |
 | xdrclose - close a xdr file
 |
 | This will flush the xdr buffers, and destroy the xdr stream.
 | It also closes the associated file descriptor (this is *not*
 | done by xdr_destroy).
 |
*/
 
int xdrclose(XDR *xdrs) {
    int xdrid;
    
    if (xdrs == NULL) {
    fprintf(stderr, "xdrclose: passed a NULL pointer\n");
    exit(1);
    }
    for (xdrid = 1; xdrid < MAXID; xdrid++) {
    if (xdridptr[xdrid] == xdrs) {
        
        xdr_destroy(xdrs);
        fclose(xdrfiles[xdrid]);
        xdridptr[xdrid] = NULL;
        return 1;
    }
    } 
    fprintf(stderr, "xdrclose: no such open xdr file\n");
    exit(1);
    
    /* to make some compilers happy: */
    return 0;    
}

/*____________________________________________________________________________
 |
 | sendbits - 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
 | better make sure that this number of bits is enough to hold the value
 | Also num must be positive.
 |
*/

static void sendbits(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);
    }
}

/*_________________________________________________________________________
 |
 | sizeofint - calculate bitsize of an integer
 |
 | return the number of bits needed to store an integer with given max size
 |
*/

static int sizeofint(const 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 the number of small unsigned integers and the maximum value
 | return the number of bits needed to read or write them with the
 | routines receiveints and sendints. You need this parameter when
 | calling these routines. Note that for many calls I can use
 | the variable 'smallidx' which is exactly the number of bits, and
 | So I don't need to call 'sizeofints for those calls.
*/

static int sizeofints( const 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;

}
    
/*____________________________________________________________________________
 |
 | sendints - send a small set of small integers in compressed format
 |
 | this routine is used internally by xdr3dfcoord, to send a set of
 | small integers to the buffer. 
 | 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 causes bad results.
 |
 */
 
static void sendints(int buf[], const int num_of_ints, const 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 sendints num %u doesn't "
            "match size %u\n", nums[i], sizes[i]);
        exit(1);
    }
    /* 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++) {
        sendbits(buf, 8, bytes[i]);
    }
    sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
    } else {
    for (i = 0; i < num_of_bytes-1; i++) {
        sendbits(buf, 8, bytes[i]);
    }
    sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
    }
}


/*___________________________________________________________________________
 |
 | receivebits - 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 receivebits(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; 
}

/*____________________________________________________________________________
 |
 | receiveints - decode 'small' integers from the buf array
 |
 | this routine is the inverse from sendints() 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 receiveints(int buf[], const 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++] = receivebits(buf, 8);
    num_of_bits -= 8;
    }
    if (num_of_bits > 0) {
    bytes[num_of_bytes++] = receivebits(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);
}
    
/*____________________________________________________________________________
 |
 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
 |
 | this routine reads or writes (depending on how you opened the file with
 | xdropen() ) a large number of 3d coordinates (stored in *fp).
 | The number of coordinates triplets to write is given by *size. On
 | read this number may be zero, in which case it reads as many as were written
 | or it may specify the number if triplets to read (which should match the
 | number written).
 | Compression is achieved by first converting all floating numbers to integer
 | using multiplication by *precision and rounding to the nearest integer.
 | Then the minimum and maximum value are calculated to determine the range.
 | The limited range of integers so found, is used to compress the coordinates.
 | In addition the differences between succesive coordinates is calculated.
 | If the difference happens to be 'small' then only the difference is saved,
 | compressing the data even more. The notion of 'small' is changed dynamically
 | and is enlarged or reduced whenever needed or possible.
 | Extra compression is achieved in the case of GROMOS and coordinates of
 | water molecules. GROMOS first writes out the Oxygen position, followed by
 | the two hydrogens. In order to make the differences smaller (and thereby
 | compression the data better) the order is changed into first one hydrogen
 | then the oxygen, followed by the other hydrogen. This is rather special, but
 | it shouldn't harm in the general case.
 |
 */
 
int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
    

    static int *ip = NULL;
    static int oldsize;
    static int *buf;

    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 flag, k;
    int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun;
    float *lfp, lf;
    int tmp, *thiscoord,  prevcoord[3];
    unsigned int tmpcoord[30];

    int bufsize, xdrid, lsize;
    unsigned int bitsize;
    float inv_precision;
    int errval = 1;

    bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0;
    prevcoord[0]  = prevcoord[1]  = prevcoord[2]  = 0;
    
    /* find out if xdrs is opened for reading or for writing */
    xdrid = 0;
    while (xdridptr[xdrid] != xdrs) {
    xdrid++;
    if (xdrid >= MAXID) {
        fprintf(stderr, "xdr error. no open xdr stream\n");
        exit (1);
    }
    }
    if ((xdrmodes[xdrid] == 'w') || (xdrmodes[xdrid] == 'a')) {

    /* xdrs is open for writing */

    if (xdr_int(xdrs, size) == 0)
        return 0;
    size3 = *size * 3;
    /* when the number of coordinates is small, don't try to compress; just
     * write them as floats using xdr_vector
     */
    if (*size <= 9 ) {
            return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3, 
                    (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
    }
    
    xdr_float(xdrs, precision);
    if (ip == NULL) {
        ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
        if (ip == NULL) {
        fprintf(stderr,"malloc failed\n");
        exit(1);
        }
        bufsize = size3 * 1.2;
        buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
        if (buf == NULL) {
        fprintf(stderr,"malloc failed\n");
        exit(1);
        }
        oldsize = *size;
    } else if (*size > oldsize) {
        ip = (int *)realloc(ip, (size_t)(size3 * sizeof(*ip)));
        if (ip == NULL) {
        fprintf(stderr,"malloc failed\n");
        exit(1);
        }
        bufsize = size3 * 1.2;
        buf = (int *)realloc(buf, (size_t)(bufsize * sizeof(*buf)));
        if (buf == NULL) {
        fprintf(stderr,"realloc failed\n");
        exit(1);
        }
        oldsize = *size;
    }
    /* buf[0-2] are special and do not contain actual data */
    buf[0] = buf[1] = buf[2] = 0;
    minint[0] = minint[1] = minint[2] = INT_MAX;
    maxint[0] = maxint[1] = maxint[2] = INT_MIN;
    prevrun = -1;
    lfp = fp;
    lip = ip;
    mindiff = INT_MAX;
    oldlint1 = oldlint2 = oldlint3 = 0;
    while(lfp < fp + size3 ) {
        /* find nearest integer */
        if (*lfp >= 0.0)
        lf = *lfp * *precision + 0.5;
        else
        lf = *lfp * *precision - 0.5;
        if (fabs(lf) > MAXABS) {
        /* scaling would cause overflow */
        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) > MAXABS) {
        /* scaling would cause overflow */
        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) > MAXABS) {
        /* scaling would cause overflow */
        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 > fp + 3)
        mindiff = diff;
        oldlint1 = lint1;
        oldlint2 = lint2;
        oldlint3 = lint3;
    }
    xdr_int(xdrs, &(minint[0]));
    xdr_int(xdrs, &(minint[1]));
    xdr_int(xdrs, &(minint[2]));
    
    xdr_int(xdrs, &(maxint[0]));
    xdr_int(xdrs, &(maxint[1]));
    xdr_int(xdrs, &(maxint[2]));
    
    if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
        (float)maxint[1] - (float)minint[1] >= MAXABS ||
        (float)maxint[2] - (float)minint[2] >= MAXABS) {
        /* turning value in unsigned by subtracting minint
         * would cause overflow
         */
        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 = ip;
    luip = (unsigned int *) ip;
    smallidx = FIRSTIDX;
    while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
        smallidx++;
    }
    xdr_int(xdrs, &smallidx);
    maxidx = MIN(LASTIDX, smallidx + 8) ;
    minidx = maxidx - 8; /* often this equal smallidx */
    smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 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) {
        sendbits(buf, bitsizeint[0], tmpcoord[0]);
        sendbits(buf, bitsizeint[1], tmpcoord[1]);
        sendbits(buf, bitsizeint[2], tmpcoord[2]);
        } else {
        sendints(buf, 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) {
        if (is_smaller == -1 && (
            SQR(thiscoord[0] - prevcoord[0]) +
            SQR(thiscoord[1] - prevcoord[1]) +
            SQR(thiscoord[2] - prevcoord[2]) >= 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;
        sendbits(buf, 1, 1); /* flag the change in run-length */
        sendbits(buf, 5, run+is_smaller+1);
        } else {
        sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
        }
        for (k=0; k < run; k+=3) {
        sendints(buf, 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 (buf[1] != 0) buf[0]++;
    xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */
        return errval * (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]));
    } else {
    
    /* xdrs is open for reading */
    
    if (xdr_int(xdrs, &lsize) == 0) 
        return 0;
    if (*size != 0 && lsize != *size) {
        fprintf(stderr, "wrong number of coordinates in xdr3dfcoord; "
            "%d arg vs %d in file", *size, lsize);
    }
    *size = lsize;
    size3 = *size * 3;
    if (*size <= 9) {
        *precision = -1;
            return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3, 
                    (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
    }
    xdr_float(xdrs, precision);
    if (ip == NULL) {
        ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
        if (ip == NULL) {
        fprintf(stderr,"malloc failed\n");
        exit(1);
        }
        bufsize = size3 * 1.2;
        buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
        if (buf == NULL) {
        fprintf(stderr,"malloc failed\n");
        exit(1);
        }
        oldsize = *size;
    } else if (*size > oldsize) {
        ip = (int *)realloc(ip, (size_t)(size3 * sizeof(*ip)));
        if (ip == NULL) {
        fprintf(stderr,"malloc failed\n");
        exit(1);
        }
        bufsize = size3 * 1.2;
        buf = (int *)realloc(buf, (size_t)(bufsize * sizeof(*buf)));
        if (buf == NULL) {
        fprintf(stderr,"malloc failed\n");
        exit(1);
        }
        oldsize = *size;
    }
    buf[0] = buf[1] = buf[2] = 0;
    
    xdr_int(xdrs, &(minint[0]));
    xdr_int(xdrs, &(minint[1]));
    xdr_int(xdrs, &(minint[2]));

    xdr_int(xdrs, &(maxint[0]));
    xdr_int(xdrs, &(maxint[1]));
    xdr_int(xdrs, &(maxint[2]));
        
    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 (xdr_int(xdrs, &smallidx) == 0)    
        return 0;
    maxidx = MIN(LASTIDX, smallidx + 8) ;
    minidx = maxidx - 8; /* often this equal smallidx */
    smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
    smallnum = magicints[smallidx] / 2;
    sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
    larger = magicints[maxidx];

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

    if (xdr_int(xdrs, &(buf[0])) == 0)
        return 0;
        if (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]) == 0)
        return 0;
    buf[0] = buf[1] = buf[2] = 0;
    
    lfp = fp;
    inv_precision = 1.0 / * precision;
    run = 0;
    i = 0;
    lip = ip;
    while ( i < lsize ) {
        thiscoord = (int *)(lip) + i * 3;

        if (bitsize == 0) {
        thiscoord[0] = receivebits(buf, bitsizeint[0]);
        thiscoord[1] = receivebits(buf, bitsizeint[1]);
        thiscoord[2] = receivebits(buf, bitsizeint[2]);
        } else {
        receiveints(buf, 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 = receivebits(buf, 1);
        is_smaller = 0;
        if (flag == 1) {
        run = receivebits(buf, 5);
        is_smaller = run % 3;
        run -= is_smaller;
        is_smaller--;
        }
        if (run > 0) {
        thiscoord += 3;
        for (k = 0; k < run; k+=3) {
            receiveints(buf, 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 1;
}



/******************************************************************

  XTC files have a relatively simple structure.
  They have a header of 16 bytes and the rest are
  the compressed coordinates of the files. Due to the
  compression 00 is not present in the coordinates.
  The first 4 bytes of the header are the magic number
  1995 (0x000007CB). If we find this number we are guaranteed
  to be in the header, due to the presence of so many zeros.
  The second 4 bytes are the number of atoms in the frame, and is
  assumed to be constant. The third 4 bytes are the frame number.
  The last 4 bytes are a floating point representation of the time.

********************************************************************/

/* Must match definition in xtcio.c */
#ifndef XTC_MAGIC
#define XTC_MAGIC 1995
#endif

static const int header_size = 16;


/* Check if we are at the header start.
   At the same time it will also read 1 int */
static int xtc_at_header_start(int fp, int natoms, int * timestep, float * time){
  int i_inp[3];
  float f_inp[10];
  int i;
  gmx_off_t off;


  if((off = gmx_ftell(xdrfiles[fp+1])) < 0){
    return -1;
  }
  /* read magic natoms and timestep */
  for(i = 0;i<3;i++){
    if(!xdr_int(xdridptr[fp+1], &(i_inp[i]))){
      gmx_fseek(xdrfiles[fp+1],off+sizeof(int),SEEK_SET);
      return -1;
    }    
  }
  /* quick return */
  if(i_inp[0] != XTC_MAGIC){
    if(gmx_fseek(xdrfiles[fp+1],off+sizeof(int),SEEK_SET)){
      return -1;
    }
    return 0;
  }
  /* read time and box */
  for(i = 0;i<10;i++){
    if(!xdr_float(xdridptr[fp+1], &(f_inp[i]))){
      gmx_fseek(xdrfiles[fp+1],off+sizeof(int),SEEK_SET);
      return -1;
    }    
  }
  /* Make a rigourous check to see if we are in the beggining of a header
     Hopefully there are no ambiguous cases */
  /* This check makes use of the fact that the box matrix has 3 zeroes on the upper
     right triangle and that the first element must be nonzero unless the entire matrix is zero
  */
  if(i_inp[1] == natoms && 
     ((f_inp[1] != 0 && f_inp[6] == 0) ||
      (f_inp[1] == 0 && f_inp[5] == 0 && f_inp[9] == 0))){
    if(gmx_fseek(xdrfiles[fp+1],off+sizeof(int),SEEK_SET)){
      return -1;
    }
    *time = f_inp[0];
    *timestep = i_inp[2];
    return 1;
  }
  if(gmx_fseek(xdrfiles[fp+1],off+sizeof(int),SEEK_SET)){
    return -1;
  }
  return 0;         
}

static int 
xtc_get_next_frame_number(int fp,int natoms)
{
    gmx_off_t off;
    int step;  
    float time;
    int ret;

    if((off = gmx_ftell(xdrfiles[fp+1])) < 0){
      return -1;
    }

    /* read one int just to make sure we dont read this frame but the next */
    xdr_int(xdridptr[fp+1],&step);
    while(1){
      ret = xtc_at_header_start(fp,natoms,&step,&time);
      if(ret == 1){
    if(gmx_fseek(xdrfiles[fp+1],off,SEEK_SET)){
      return -1;
    }
    return step;
      }else if(ret == -1){
    if(gmx_fseek(xdrfiles[fp+1],off,SEEK_SET)){
      return -1;
    }
      }
    }
    return -1;
}


static float 
xtc_get_next_frame_time(int fp,int natoms, bool * bOK)
{
    gmx_off_t off;
    float time;
    int step;
    int ret;
    *bOK = 0;
    
    if((off = gmx_ftell(xdrfiles[fp+1])) < 0){
      return -1;
    }
    /* read one int just to make sure we dont read this frame but the next */
    xdr_int(xdridptr[fp+1],&step);
    while(1){
      ret = xtc_at_header_start(fp,natoms,&step,&time);
      if(ret == 1){
    *bOK = 1;
    if(gmx_fseek(xdrfiles[fp+1],off,SEEK_SET)){
      *bOK = 0;
      return -1;
    }
    return time;
      }else if(ret == -1){
    if(gmx_fseek(xdrfiles[fp+1],off,SEEK_SET)){
      return -1;
    }
    return -1;
      }
    }
    return -1;
}


static float 
xtc_get_current_frame_time(int fp,int natoms, bool * bOK)
{
    gmx_off_t off;
    int step;  
    float time;
    int ret;
    *bOK = 0;

    if((off = gmx_ftell(xdrfiles[fp+1])) < 0){
      return -1;
    }

    
    while(1){
      ret = xtc_at_header_start(fp,natoms,&step,&time);
      if(ret == 1){
    *bOK = 1;
    if(gmx_fseek(xdrfiles[fp+1],off,SEEK_SET)){
      *bOK = 0;
      return -1;
    }
    return time;
      }else if(ret == -1){
    if(gmx_fseek(xdrfiles[fp+1],off,SEEK_SET)){
      return -1;
    }
    return -1;
      }else if(ret == 0){
    /*Go back.*/
    if(gmx_fseek(xdrfiles[fp+1],-8,SEEK_CUR)){
      return -1;
    }
      }
    }
    return -1;
}

/* Currently not used, just for completeness */
static int 
xtc_get_current_frame_number(int fp,int natoms, bool * bOK)
{
    gmx_off_t off;
    int ret;  
    int step;
    float time;
    *bOK = 0;
    
    if((off = gmx_ftell(xdrfiles[fp+1])) < 0){
      return -1;
    }


    while(1){
      ret = xtc_at_header_start(fp,natoms,&step,&time);
      if(ret == 1){
    *bOK = 1;
    if(gmx_fseek(xdrfiles[fp+1],off,SEEK_SET)){
      return -1;
    }
    return step;
      }else if(ret == -1){
    if(gmx_fseek(xdrfiles[fp+1],off,SEEK_SET)){
      return -1;
    }
    return -1;
      }
    }
    return -1;
}


static gmx_off_t xtc_get_next_frame_start(int fp, int natoms)
{
    int inp;
    gmx_off_t res;
    int ret;
    int step;
    float time;
    /* read one int just to make sure we dont read this frame but the next */
    xdr_int(xdridptr[fp+1],&step);
    while(1)
    {
      ret = xtc_at_header_start(fp,natoms,&step,&time);
      if(ret == 1){
    if((res = gmx_ftell(xdrfiles[fp+1])) >= 0){
      return res - sizeof(int);
    }else{
      return res;
    }
      }else if(ret == -1){
    return -1;
      }
    }
    return -1;
}



static float 
xtc_estimate_dt(int fp, int natoms, bool * bOK)
{
  float  res;
  float  tinit;
  gmx_off_t off;
  
  if((off   = gmx_ftell(xdrfiles[fp+1])) < 0){
    return -1;
  }
  
    tinit = xtc_get_current_frame_time(fp,natoms,bOK);
    
    *bOK = 1;
    
    if(!bOK)
    {
        return -1;
    }
    
    res = xtc_get_next_frame_time(fp,natoms,bOK);
    
    if(!bOK)
    {
        return -1;
    }
    
    res -= tinit;
    if(gmx_fseek(xdrfiles[fp+1],off,SEEK_SET)){
      *bOK = 0;
      return -1;
    }
    return res;
}


int 
xtc_seek_frame(int frame, int fp, int natoms)
{
    gmx_off_t low = 0;
    gmx_off_t high,pos;

    
    /* round to 4 bytes */
    int fr;
    gmx_off_t  offset;
    if(gmx_fseek(xdrfiles[fp+1],0,SEEK_END)){
      return -1;
    }

    if((high = gmx_ftell(xdrfiles[fp+1])) < 0){
      return -1;
    }
    
    /* round to 4 bytes  */
    high /= sizeof(int);
    high *= sizeof(int);
    offset = ((high/2)/sizeof(int))*sizeof(int);
    
    if(gmx_fseek(xdrfiles[fp+1],offset,SEEK_SET)){
      return -1;
    }
    
    while(1)
    {
        fr = xtc_get_next_frame_number(fp,natoms);
        if(fr < 0)
        {
            return -1;
        }
        if(fr != frame && abs(low-high) > header_size)
        {
            if(fr < frame)
            {
                low = offset;      
            }
            else
            {
                high = offset;      
            }
            /* round to 4 bytes */
            offset = (((high+low)/2)/4)*4;
            
            if(gmx_fseek(xdrfiles[fp+1],offset,SEEK_SET)){
          return -1;
        }
        }
        else
        {
            break;
        }
    }
    if(offset <= header_size)
    {
        offset = low;
    }
    
    if(gmx_fseek(xdrfiles[fp+1],offset,SEEK_SET)){
      return -1;
    }

    if((pos = xtc_get_next_frame_start(fp,natoms))< 0){
    /* we probably hit an end of file */
      return -1;
    }
    
    if(gmx_fseek(xdrfiles[fp+1],pos,SEEK_SET)){
      return -1;
    }
    
    return 0;
}

     

int 
xtc_seek_time(real time, int fp, int natoms)
{
    float t;
    float dt;
    bool bOK;
    gmx_off_t low = 0;
    gmx_off_t high,offset,pos;
    int res;
    int dt_sign = 0;
  
    if(gmx_fseek(xdrfiles[fp+1],0,SEEK_END)){
      return -1;
    }
    
    if((high = gmx_ftell(xdrfiles[fp+1])) < 0){
      return -1;
    }
    /* round to int  */
    high /= sizeof(int);
    high *= sizeof(int);
    offset = ((high/2)/sizeof(int))*sizeof(int);

    if(gmx_fseek(xdrfiles[fp+1],offset,SEEK_SET)){
      return -1;    
   }

    dt = xtc_estimate_dt(fp,natoms,&bOK);
      
      
    
    if(!bOK)
    {
        return -1;
    }
    
    while(1)
    {
    dt = xtc_estimate_dt(fp,natoms,&bOK);
        if(!bOK)
        {
            return -1;
        }else{
      if(dt > 0){
        if(dt_sign == -1){
          /* Found a place in the trajectory that has positive time step while
         other has negative time step */
          return -2;
        }
        dt_sign = 1;
      }else if(dt < 0){
        if(dt_sign == 1){
          /* Found a place in the trajectory that has positive time step while
         other has negative time step */
          return -2;
        }
        dt_sign = -1;
      }      
    }
        t = xtc_get_next_frame_time(fp,natoms,&bOK);
        if(!bOK)
        {
            return -1;
        }

    /* If we are before the target time and the time step is positive or 0, or we have
     after the target time and the time step is negative, or the difference between 
    the current time and the target time is bigger than dt and above all the distance between high
    and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
    if we reached the solution */
        if((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) || ((t-time) >= dt && dt_sign >= 0) || ((time-t) >= -dt && dt_sign < 0))  && (abs(low-high) > header_size))
        {
      if(dt >= 0 && dt_sign != -1)
            {
                if(t < time)
                {
                    low = offset;      
                }
                else
                {
                    high = offset;      
                }
        }
      else if(dt <= 0 && dt_sign == -1)
            {
          if(t >= time)
                {
          low = offset;      
                }
          else
                {
          high = offset;      
                }
        }else{
          /* We should never reach here */
          return -1;
        }
            /* round to 4 bytes and subtract header*/
            offset = (((high+low)/2)/sizeof(int))*sizeof(int);
            if(gmx_fseek(xdrfiles[fp+1],offset,SEEK_SET)){
          return -1;
        }
        }
        else
        {
            if(abs(low-high) <= header_size)
            {
                break;
            }
            /* reestimate dt */
            if(xtc_estimate_dt(fp,natoms,&bOK) != dt)
            {
                if(bOK)
                {
                    dt = xtc_estimate_dt(fp,natoms,&bOK);
                }
            }
            if(t >= time && t-time < dt)
            {
                break;
            }
        }        
    }
    
    if(offset <= header_size)
    {
        offset = low;
    }
    
    gmx_fseek(xdrfiles[fp+1],offset,SEEK_SET);

    if((pos = xtc_get_next_frame_start(fp,natoms)) < 0){
      return -1;
    }
    
    if(gmx_fseek(xdrfiles[fp+1],pos,SEEK_SET)){
      return -1;
    }
    return 0;
}

float 
xtc_get_last_frame_time(int fp, int natoms, bool * bOK)
{
    float  time;
    gmx_off_t  off;
    int res;
    *bOK = 1;
    off = gmx_ftell(xdrfiles[fp+1]);  
    if(off < 0){
      *bOK = 0;
      return -1;
    }
    
    if( (res = gmx_fseek(xdrfiles[fp+1],-4,SEEK_END)) != 0){
      *bOK = 0;
      return -1;
    }

    time = xtc_get_current_frame_time(fp, natoms, bOK);
    if(!(*bOK)){
      return -1;
    }
    
    if( (res = gmx_fseek(xdrfiles[fp+1],off,SEEK_SET)) != 0){
      *bOK = 0;
      return -1;
    } 
    return time;
}


int
xtc_get_last_frame_number(int fp, int natoms, bool * bOK)
{
    int    frame;
    gmx_off_t  off;
    int res;
    *bOK = 1;
    
    if((off = gmx_ftell(xdrfiles[fp+1])) < 0){
      *bOK = 0;
      return -1;
    }

    
    if(gmx_fseek(xdrfiles[fp+1],-4,SEEK_END)){
      *bOK = 0;
      return -1;
    }

    frame = xtc_get_current_frame_number(fp, natoms, bOK);
    if(!bOK){
      return -1;
    }


    if(gmx_fseek(xdrfiles[fp+1],off,SEEK_SET)){
      *bOK = 0;
      return -1;
    }    

    return frame;
}

/*
 * This file is part of the GROMACS molecular simulation package.
 *
 * Copyright 1991- The GROMACS Authors
 * and the project initiators Erik Lindahl, Berk Hess and David van der Spoel.
 * Consult the AUTHORS/COPYING files and https://www.gromacs.org for details.
 *
 * GROMACS is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 *
 * GROMACS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with GROMACS; if not, see
 * https://www.gnu.org/licenses, or write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 *
 * If you want to redistribute modifications to GROMACS, please
 * consider that scientific software is very special. Version
 * control is crucial - bugs must be traceable. We will be happy to
 * consider code for inclusion in the official distribution, but
 * derived work must not be called official GROMACS. Details are found
 * in the README & COPYING files - if they are missing, get the
 * official version at https://www.gromacs.org.
 *
 * To help us fund GROMACS development, we humbly ask that you cite
 * the research papers on the package. Check out https://www.gromacs.org.
 */
#include "gmxpre.h"

#include "config.h"

#if GMX_INTERNAL_XDR


#    include <cstdlib>
#    include <cstring>

#    include "gmx_internal_xdr.h"


/* NB - THIS FILE IS ONLY USED ON MICROSOFT WINDOWS, since that
 * system doesn't provide any standard XDR system libraries. It will
 * most probably work on other platforms too, but make sure you
 * test that the xtc files produced are ok before using it.
 *
 * This header file contains Gromacs versions of the definitions for
 * Sun External Data Representation (XDR) headers and routines.
 *
 * On most UNIX systems this is already present as part of your
 * system libraries, but since we want to make Gromacs portable to
 * platforms like Microsoft Windows we have created a private version
 * of the necessary routines and distribute them with the Gromacs source.
 *
 * Although the rest of Gromacs is LGPL, you can copy and use the XDR
 * routines in any way you want as long as you obey Sun's 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
 */


/*
 * for unit alignment
 */
static char xdr_zero[BYTES_PER_XDR_UNIT] = { 0, 0, 0, 0 };

static xdr_uint32_t xdr_swapbytes(xdr_uint32_t x)
{
    xdr_uint32_t y;
    int          i;
    char*        px = reinterpret_cast<char*>(&x);
    char*        py = reinterpret_cast<char*>(&y);

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

    return y;
}

static xdr_uint32_t xdr_htonl(xdr_uint32_t x)
{
    short s = 0x0F00;
    if (*(reinterpret_cast<char*>(&s)) == static_cast<char>(0x0F))
    {
        /* bigendian, do nothing */
        return x;
    }
    else
    {
        /* smallendian,swap bytes */
        return xdr_swapbytes(x);
    }
}

static xdr_uint32_t xdr_ntohl(xdr_uint32_t x)
{
    short s = 0x0F00;
    if (*(reinterpret_cast<char*>(&s)) == static_cast<char>(0x0F))
    {
        /* bigendian, do nothing */
        return x;
    }
    else
    {
        /* smallendian, swap bytes */
        return xdr_swapbytes(x);
    }
}


/*
 * Free a data structure using XDR
 * Not a filter, but a convenient utility nonetheless
 */
void xdr_free(xdrproc_t proc, char* objp)
{
    XDR x;

    x.x_op = XDR_FREE;
    (*proc)(&x, objp);
}

/*
 * XDR nothing
 */
bool_t xdr_void()
{
    return TRUE;
}

/*
 * XDR integers
 */
bool_t xdr_int(XDR* xdrs, int* ip)
{
    xdr_int32_t l;

    switch (xdrs->x_op)
    {
        case XDR_ENCODE: l = static_cast<xdr_int32_t>(*ip); return xdr_putint32(xdrs, &l);

        case XDR_DECODE:
            if (!xdr_getint32(xdrs, &l))
            {
                return FALSE;
            }
            *ip = static_cast<int>(l);
            return TRUE;

        case XDR_FREE: return TRUE;
    }
    return FALSE;
}


/*
 * XDR unsigned integers
 */
bool_t xdr_u_int(XDR* xdrs, unsigned int* up)
{
    xdr_uint32_t l;

    switch (xdrs->x_op)
    {
        case XDR_ENCODE: l = static_cast<xdr_uint32_t>(*up); return xdr_putuint32(xdrs, &l);

        case XDR_DECODE:
            if (!xdr_getuint32(xdrs, &l))
            {
                return FALSE;
            }
            *up = static_cast<unsigned int>(l);
            return TRUE;

        case XDR_FREE: return TRUE;
    }
    return FALSE;
}


/*
 * XDR short integers
 */
bool_t xdr_short(XDR* xdrs, short* sp)
{
    xdr_int32_t l;

    switch (xdrs->x_op)
    {
        case XDR_ENCODE: l = static_cast<xdr_int32_t>(*sp); return xdr_putint32(xdrs, &l);

        case XDR_DECODE:
            if (!xdr_getint32(xdrs, &l))
            {
                return FALSE;
            }
            *sp = static_cast<short>(l);
            return TRUE;

        case XDR_FREE: return TRUE;
    }
    return FALSE;
}


/*
 * XDR unsigned short integers
 */
bool_t xdr_u_short(XDR* xdrs, unsigned short* usp)
{
    xdr_uint32_t l;

    switch (xdrs->x_op)
    {
        case XDR_ENCODE: l = static_cast<xdr_uint32_t>(*usp); return xdr_putuint32(xdrs, &l);

        case XDR_DECODE:
            if (!xdr_getuint32(xdrs, &l))
            {
                return FALSE;
            }
            *usp = static_cast<unsigned short>(l);
            return TRUE;

        case XDR_FREE: return TRUE;
    }
    return FALSE;
}


/*
 * XDR a char
 */
bool_t xdr_char(XDR* xdrs, char* cp)
{
    int i;

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

/*
 * XDR an unsigned char
 */
bool_t xdr_u_char(XDR* xdrs, unsigned char* cp)
{
    unsigned int u;

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

/*
 * XDR booleans
 */
bool_t xdr_bool(XDR* xdrs, int* bp)
{
#    define XDR_FALSE static_cast<xdr_int32_t>(0)
#    define XDR_TRUE static_cast<xdr_int32_t>(1)

    xdr_int32_t lb;

    switch (xdrs->x_op)
    {
        case XDR_ENCODE: lb = *bp ? XDR_TRUE : XDR_FALSE; return xdr_putint32(xdrs, &lb);

        case XDR_DECODE:
            if (!xdr_getint32(xdrs, &lb))
            {
                return FALSE;
            }
            *bp = (lb == XDR_FALSE) ? FALSE : TRUE;
            return TRUE;

        case XDR_FREE: return TRUE;
    }
    return FALSE;
#    undef XDR_FALSE
#    undef XDR_TRUE
}


/*
 * 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.
 */
bool_t xdr_opaque(XDR* xdrs, char* cp, unsigned int cnt)
{
    unsigned int rndup;
    char         crud[BYTES_PER_XDR_UNIT];

    /*
     * if no data we are done
     */
    if (cnt == 0)
    {
        return TRUE;
    }

    /*
     * 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 FALSE;
            }
            if (rndup == 0)
            {
                return TRUE;
            }
            return xdr_getbytes(xdrs, crud, rndup);

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

        case XDR_FREE: return TRUE;
    }
    return FALSE;
}


/*
 * XDR null terminated ASCII strings
 * xdr_string deals with "C strings" - arrays of bytes that are
 * terminated by a NULL character.  The parameter cpp references a
 * pointer to storage; If the pointer is null, then the necessary
 * storage is allocated.  The last parameter is the max allowed length
 * of the string as specified by a protocol.
 */
bool_t 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 == nullptr)
            {
                return TRUE; /* already free */
            }
            size = std::strlen(sp);
            break;

        case XDR_ENCODE:
            if (sp == nullptr)
            {
                return FALSE;
            }
            size = std::strlen(sp);
            break;
        case XDR_DECODE: break;
    }

    if (!xdr_u_int(xdrs, &size))
    {
        return FALSE;
    }
    if (size > maxsize)
    {
        return FALSE;
    }
    nodesize = size + 1;

    /*
     * now deal with the actual bytes
     */
    switch (xdrs->x_op)
    {
        case XDR_DECODE:
            if (nodesize == 0)
            {
                return TRUE;
            }
            if (sp == nullptr)
            {
                *cpp = sp = static_cast<char*>(std::malloc(nodesize));
            }
            if (sp == nullptr)
            {
                (void)fputs("xdr_string: out of memory\n", stderr);
                return FALSE;
            }
            sp[size] = 0;
            return xdr_opaque(xdrs, sp, size);

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

        case XDR_FREE:
            free(sp);
            *cpp = nullptr;
            return TRUE;
    }
    return FALSE;
}


/* Floating-point stuff */

bool_t xdr_float(XDR* xdrs, float* fp)
{
    xdr_int32_t tmp;

    switch (xdrs->x_op)
    {

        case XDR_ENCODE:
            tmp = *(reinterpret_cast<xdr_int32_t*>(fp));
            return (xdr_putint32(xdrs, &tmp));

        case XDR_DECODE:
            if (xdr_getint32(xdrs, &tmp))
            {
                *(reinterpret_cast<xdr_int32_t*>(fp)) = tmp;
                return (TRUE);
            }

            break;

        case XDR_FREE: return (TRUE);
    }
    return (FALSE);
}


bool_t xdr_double(XDR* xdrs, double* dp)
{

    /* Windows and some other systems don't define double-precision
     * word order in the header files, so unfortunately we have
     * to calculate it!
     *
     * For thread safety, we calculate it every time: locking this would
     * be more expensive.
     */
    /*static int LSW=-1;*/ /* Least significant fp word */
    int LSW = -1;          /* Least significant fp word */


    int*        ip;
    xdr_int32_t tmp[2];

    if (LSW < 0)
    {
        double x = 0.987654321; /* Just a number */

        /* 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
         */

        unsigned char ix = *(reinterpret_cast<char*>(&x));

        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 */
        {
            printf("Error when detecting floating-point word order.\n"
                   "Do you have a non-IEEE system?\n"
                   "If possible, use the XDR libraries provided with your system,\n"
                   "instead of the GROMACS fallback XDR source.\n");
            exit(0);
        }
    }

    switch (xdrs->x_op)
    {

        case XDR_ENCODE:
            ip     = reinterpret_cast<int*>(dp);
            tmp[0] = ip[bool(LSW == 0)];
            tmp[1] = ip[LSW];
            return static_cast<bool_t>(bool(xdr_putint32(xdrs, tmp)) && bool(xdr_putint32(xdrs, tmp + 1)));

        case XDR_DECODE:
            ip = reinterpret_cast<int*>(dp);
            if (xdr_getint32(xdrs, tmp + !LSW) && xdr_getint32(xdrs, tmp + LSW))
            {
                ip[0] = tmp[0];
                ip[1] = tmp[1];
                return (TRUE);
            }

            break;

        case XDR_FREE: return (TRUE);
    }
    return (FALSE);
}


/* Array routines */

/*
 * xdr_vector():
 *
 * XDR a fixed length array. Unlike variable-length arrays,
 * the storage of fixed length arrays is static and unfreeable.
 * > basep: base of the array
 * > size: size of the array
 * > elemsize: size of each element
 * > xdr_elem: routine to XDR each element
 */
bool_t xdr_vector(XDR* xdrs, char* basep, unsigned int nelem, unsigned int elemsize, xdrproc_t xdr_elem)
{
#    define LASTUNSIGNED (static_cast<unsigned int>(0) - 1)
    unsigned int i;
    char*        elptr;

    elptr = basep;
    for (i = 0; i < nelem; i++)
    {
        if (!(*xdr_elem)(xdrs, elptr, LASTUNSIGNED))
        {
            return FALSE;
        }
        elptr += elemsize;
    }
    return TRUE;
#    undef LASTUNSIGNED
}


static bool_t       xdrstdio_getbytes(XDR* /*xdrs*/, char* /*addr*/, unsigned int /*len*/);
static bool_t       xdrstdio_putbytes(XDR* /*xdrs*/, char* /*addr*/, unsigned int /*len*/);
static unsigned int xdrstdio_getpos(XDR* /*xdrs*/);
static bool_t       xdrstdio_setpos(XDR* /*xdrs*/, unsigned int /*pos*/);
static xdr_int32_t* xdrstdio_inline(XDR* /*xdrs*/, int /*len*/);
static void         xdrstdio_destroy(XDR* /*xdrs*/);
static bool_t       xdrstdio_getint32(XDR* /*xdrs*/, xdr_int32_t* /*ip*/);
static bool_t       xdrstdio_putint32(XDR* /*xdrs*/, xdr_int32_t* /*ip*/);
static bool_t       xdrstdio_getuint32(XDR* /*xdrs*/, xdr_uint32_t* /*ip*/);
static bool_t       xdrstdio_putuint32(XDR* /*xdrs*/, xdr_uint32_t* /*ip*/);

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


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

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

static unsigned int xdrstdio_getpos(XDR* xdrs)
{
    return static_cast<int>(ftell(reinterpret_cast<FILE*>(xdrs->x_private)));
}

static bool_t xdrstdio_setpos(XDR* xdrs, unsigned int pos)
{
    return fseek(reinterpret_cast<FILE*>(xdrs->x_private), static_cast<xdr_int32_t>(pos), 0) < 0 ? FALSE
                                                                                                 : TRUE;
}

static xdr_int32_t* xdrstdio_inline(XDR* xdrs, int len)
{
    (void)xdrs;
    (void)len;
    /*
     * Must do some work to implement this: must insure
     * enough data in the underlying stdio buffer,
     * that the buffer is aligned so that we can indirect through a
     * long *, and stuff this pointer in xdrs->x_buf.  Doing
     * a fread or fwrite to a scratch buffer would defeat
     * most of the gains to be had here and require storage
     * management on this buffer, so we don't do this.
     */
    return nullptr;
}

static bool_t xdrstdio_getint32(XDR* xdrs, xdr_int32_t* ip)
{
    xdr_int32_t mycopy;

    if (fread(&mycopy, 4, 1, reinterpret_cast<FILE*>(xdrs->x_private)) != 1)
    {
        return FALSE;
    }
    *ip = xdr_ntohl(mycopy);
    return TRUE;
}

static bool_t xdrstdio_putint32(XDR* xdrs, xdr_int32_t* ip)
{
    xdr_int32_t mycopy = xdr_htonl(*ip);

    ip = &mycopy;
    if (fwrite(ip, 4, 1, reinterpret_cast<FILE*>(xdrs->x_private)) != 1)
    {
        return FALSE;
    }
    return TRUE;
}

static bool_t xdrstdio_getuint32(XDR* xdrs, xdr_uint32_t* ip)
{
    xdr_uint32_t mycopy;

    if (fread(&mycopy, 4, 1, reinterpret_cast<FILE*>(xdrs->x_private)) != 1)
    {
        return FALSE;
    }
    *ip = xdr_ntohl(mycopy);
    return TRUE;
}

static bool_t xdrstdio_putuint32(XDR* xdrs, xdr_uint32_t* ip)
{
    xdr_uint32_t mycopy = xdr_htonl(*ip);

    ip = &mycopy;
    if (fwrite(ip, 4, 1, reinterpret_cast<FILE*>(xdrs->x_private)) != 1)
    {
        return FALSE;
    }
    return TRUE;
}

/*
 * Ops vector for stdio type XDR
 */
static struct XDR::xdr_ops xdrstdio_ops = {
    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_inline,    /* prime stream for inline macros */
    xdrstdio_destroy,   /* destroy stream */
    xdrstdio_getint32,  /* deserialize a int */
    xdrstdio_putint32,  /* serialize a int */
    xdrstdio_getuint32, /* deserialize a int */
    xdrstdio_putuint32  /* serialize a int */
};

/*
 * Initialize a stdio xdr stream.
 * Sets the xdr stream handle xdrs for use on the stream file.
 * Operation flag is set to op.
 */
void xdrstdio_create(XDR* xdrs, FILE* file, enum xdr_op op)
{
    xdrs->x_op      = op;
    xdrs->x_ops     = &xdrstdio_ops;
    xdrs->x_private = reinterpret_cast<char*>(file);
    xdrs->x_handy   = 0;
    xdrs->x_base    = nullptr;
}
#endif /* GMX_INTERNAL_XDR */

/*
 * This file is part of the GROMACS molecular simulation package.
 *
 * Copyright 1991- The GROMACS Authors
 * and the project initiators Erik Lindahl, Berk Hess and David van der Spoel.
 * Consult the AUTHORS/COPYING files and https://www.gromacs.org for details.
 *
 * GROMACS is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 *
 * GROMACS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with GROMACS; if not, see
 * https://www.gnu.org/licenses, or write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 *
 * If you want to redistribute modifications to GROMACS, please
 * consider that scientific software is very special. Version
 * control is crucial - bugs must be traceable. We will be happy to
 * consider code for inclusion in the official distribution, but
 * derived work must not be called official GROMACS. Details are found
 * in the README & COPYING files - if they are missing, get the
 * official version at https://www.gromacs.org.
 *
 * To help us fund GROMACS development, we humbly ask that you cite
 * the research papers on the package. Check out https://www.gromacs.org.
 */
#include "gmxpre.h"

#include "xtcio.h"

#include <cstring>

#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/gmxfio_xdr.h"
#include "gromacs/fileio/xdrf.h"
#include "gromacs/math/vec.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
#include "gromacs/utility/smalloc.h"

static int xdr_r2f(XDR* xdrs, real* r, gmx_bool gmx_unused bRead)
{
#if GMX_DOUBLE
    float f;
    int   ret;

    if (!bRead)
    {
        f = *r;
    }
    ret = xdr_float(xdrs, &f);
    if (bRead)
    {
        *r = f;
    }

    return ret;
#else
    return xdr_float(xdrs, static_cast<float*>(r));
#endif
}


t_fileio* open_xtc(const std::filesystem::path& fn, const char* mode)
{
    return gmx_fio_open(fn, mode);
}

void close_xtc(t_fileio* fio)
{
    gmx_fio_close(fio);
}

static void check_xtc_magic(int magic)
{
    if (magic != XTC_MAGIC && magic != XTC_NEW_MAGIC)
    {
        gmx_fatal(FARGS, "Magic Number Error in XTC file (read %d, should be %d or %d)", magic, XTC_MAGIC, XTC_NEW_MAGIC);
    }
}

static int xtc_check(const char* str, gmx_bool bResult, const char* file, int line)
{
    if (!bResult)
    {
        if (debug)
        {
            fprintf(debug,
                    "\nXTC error: read/write of %s failed, "
                    "source file %s, line %d\n",
                    str,
                    file,
                    line);
        }
        return 0;
    }
    return 1;
}

#define XTC_CHECK(s, b) xtc_check(s, b, __FILE__, __LINE__)

static int xtc_header(XDR* xd, int* magic, int* natoms, int64_t* step, real* time, gmx_bool bRead, gmx_bool* bOK)
{
    int result;

    if (xdr_int(xd, magic) == 0)
    {
        return 0;
    }
    result = XTC_CHECK("natoms", xdr_int(xd, natoms)); /* number of atoms */
    if (result)
    {
        /* Note that XTC wasn't defined to be extensible, so we can't
         * fix the fact that we used xdr_int for the step number,
         * which is defined to be signed and 32 bit. */
        int intStep = *step;
        result      = XTC_CHECK("step", xdr_int(xd, &intStep)); /* frame number    */
        *step       = intStep;
    }
    if (result)
    {
        result = XTC_CHECK("time", xdr_r2f(xd, time, bRead)); /* time */
    }
    *bOK = (result != 0);

    return result;
}

static int xtc_coord(XDR* xd, int* natoms, rvec* box, rvec* x, real* prec, int magic_number, gmx_bool bRead)
{
    int i, j, result;
#if GMX_DOUBLE
    float* ftmp;
    float  fprec;
#endif

    /* box */
    result = 1;
    for (i = 0; ((i < DIM) && result); i++)
    {
        for (j = 0; ((j < DIM) && result); j++)
        {
            result = XTC_CHECK("box", xdr_r2f(xd, &(box[i][j]), bRead));
        }
    }

    if (!result)
    {
        return result;
    }

#if GMX_DOUBLE
    /* allocate temp. single-precision array */
    snew(ftmp, static_cast<std::size_t>(*natoms) * DIM);

    /* Copy data to temp. array if writing */
    if (!bRead)
    {
        for (i = 0; (i < *natoms); i++)
        {
            ftmp[DIM * i + XX] = x[i][XX];
            ftmp[DIM * i + YY] = x[i][YY];
            ftmp[DIM * i + ZZ] = x[i][ZZ];
        }
        fprec = *prec;
    }
    result = XTC_CHECK("x", xdr3dfcoord(xd, ftmp, natoms, &fprec, magic_number));

    /* Copy from temp. array if reading */
    if (bRead)
    {
        for (i = 0; (i < *natoms); i++)
        {
            x[i][XX] = ftmp[DIM * i + XX];
            x[i][YY] = ftmp[DIM * i + YY];
            x[i][ZZ] = ftmp[DIM * i + ZZ];
        }
        *prec = fprec;
    }
    sfree(ftmp);
#else
    result = XTC_CHECK("x", xdr3dfcoord(xd, x[0], natoms, prec, magic_number));
#endif

    return result;
}


int write_xtc(t_fileio* fio, int natoms, int64_t step, real time, const rvec* box, const rvec* x, real prec)
{
    // By default we only write the new format for very large systems, but since the reading code
    // will adapt to whatever magic number is present in the header you could generate frames
    // for small systems that use the new format (which is useful for testing), and those should
    // be readable by normal implementations no matter how many atoms are present in the file.
    int      magic_number = (natoms > XTC_1995_MAX_NATOMS) ? XTC_NEW_MAGIC : XTC_MAGIC;
    XDR*     xd;
    gmx_bool bDum;
    int      bOK;

    if (!fio)
    {
        /* This means the fio object is not being used, e.g. because
           we are actually writing TNG output. We still have to return
           a pseudo-success value, to keep some callers happy. */
        return 1;
    }

    xd = gmx_fio_getxdr(fio);
    /* write magic number and xtc identidier */
    if (xtc_header(xd, &magic_number, &natoms, &step, &time, FALSE, &bDum) == 0)
    {
        return 0;
    }

    /* write data */
    bOK = xtc_coord(xd, &natoms, const_cast<rvec*>(box), const_cast<rvec*>(x), &prec, magic_number, FALSE); /* bOK will be 1 if writing went well */

    if (bOK)
    {
        if (gmx_fio_flush(fio) != 0)
        {
            bOK = 0;
        }
    }
    return bOK; /* 0 if bad, 1 if writing went well */
}

int read_first_xtc(t_fileio* fio, int* natoms, int64_t* step, real* time, matrix box, rvec** x, real* prec, gmx_bool* bOK)
{
    int  magic;
    XDR* xd;

    *bOK = TRUE;
    xd   = gmx_fio_getxdr(fio);

    /* read header and malloc x */
    if (!xtc_header(xd, &magic, natoms, step, time, TRUE, bOK))
    {
        return 0;
    }

    /* Check magic number */
    check_xtc_magic(magic);

    snew(*x, *natoms);

    *bOK = (xtc_coord(xd, natoms, box, *x, prec, magic, TRUE) != 0);

    return static_cast<int>(*bOK);
}

int read_next_xtc(t_fileio* fio, int natoms, int64_t* step, real* time, matrix box, rvec* x, real* prec, gmx_bool* bOK)
{
    int  magic;
    int  n;
    XDR* xd;

    *bOK = TRUE;
    xd   = gmx_fio_getxdr(fio);

    /* read header */
    if (!xtc_header(xd, &magic, &n, step, time, TRUE, bOK))
    {
        return 0;
    }

    /* Check magic number */
    check_xtc_magic(magic);

    if (n > natoms)
    {
        gmx_fatal(FARGS, "Frame contains more atoms (%d) than expected (%d)", n, natoms);
    }

    *bOK = (xtc_coord(xd, &natoms, box, x, prec, magic, TRUE) != 0);

    return static_cast<int>(*bOK);
}

int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
    

    static int *ip = NULL;
    static int oldsize;
    static int *buf;

    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 flag, k;
    int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun;
    float *lfp, lf;
    int tmp, *thiscoord,  prevcoord[3];
    unsigned int tmpcoord[30];

    int bufsize, xdrid, lsize;
    unsigned int bitsize;
    float inv_precision;
    int errval = 1;

    bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0;
    prevcoord[0]  = prevcoord[1]  = prevcoord[2]  = 0;
    
    /* find out if xdrs is opened for reading or for writing */
    xdrid = 0;
    while (xdridptr[xdrid] != xdrs) {
    xdrid++;
    if (xdrid >= MAXID) {
        fprintf(stderr, "xdr error. no open xdr stream\n");
        exit (1);
    }
    }
    if ((xdrmodes[xdrid] == 'w') || (xdrmodes[xdrid] == 'a')) {
    /* xdrs is open for writing – code skipped */

    } else {
    
    /* xdrs is open for reading */
    
    if (xdr_int(xdrs, &lsize) == 0) 
        return 0;
    if (*size != 0 && lsize != *size) {
        fprintf(stderr, "wrong number of coordinates in xdr3dfcoord; "
            "%d arg vs %d in file", *size, lsize);
    }
    *size = lsize;
    size3 = *size * 3;
    if (*size <= 9) {
        *precision = -1;
            return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3, 
                    (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
    }
    xdr_float(xdrs, precision);
    if (ip == NULL) {
        ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
        if (ip == NULL) {
        fprintf(stderr,"malloc failed\n");
        exit(1);
        }
        bufsize = size3 * 1.2;
        buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
        if (buf == NULL) {
        fprintf(stderr,"malloc failed\n");
        exit(1);
        }
        oldsize = *size;
    } else if (*size > oldsize) {
        ip = (int *)realloc(ip, (size_t)(size3 * sizeof(*ip)));
        if (ip == NULL) {
        fprintf(stderr,"malloc failed\n");
        exit(1);
        }
        bufsize = size3 * 1.2;
        buf = (int *)realloc(buf, (size_t)(bufsize * sizeof(*buf)));
        if (buf == NULL) {
        fprintf(stderr,"malloc failed\n");
        exit(1);
        }
        oldsize = *size;
    }
    buf[0] = buf[1] = buf[2] = 0;
    
    xdr_int(xdrs, &(minint[0]));
    xdr_int(xdrs, &(minint[1]));
    xdr_int(xdrs, &(minint[2]));

    xdr_int(xdrs, &(maxint[0]));
    xdr_int(xdrs, &(maxint[1]));
    xdr_int(xdrs, &(maxint[2]));
        
    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 (xdr_int(xdrs, &smallidx) == 0)    
        return 0;
    maxidx = MIN(LASTIDX, smallidx + 8) ;
    minidx = maxidx - 8; /* often this equal smallidx */
    smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
    smallnum = magicints[smallidx] / 2;
    sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
    larger = magicints[maxidx];

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

    if (xdr_int(xdrs, &(buf[0])) == 0)
        return 0;
        if (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]) == 0)
        return 0;
    buf[0] = buf[1] = buf[2] = 0;
    
    lfp = fp;
    inv_precision = 1.0 / * precision;
    run = 0;
    i = 0;
    lip = ip;
    }
    return 1;
}

int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
    

    static int *ip = NULL;
    static int oldsize;
    static int *buf;

    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 flag, k;
    int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun;
    float *lfp, lf;
    int tmp, *thiscoord,  prevcoord[3];
    unsigned int tmpcoord[30];

    int bufsize, xdrid, lsize;
    unsigned int bitsize;
    float inv_precision;
    int errval = 1;

    bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0;
    prevcoord[0]  = prevcoord[1]  = prevcoord[2]  = 0;
    
    /* find out if xdrs is opened for reading or for writing */
    xdrid = 0;
    while (xdridptr[xdrid] != xdrs) {
    xdrid++;
    if (xdrid >= MAXID) {
        fprintf(stderr, "xdr error. no open xdr stream\n");
        exit (1);
    }
    }
    if ((xdrmodes[xdrid] == 'w') || (xdrmodes[xdrid] == 'a')) {
    /* xdrs is open for writing – code skipped */

    } else {
    while ( i < lsize ) {
        thiscoord = (int *)(lip) + i * 3;

        if (bitsize == 0) {
        thiscoord[0] = receivebits(buf, bitsizeint[0]);
        thiscoord[1] = receivebits(buf, bitsizeint[1]);
        thiscoord[2] = receivebits(buf, bitsizeint[2]);
        } else {
        receiveints(buf, 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 = receivebits(buf, 1);
        is_smaller = 0;
        if (flag == 1) {
        run = receivebits(buf, 5);
        is_smaller = run % 3;
        run -= is_smaller;
        is_smaller--;
        }
        if (run > 0) {
        thiscoord += 3;
        for (k = 0; k < run; k+=3) {
            receiveints(buf, 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 1;
}

static int sizeofints( const 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;

}

static int receivebits(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; 
}

static void receiveints(int buf[], const 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++] = receivebits(buf, 8);
    num_of_bits -= 8;
    }
    if (num_of_bits > 0) {
    bytes[num_of_bytes++] = receivebits(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);
}

bool_t xdr_opaque(XDR* xdrs, char* cp, unsigned int cnt)
{
    unsigned int rndup;
    char         crud[BYTES_PER_XDR_UNIT];

    /*
     * if no data we are done
     */
    if (cnt == 0)
    {
        return TRUE;
    }

    /*
     * 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 FALSE;
            }
            if (rndup == 0)
            {
                return TRUE;
            }
            return xdr_getbytes(xdrs, crud, rndup);

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

        case XDR_FREE: return TRUE;
    }
    return FALSE;
}