src/external/tng_io/src/compression/tng_compress.c

   1 /* This code is part of the tng compression routines.
   2  *
   3  * Written by Daniel Spangberg
   4  * Copyright (c) 2010, 2013, The GROMACS development team.
   5  *
   6  *
   7  * This program is free software; you can redistribute it and/or
   8  * modify it under the terms of the Revised BSD License.
   9  */
  10
  11 #include <stdlib.h>
  12 #include <string.h>
  13 #include <math.h>
  14 #include "../../include/compression/tng_compress.h"
  15 #include "../../include/compression/coder.h"
  16 #include "../../include/compression/fixpoint.h"
  17
  18 /* Please see tng_compress.h for info on how to call these routines. */
  19
  20 /* This becomes TNGP for positions (little endian) and TNGV for velocities. In ASCII. */
  21 #define MAGIC_INT_POS 0x50474E54
  22 #define MAGIC_INT_VEL 0x56474E54
  23
  24 #define SPEED_DEFAULT 2 /* Default to relatively fast compression. For very good compression it makes sense to
  25                            choose speed=4 or speed=5 */
  26
  27 #define PRECISION(hi,lo) (Ptngc_i32x2_to_d(hi,lo))
  28
  29 #define MAX_FVAL 2147483647.
  30
  31 static int verify_input_data(double *x, int natoms, int nframes, double precision)
  32 {
  33   int iframe, i, j;
  34   for (iframe=0; iframe<nframes; iframe++)
  35     for (i=0; i<natoms; i++)
  36       for (j=0; j<3; j++)
  37         if (fabs(x[iframe*natoms*3+i*3+j]/precision+0.5)>=MAX_FVAL)
  38           goto error;
  39   return 0;
  40  error:
  41 #if 0
  42   for (iframe=0; iframe<nframes; iframe++)
  43     for (i=0; i<natoms; i++)
  44       for (j=0; j<3; j++)
  45         if (fabs(x[iframe*natoms*3+i*3+j]/precision+0.5)>=MAX_FVAL)
  46           printf("ERROR. Too large value: %d %d %d: %g %g %g\n",iframe,i,j,x[iframe*natoms*3+i*3+j],precision,x[iframe*natoms*3+i*3+j]/precision/MAX_FVAL);
  47 #endif
  48   return 1;
  49 }
  50
  51 static int verify_input_data_float(float *x, int natoms, int nframes, float precision)
  52 {
  53   int iframe, i, j;
  54   for (iframe=0; iframe<nframes; iframe++)
  55     for (i=0; i<natoms; i++)
  56       for (j=0; j<3; j++)
  57         if (fabs(x[iframe*natoms*3+i*3+j]/precision+0.5)>=MAX_FVAL)
  58           goto error;
  59   return 0;
  60  error:
  61 #if 0
  62   for (iframe=0; iframe<nframes; iframe++)
  63     for (i=0; i<natoms; i++)
  64       for (j=0; j<3; j++)
  65         if (fabs(x[iframe*natoms*3+i*3+j]/precision+0.5)>=MAX_FVAL)
  66           printf("ERROR. Too large value: %d %d %d: %g %g %g\n",iframe,i,j,x[iframe*natoms*3+i*3+j],precision,x[iframe*natoms*3+i*3+j]/precision/MAX_FVAL);
  67 #endif
  68   return 1;
  69 }
  70
  71 static int quantize(double *x, int natoms, int nframes,
  72                      double precision,
  73                      int *quant)
  74 {
  75   int iframe, i, j;
  76   for (iframe=0; iframe<nframes; iframe++)
  77     for (i=0; i<natoms; i++)
  78       for (j=0; j<3; j++)
  79         quant[iframe*natoms*3+i*3+j]=(int)floor((x[iframe*natoms*3+i*3+j]/precision)+0.5);
  80   return verify_input_data(x,natoms,nframes,precision);
  81 }
  82
  83 static int quantize_float(float *x, int natoms, int nframes,
  84                            float precision,
  85                            int *quant)
  86 {
  87   int iframe, i, j;
  88   for (iframe=0; iframe<nframes; iframe++)
  89     for (i=0; i<natoms; i++)
  90       for (j=0; j<3; j++)
  91         quant[iframe*natoms*3+i*3+j]=(int)floor((x[iframe*natoms*3+i*3+j]/precision)+0.5);
  92   return verify_input_data_float(x,natoms,nframes,precision);
  93 }
  94
  95 static void quant_inter_differences(int *quant, int natoms, int nframes,
  96                                     int *quant_inter)
  97 {
  98   int iframe, i, j;
  99   /* The first frame is used for absolute positions. */
 100   for (i=0; i<natoms; i++)
 101     for (j=0; j<3; j++)
 102       quant_inter[i*3+j]=quant[i*3+j];
 103   /* For all other frames, the difference to the previous frame is used. */
 104   for (iframe=1; iframe<nframes; iframe++)
 105     for (i=0; i<natoms; i++)
 106       for (j=0; j<3; j++)
 107         quant_inter[iframe*natoms*3+i*3+j]=quant[iframe*natoms*3+i*3+j]-quant[(iframe-1)*natoms*3+i*3+j];
 108 }
 109
 110 static void quant_intra_differences(int *quant, int natoms, int nframes,
 111                                     int *quant_intra)
 112 {
 113   int iframe, i, j;
 114   for (iframe=0; iframe<nframes; iframe++)
 115     {
 116       /* The first atom is used with its absolute position. */
 117       for (j=0; j<3; j++)
 118         quant_intra[iframe*natoms*3+j]=quant[iframe*natoms*3+j];
 119       /* For all other atoms the intraframe differences are computed. */
 120       for (i=1; i<natoms; i++)
 121         for (j=0; j<3; j++)
 122           quant_intra[iframe*natoms*3+i*3+j]=quant[iframe*natoms*3+i*3+j]-quant[iframe*natoms*3+(i-1)*3+j];
 123     }
 124 }
 125
 126 static void unquantize(double *x, int natoms, int nframes,
 127                        double precision,
 128                        int *quant)
 129 {
 130   int iframe, i, j;
 131   for (iframe=0; iframe<nframes; iframe++)
 132     for (i=0; i<natoms; i++)
 133       for (j=0; j<3; j++)
 134         x[iframe*natoms*3+i*3+j]=(double)quant[iframe*natoms*3+i*3+j]*precision;
 135 }
 136
 137 static void unquantize_float(float *x, int natoms, int nframes,
 138                              float precision,
 139                              int *quant)
 140 {
 141   int iframe, i, j;
 142   for (iframe=0; iframe<nframes; iframe++)
 143     for (i=0; i<natoms; i++)
 144       for (j=0; j<3; j++)
 145         x[iframe*natoms*3+i*3+j]=(float)quant[iframe*natoms*3+i*3+j]*precision;
 146 }
 147
 148 static void unquantize_inter_differences(double *x, int natoms, int nframes,
 149                                          double precision,
 150                                          int *quant)
 151 {
 152   int iframe, i, j;
 153   for (i=0; i<natoms; i++)
 154     for (j=0; j<3; j++)
 155       {
 156         int q=quant[i*3+j]; /* First value. */
 157         x[i*3+j]=(double)q*precision;
 158         for (iframe=1; iframe<nframes; iframe++)
 159           {
 160             q+=quant[iframe*natoms*3+i*3+j];
 161             x[iframe*natoms*3+i*3+j]=(double)q*precision;
 162           }
 163       }
 164 }
 165
 166 static void unquantize_inter_differences_float(float *x, int natoms, int nframes,
 167                                                float precision,
 168                                                int *quant)
 169 {
 170   int iframe, i, j;
 171   for (i=0; i<natoms; i++)
 172     for (j=0; j<3; j++)
 173       {
 174         int q=quant[i*3+j]; /* First value. */
 175         x[i*3+j]=(float)q*precision;
 176         for (iframe=1; iframe<nframes; iframe++)
 177           {
 178             q+=quant[iframe*natoms*3+i*3+j];
 179             x[iframe*natoms*3+i*3+j]=(float)q*precision;
 180           }
 181       }
 182 }
 183
 184 static void unquantize_inter_differences_int(int *x, int natoms, int nframes,
 185                                              int *quant)
 186 {
 187   int iframe, i, j;
 188   for (i=0; i<natoms; i++)
 189     for (j=0; j<3; j++)
 190       {
 191         int q=quant[i*3+j]; /* First value. */
 192         x[i*3+j]=q;
 193         for (iframe=1; iframe<nframes; iframe++)
 194           {
 195             q+=quant[iframe*natoms*3+i*3+j];
 196             x[iframe*natoms*3+i*3+j]=q;
 197           }
 198       }
 199 }
 200
 201 /* In frame update required for the initial frame if intra-frame
 202    compression was used. */
 203 static void unquant_intra_differences_first_frame(int *quant, int natoms)
 204 {
 205   int i,j;
 206   for (j=0; j<3; j++)
 207     {
 208       int q=quant[j];
 209       for (i=1; i<natoms; i++)
 210         {
 211           q+=quant[i*3+j];
 212           quant[i*3+j]=q;
 213         }
 214     }
 215 #if 0
 216   for (j=0; j<3; j++)
 217     for (i=0; i<natoms; i++)
 218       {
 219         printf("UQ: %d %d %d: %d\n",0,j,i,quant[i*3+j]);
 220       }
 221 #endif
 222 }
 223
 224 static void unquantize_intra_differences(double *x, int natoms, int nframes,
 225                                          double precision,
 226                                          int *quant)
 227 {
 228   int iframe, i, j;
 229 #if 0
 230   printf("UQ precision=%g\n",precision);
 231 #endif
 232   for (iframe=0; iframe<nframes; iframe++)
 233     for (j=0; j<3; j++)
 234       {
 235         int q=quant[iframe*natoms*3+j];
 236         x[iframe*natoms*3+j]=(double)q*precision;
 237         for (i=1; i<natoms; i++)
 238           {
 239             q+=quant[iframe*natoms*3+i*3+j];
 240             x[iframe*natoms*3+i*3+j]=(double)q*precision;
 241           }
 242       }
 243 }
 244
 245 static void unquantize_intra_differences_float(float *x, int natoms, int nframes,
 246                                                float precision,
 247                                                int *quant)
 248 {
 249   int iframe, i, j;
 250   for (iframe=0; iframe<nframes; iframe++)
 251     for (j=0; j<3; j++)
 252       {
 253         int q=quant[iframe*natoms*3+j];
 254         x[iframe*natoms*3+j]=(float)q*precision;
 255         for (i=1; i<natoms; i++)
 256           {
 257             q+=quant[iframe*natoms*3+i*3+j];
 258             x[iframe*natoms*3+i*3+j]=(float)q*precision;
 259           }
 260       }
 261 }
 262
 263 static void unquantize_intra_differences_int(int *x, int natoms, int nframes,
 264                                              int *quant)
 265 {
 266   int iframe, i, j;
 267   for (iframe=0; iframe<nframes; iframe++)
 268     for (j=0; j<3; j++)
 269       {
 270         int q=quant[iframe*natoms*3+j];
 271         x[iframe*natoms*3+j]=q;
 272         for (i=1; i<natoms; i++)
 273           {
 274             q+=quant[iframe*natoms*3+i*3+j];
 275             x[iframe*natoms*3+i*3+j]=q;
 276           }
 277       }
 278 }
 279
 280 /* Buffer num 8 bit bytes into buffer location buf */
 281 static void bufferfix(unsigned char *buf, fix_t v, int num)
 282 {
 283   /* Store in little endian format. */
 284   unsigned char c; /* at least 8 bits. */
 285   c=(unsigned char)(v & 0xFFU);
 286   while (num--)
 287     {
 288       *buf++=c;
 289       v >>= 8;
 290       c=(unsigned char)(v & 0xFFU);
 291     }
 292 }
 293
 294 static fix_t readbufferfix(unsigned char *buf, int num)
 295 {
 296   unsigned char b;
 297   int shift=0;
 298   fix_t f=0UL;
 299   int cnt=0;
 300   do
 301     {
 302       b=buf[cnt++];
 303       f |= ((fix_t)b & 0xFF)<<shift;
 304       shift+=8;
 305     } while (--num);
 306   return f;
 307 }
 308
 309 /* Perform position compression from the quantized data. */
 310 static void compress_quantized_pos(int *quant, int *quant_inter, int *quant_intra,
 311                                    int natoms, int nframes,
 312                                    int speed,
 313                                    int initial_coding, int initial_coding_parameter,
 314                                    int coding, int coding_parameter,
 315                                    fix_t prec_hi, fix_t prec_lo,
 316                                    int *nitems,
 317                                    char *data)
 318 {
 319   int bufloc=0;
 320   char *datablock=NULL;
 321   int length=0;
 322   /* Information needed for decompression. */
 323   if (data)
 324     bufferfix((unsigned char*)data+bufloc,(fix_t)MAGIC_INT_POS,4);
 325   bufloc+=4;
 326   /* Number of atoms. */
 327   if (data)
 328     bufferfix((unsigned char*)data+bufloc,(fix_t)natoms,4);
 329   bufloc+=4;
 330   /* Number of frames. */
 331   if (data)
 332     bufferfix((unsigned char*)data+bufloc,(fix_t)nframes,4);
 333   bufloc+=4;
 334   /* Initial coding. */
 335   if (data)
 336     bufferfix((unsigned char*)data+bufloc,(fix_t)initial_coding,4);
 337   bufloc+=4;
 338   /* Initial coding parameter. */
 339   if (data)
 340     bufferfix((unsigned char*)data+bufloc,(fix_t)initial_coding_parameter,4);
 341   bufloc+=4;
 342   /* Coding. */
 343   if (data)
 344     bufferfix((unsigned char*)data+bufloc,(fix_t)coding,4);
 345   bufloc+=4;
 346   /* Coding parameter. */
 347   if (data)
 348     bufferfix((unsigned char*)data+bufloc,(fix_t)coding_parameter,4);
 349   bufloc+=4;
 350   /* Precision. */
 351   if (data)
 352     bufferfix((unsigned char*)data+bufloc,prec_lo,4);
 353   bufloc+=4;
 354   if (data)
 355     bufferfix((unsigned char*)data+bufloc,prec_hi,4);
 356   bufloc+=4;
 357   /* The initial frame */
 358   if ((initial_coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
 359       (initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE) ||
 360       (initial_coding==TNG_COMPRESS_ALGO_POS_XTC3))
 361     {
 362       struct coder *coder=Ptngc_coder_init();
 363       length=natoms*3;
 364       datablock=(char*)Ptngc_pack_array(coder,quant,&length,
 365                                             initial_coding,initial_coding_parameter,natoms,speed);
 366       Ptngc_coder_deinit(coder);
 367     }
 368   else if ((initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) ||
 369            (initial_coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
 370     {
 371       struct coder *coder=Ptngc_coder_init();
 372       length=natoms*3;
 373       datablock=(char*)Ptngc_pack_array(coder,quant_intra,&length,
 374                                             initial_coding,initial_coding_parameter,natoms,speed);
 375       Ptngc_coder_deinit(coder);
 376     }
 377   /* Block length. */
 378   if (data)
 379     bufferfix((unsigned char*)data+bufloc,(fix_t)length,4);
 380   bufloc+=4;
 381   /* The actual data block. */
 382   if (data)
 383     memcpy(data+bufloc,datablock,length);
 384   free(datablock);
 385   bufloc+=length;
 386   /* The remaining frames */
 387   if (nframes>1)
 388     {
 389       datablock=NULL;
 390       /* Inter-frame compression? */
 391       if ((coding==TNG_COMPRESS_ALGO_POS_STOPBIT_INTER) ||
 392           (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTER) ||
 393           (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTER))
 394         {
 395           struct coder *coder=Ptngc_coder_init();
 396           length=natoms*3*(nframes-1);
 397           datablock=(char*)Ptngc_pack_array(coder,quant_inter+natoms*3,&length,
 398                                                 coding,coding_parameter,natoms,speed);
 399           Ptngc_coder_deinit(coder);
 400         }
 401       /* One-to-one compression? */
 402       else if ((coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
 403                (coding==TNG_COMPRESS_ALGO_POS_XTC3) ||
 404                (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE))
 405         {
 406           struct coder *coder=Ptngc_coder_init();
 407           length=natoms*3*(nframes-1);
 408           datablock=(char*)Ptngc_pack_array(coder,quant+natoms*3,&length,
 409                                                 coding,coding_parameter,natoms,speed);
 410           Ptngc_coder_deinit(coder);
 411         }
 412       /* Intra-frame compression? */
 413       else if ((coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) ||
 414                (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
 415         {
 416           struct coder *coder=Ptngc_coder_init();
 417           length=natoms*3*(nframes-1);
 418           datablock=(char*)Ptngc_pack_array(coder,quant_intra+natoms*3,&length,
 419                                                 coding,coding_parameter,natoms,speed);
 420           Ptngc_coder_deinit(coder);
 421         }
 422       /* Block length. */
 423       if (data)
 424         bufferfix((unsigned char*)data+bufloc,(fix_t)length,4);
 425       bufloc+=4;
 426       if (datablock)
 427         {
 428           if (data)
 429             memcpy(data+bufloc,datablock,length);
 430           free(datablock);
 431         }
 432       bufloc+=length;
 433     }
 434   *nitems=bufloc;
 435 }
 436
 437 /* Perform velocity compression from vel into the data block */
 438 static void compress_quantized_vel(int *quant, int *quant_inter,
 439                                    int natoms, int nframes,
 440                                    int speed,
 441                                    int initial_coding, int initial_coding_parameter,
 442                                    int coding, int coding_parameter,
 443                                    fix_t prec_hi, fix_t prec_lo,
 444                                    int *nitems,
 445                                    char *data)
 446 {
 447   int bufloc=0;
 448   char *datablock=NULL;
 449   int length;
 450   /* Information needed for decompression. */
 451   if (data)
 452     bufferfix((unsigned char*)data+bufloc,(fix_t)MAGIC_INT_VEL,4);
 453   bufloc+=4;
 454   /* Number of atoms. */
 455   if (data)
 456     bufferfix((unsigned char*)data+bufloc,(fix_t)natoms,4);
 457   bufloc+=4;
 458   /* Number of frames. */
 459   if (data)
 460     bufferfix((unsigned char*)data+bufloc,(fix_t)nframes,4);
 461   bufloc+=4;
 462   /* Initial coding. */
 463   if (data)
 464     bufferfix((unsigned char*)data+bufloc,(fix_t)initial_coding,4);
 465   bufloc+=4;
 466   /* Initial coding parameter. */
 467   if (data)
 468     bufferfix((unsigned char*)data+bufloc,(fix_t)initial_coding_parameter,4);
 469   bufloc+=4;
 470   /* Coding. */
 471   if (data)
 472     bufferfix((unsigned char*)data+bufloc,(fix_t)coding,4);
 473   bufloc+=4;
 474   /* Coding parameter. */
 475   if (data)
 476     bufferfix((unsigned char*)data+bufloc,(fix_t)coding_parameter,4);
 477   bufloc+=4;
 478   /* Precision. */
 479   if (data)
 480     bufferfix((unsigned char*)data+bufloc,prec_lo,4);
 481   bufloc+=4;
 482   if (data)
 483     bufferfix((unsigned char*)data+bufloc,prec_hi,4);
 484   bufloc+=4;
 485
 486   length=natoms*3;
 487   /* The initial frame */
 488   if ((initial_coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) ||
 489       (initial_coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) ||
 490       (initial_coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE))
 491     {
 492       struct coder *coder=Ptngc_coder_init();
 493       datablock=(char*)Ptngc_pack_array(coder,quant,&length,
 494                                             initial_coding,initial_coding_parameter,natoms,speed);
 495       Ptngc_coder_deinit(coder);
 496     }
 497   /* Block length. */
 498   if (data)
 499     bufferfix((unsigned char*)data+bufloc,(fix_t)length,4);
 500   bufloc+=4;
 501   /* The actual data block. */
 502   if (data && datablock)
 503     {
 504       memcpy(data+bufloc,datablock,length);
 505       free(datablock);
 506       bufloc+=length;
 507     }
 508   /* The remaining frames */
 509   if (nframes>1)
 510     {
 511       datablock=NULL;
 512       /* Inter-frame compression? */
 513       if ((coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER) ||
 514           (coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER) ||
 515           (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_INTER))
 516         {
 517           struct coder *coder=Ptngc_coder_init();
 518           length=natoms*3*(nframes-1);
 519           datablock=(char*)Ptngc_pack_array(coder,quant_inter+natoms*3,&length,
 520                                                 coding,coding_parameter,natoms,speed);
 521           Ptngc_coder_deinit(coder);
 522         }
 523       /* One-to-one compression? */
 524       else if ((coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) ||
 525                (coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) ||
 526                (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE))
 527         {
 528           struct coder *coder=Ptngc_coder_init();
 529           length=natoms*3*(nframes-1);
 530           datablock=(char*)Ptngc_pack_array(coder,quant+natoms*3,&length,
 531                                                 coding,coding_parameter,natoms,speed);
 532           Ptngc_coder_deinit(coder);
 533         }
 534       /* Block length. */
 535       if (data)
 536         bufferfix((unsigned char*)data+bufloc,(fix_t)length,4);
 537       bufloc+=4;
 538       if (data)
 539         memcpy(data+bufloc,datablock,length);
 540       free(datablock);
 541       bufloc+=length;
 542     }
 543   *nitems=bufloc;
 544 }
 545
 546 static int determine_best_coding_stop_bits(struct coder *coder,int *input, int *length,
 547                                            int *coding_parameter, int natoms)
 548 {
 549   int bits;
 550   unsigned char *packed;
 551   int best_length=0;
 552   int new_parameter=-1;
 553   int io_length;
 554   for (bits=1; bits<20; bits++)
 555     {
 556       io_length=*length;
 557       packed=Ptngc_pack_array(coder,input,&io_length,
 558                                   TNG_COMPRESS_ALGO_STOPBIT,bits,natoms,0);
 559       if (packed)
 560         {
 561           if ((new_parameter==-1) || (io_length<best_length))
 562             {
 563               new_parameter=bits;
 564               best_length=io_length;
 565             }
 566           free(packed);
 567         }
 568     }
 569   if (new_parameter==-1)
 570     return 1;
 571
 572   *coding_parameter=new_parameter;
 573   *length=best_length;
 574   return 0;
 575 }
 576
 577 static int determine_best_coding_triple(struct coder *coder,int *input, int *length,
 578                                         int *coding_parameter, int natoms)
 579 {
 580   int bits;
 581   unsigned char *packed;
 582   int best_length=0;
 583   int new_parameter=-1;
 584   int io_length;
 585   for (bits=1; bits<20; bits++)
 586     {
 587       io_length=*length;
 588       packed=Ptngc_pack_array(coder,input,&io_length,
 589                                   TNG_COMPRESS_ALGO_TRIPLET,bits,natoms,0);
 590       if (packed)
 591         {
 592           if ((new_parameter==-1) || (io_length<best_length))
 593             {
 594               new_parameter=bits;
 595               best_length=io_length;
 596             }
 597           free(packed);
 598         }
 599     }
 600   if (new_parameter==-1)
 601     return 1;
 602
 603   *coding_parameter=new_parameter;
 604   *length=best_length;
 605   return 0;
 606 }
 607
 608 static void determine_best_pos_initial_coding(int *quant, int *quant_intra, int natoms, int speed,
 609                                               fix_t prec_hi, fix_t prec_lo,
 610                                               int *initial_coding, int *initial_coding_parameter)
 611 {
 612   if (*initial_coding==-1)
 613     {
 614       /* Determine all parameters automatically */
 615       int best_coding;
 616       int best_coding_parameter;
 617       int best_code_size;
 618       int current_coding;
 619       int current_coding_parameter;
 620       int current_code_size;
 621       struct coder *coder;
 622       /* Start with XTC2, it should always work. */
 623       current_coding=TNG_COMPRESS_ALGO_POS_XTC2;
 624       current_coding_parameter=0;
 625       compress_quantized_pos(quant,NULL,quant_intra,natoms,1,speed,
 626                              current_coding,current_coding_parameter,
 627                              0,0,prec_hi,prec_lo,&current_code_size,NULL);
 628       best_coding=current_coding;
 629       best_coding_parameter=current_coding_parameter;
 630       best_code_size=current_code_size;
 631
 632       /* Determine best parameter for triplet intra. */
 633       current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA;
 634       coder=Ptngc_coder_init();
 635       current_code_size=natoms*3;
 636       current_coding_parameter=0;
 637       if (!determine_best_coding_triple(coder,quant_intra,&current_code_size,&current_coding_parameter,natoms))
 638         {
 639           if (current_code_size<best_code_size)
 640             {
 641               best_coding=current_coding;
 642               best_coding_parameter=current_coding_parameter;
 643               best_code_size=current_code_size;
 644             }
 645         }
 646       Ptngc_coder_deinit(coder);
 647
 648       /* Determine best parameter for triplet one-to-one. */
 649       current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE;
 650       coder=Ptngc_coder_init();
 651       current_code_size=natoms*3;
 652       current_coding_parameter=0;
 653       if (!determine_best_coding_triple(coder,quant,&current_code_size,&current_coding_parameter,natoms))
 654         {
 655           if (current_code_size<best_code_size)
 656             {
 657               best_coding=current_coding;
 658               best_coding_parameter=current_coding_parameter;
 659               best_code_size=current_code_size;
 660             }
 661         }
 662       Ptngc_coder_deinit(coder);
 663
 664       if (speed>=2)
 665         {
 666           current_coding=TNG_COMPRESS_ALGO_POS_XTC3;
 667           current_coding_parameter=0;
 668           compress_quantized_pos(quant,NULL,quant_intra,natoms,1,speed,
 669                                  current_coding,current_coding_parameter,
 670                                  0,0,prec_hi,prec_lo,&current_code_size,NULL);
 671           if (current_code_size<best_code_size)
 672             {
 673               best_coding=current_coding;
 674               best_coding_parameter=current_coding_parameter;
 675               best_code_size=current_code_size;
 676             }
 677         }
 678       /* Test BWLZH intra */
 679       if (speed>=6)
 680         {
 681           current_coding=TNG_COMPRESS_ALGO_POS_BWLZH_INTRA;
 682           current_coding_parameter=0;
 683           compress_quantized_pos(quant,NULL,quant_intra,natoms,1,speed,
 684                                  current_coding,current_coding_parameter,
 685                                  0,0,prec_hi,prec_lo,&current_code_size,NULL);
 686           if (current_code_size<best_code_size)
 687             {
 688               best_coding=current_coding;
 689               best_coding_parameter=current_coding_parameter;
 690             }
 691         }
 692       *initial_coding=best_coding;
 693       *initial_coding_parameter=best_coding_parameter;
 694     }
 695   else
 696     {
 697       if (*initial_coding_parameter==-1)
 698         {
 699           if ((*initial_coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
 700               (*initial_coding==TNG_COMPRESS_ALGO_POS_XTC3) ||
 701               (*initial_coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
 702             *initial_coding_parameter=0;
 703           else if (*initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA)
 704             {
 705               struct coder *coder=Ptngc_coder_init();
 706               int current_code_size=natoms*3;
 707               determine_best_coding_triple(coder,quant_intra,&current_code_size,initial_coding_parameter,natoms);
 708               Ptngc_coder_deinit(coder);
 709             }
 710           else if (*initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE)
 711             {
 712               struct coder *coder=Ptngc_coder_init();
 713               int current_code_size=natoms*3;
 714               determine_best_coding_triple(coder,quant,&current_code_size,initial_coding_parameter,natoms);
 715               Ptngc_coder_deinit(coder);
 716             }
 717         }
 718     }
 719 }
 720
 721 static void determine_best_pos_coding(int *quant, int *quant_inter, int *quant_intra, int natoms, int nframes, int speed,
 722                                       fix_t prec_hi, fix_t prec_lo,
 723                                       int *coding, int *coding_parameter)
 724 {
 725   if (*coding==-1)
 726     {
 727       /* Determine all parameters automatically */
 728       int best_coding;
 729       int best_coding_parameter;
 730       int best_code_size;
 731       int current_coding;
 732       int current_coding_parameter;
 733       int current_code_size;
 734       int initial_code_size;
 735       struct coder *coder;
 736       /* Always use XTC2 for the initial coding. */
 737       compress_quantized_pos(quant,quant_inter,quant_intra,natoms,1,speed,
 738                              TNG_COMPRESS_ALGO_POS_XTC2,0,
 739                              0,0,
 740                              prec_hi,prec_lo,&initial_code_size,NULL);
 741       /* Start with XTC2, it should always work. */
 742       current_coding=TNG_COMPRESS_ALGO_POS_XTC2;
 743       current_coding_parameter=0;
 744       compress_quantized_pos(quant,quant_inter,quant_intra,natoms,nframes,speed,
 745                              TNG_COMPRESS_ALGO_POS_XTC2,0,
 746                              current_coding,current_coding_parameter,
 747                              prec_hi,prec_lo,&current_code_size,NULL);
 748       best_coding=current_coding;
 749       best_coding_parameter=current_coding_parameter;
 750       best_code_size=current_code_size-initial_code_size; /* Correct for the use of XTC2 for the first frame. */
 751
 752       /* Determine best parameter for stopbit interframe coding. */
 753       current_coding=TNG_COMPRESS_ALGO_POS_STOPBIT_INTER;
 754       coder=Ptngc_coder_init();
 755       current_code_size=natoms*3*(nframes-1);
 756       current_coding_parameter=0;
 757       if (!determine_best_coding_stop_bits(coder,quant_inter+natoms*3,&current_code_size,
 758                                            &current_coding_parameter,natoms))
 759         {
 760           if (current_code_size<best_code_size)
 761             {
 762               best_coding=current_coding;
 763               best_coding_parameter=current_coding_parameter;
 764               best_code_size=current_code_size;
 765             }
 766         }
 767       Ptngc_coder_deinit(coder);
 768
 769       /* Determine best parameter for triplet interframe coding. */
 770       current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_INTER;
 771       coder=Ptngc_coder_init();
 772       current_code_size=natoms*3*(nframes-1);
 773       current_coding_parameter=0;
 774       if (!determine_best_coding_triple(coder,quant_inter+natoms*3,&current_code_size,
 775                                         &current_coding_parameter,natoms))
 776         {
 777           if (current_code_size<best_code_size)
 778             {
 779               best_coding=current_coding;
 780               best_coding_parameter=current_coding_parameter;
 781               best_code_size=current_code_size;
 782             }
 783         }
 784       Ptngc_coder_deinit(coder);
 785
 786       /* Determine best parameter for triplet intraframe coding. */
 787       current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA;
 788       coder=Ptngc_coder_init();
 789       current_code_size=natoms*3*(nframes-1);
 790       current_coding_parameter=0;
 791       if (!determine_best_coding_triple(coder,quant_intra+natoms*3,&current_code_size,
 792                                         &current_coding_parameter,natoms))
 793         {
 794           if (current_code_size<best_code_size)
 795             {
 796               best_coding=current_coding;
 797               best_coding_parameter=current_coding_parameter;
 798               best_code_size=current_code_size;
 799             }
 800         }
 801       Ptngc_coder_deinit(coder);
 802
 803       /* Determine best parameter for triplet one-to-one coding. */
 804       current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE;
 805       coder=Ptngc_coder_init();
 806       current_code_size=natoms*3*(nframes-1);
 807       current_coding_parameter=0;
 808       if (!determine_best_coding_triple(coder,quant+natoms*3,&current_code_size,
 809                                         &current_coding_parameter,natoms))
 810         {
 811           if (current_code_size<best_code_size)
 812             {
 813               best_coding=current_coding;
 814               best_coding_parameter=current_coding_parameter;
 815               best_code_size=current_code_size;
 816             }
 817         }
 818       Ptngc_coder_deinit(coder);
 819
 820       /* Test BWLZH inter */
 821       if (speed>=4)
 822         {
 823           current_coding=TNG_COMPRESS_ALGO_POS_BWLZH_INTER;
 824           current_coding_parameter=0;
 825           compress_quantized_pos(quant,quant_inter,quant_intra,natoms,nframes,speed,
 826                                  TNG_COMPRESS_ALGO_POS_XTC2,0,
 827                                  current_coding,current_coding_parameter,
 828                                  prec_hi,prec_lo,&current_code_size,NULL);
 829           current_code_size-=initial_code_size; /* Correct for the use of XTC2 for the first frame. */
 830           if (current_code_size<best_code_size)
 831             {
 832               best_coding=current_coding;
 833               best_coding_parameter=current_coding_parameter;
 834               best_code_size=current_code_size;
 835             }
 836         }
 837
 838       /* Test BWLZH intra */
 839       if (speed>=6)
 840         {
 841           current_coding=TNG_COMPRESS_ALGO_POS_BWLZH_INTRA;
 842           current_coding_parameter=0;
 843           compress_quantized_pos(quant,quant_inter,quant_intra,natoms,nframes,speed,
 844                                  TNG_COMPRESS_ALGO_POS_XTC2,0,
 845                                  current_coding,current_coding_parameter,
 846                                  prec_hi,prec_lo,&current_code_size,NULL);
 847           current_code_size-=initial_code_size; /* Correct for the use of XTC2 for the first frame. */
 848           if (current_code_size<best_code_size)
 849             {
 850               best_coding=current_coding;
 851               best_coding_parameter=current_coding_parameter;
 852             }
 853         }
 854       *coding=best_coding;
 855       *coding_parameter=best_coding_parameter;
 856     }
 857   else if (*coding_parameter==-1)
 858     {
 859       if ((*coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
 860           (*coding==TNG_COMPRESS_ALGO_POS_XTC3) ||
 861           (*coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTER) ||
 862           (*coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
 863         *coding_parameter=0;
 864       else if (*coding==TNG_COMPRESS_ALGO_POS_STOPBIT_INTER)
 865         {
 866           struct coder *coder=Ptngc_coder_init();
 867           int current_code_size=natoms*3*(nframes-1);
 868           determine_best_coding_stop_bits(coder,quant_inter+natoms*3,&current_code_size,
 869                                           coding_parameter,natoms);
 870           Ptngc_coder_deinit(coder);
 871         }
 872       else if (*coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTER)
 873         {
 874           struct coder *coder=Ptngc_coder_init();
 875           int current_code_size=natoms*3*(nframes-1);
 876           determine_best_coding_triple(coder,quant_inter+natoms*3,&current_code_size,
 877                                        coding_parameter,natoms);
 878           Ptngc_coder_deinit(coder);
 879         }
 880       else if (*coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA)
 881         {
 882           struct coder *coder=Ptngc_coder_init();
 883           int current_code_size=natoms*3*(nframes-1);
 884           determine_best_coding_triple(coder,quant_intra+natoms*3,&current_code_size,
 885                                        coding_parameter,natoms);
 886           Ptngc_coder_deinit(coder);
 887         }
 888       else if (*coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE)
 889         {
 890           struct coder *coder=Ptngc_coder_init();
 891           int current_code_size=natoms*3*(nframes-1);
 892           determine_best_coding_triple(coder,quant+natoms*3,&current_code_size,
 893                                        coding_parameter,natoms);
 894           Ptngc_coder_deinit(coder);
 895         }
 896     }
 897 }
 898
 899 static void determine_best_vel_initial_coding(int *quant, int natoms, int speed,
 900                                               fix_t prec_hi, fix_t prec_lo,
 901                                               int *initial_coding, int *initial_coding_parameter)
 902 {
 903   if (*initial_coding==-1)
 904     {
 905       /* Determine all parameters automatically */
 906       int best_coding=-1;
 907       int best_coding_parameter=-1;
 908       int best_code_size=-1;
 909       int current_coding;
 910       int current_coding_parameter;
 911       int current_code_size;
 912       struct coder *coder;
 913       /* Start to determine best parameter for stopbit one-to-one. */
 914       current_coding=TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE;
 915       current_code_size=natoms*3;
 916       current_coding_parameter=0;
 917       coder=Ptngc_coder_init();
 918       if (!determine_best_coding_stop_bits(coder,quant,&current_code_size,
 919                                            &current_coding_parameter,natoms))
 920         {
 921           best_coding=current_coding;
 922           best_coding_parameter=current_coding_parameter;
 923           best_code_size=current_code_size;
 924         }
 925       Ptngc_coder_deinit(coder);
 926
 927       /* Determine best parameter for triplet one-to-one. */
 928       current_coding=TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE;
 929       coder=Ptngc_coder_init();
 930       current_code_size=natoms*3;
 931       current_coding_parameter=0;
 932       if (!determine_best_coding_triple(coder,quant,&current_code_size,&current_coding_parameter,natoms))
 933         {
 934           if ((best_coding==-1) || (current_code_size<best_code_size))
 935             {
 936               best_coding=current_coding;
 937               best_coding_parameter=current_coding_parameter;
 938               best_code_size=current_code_size;
 939             }
 940         }
 941       Ptngc_coder_deinit(coder);
 942
 943       /* Test BWLZH one-to-one */
 944       if (speed>=4)
 945         {
 946           current_coding=TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE;
 947           current_coding_parameter=0;
 948           compress_quantized_vel(quant,NULL,natoms,1,speed,
 949                                  current_coding,current_coding_parameter,
 950                                  0,0,prec_hi,prec_lo,&current_code_size,NULL);
 951           if ((best_coding==-1) || (current_code_size<best_code_size))
 952             {
 953               best_coding=current_coding;
 954               best_coding_parameter=current_coding_parameter;
 955             }
 956         }
 957       *initial_coding=best_coding;
 958       *initial_coding_parameter=best_coding_parameter;
 959     }
 960   else if (*initial_coding_parameter==-1)
 961     {
 962       if (*initial_coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE)
 963         *initial_coding_parameter=0;
 964       else if (*initial_coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE)
 965         {
 966           struct coder *coder=Ptngc_coder_init();
 967           int current_code_size=natoms*3;
 968           determine_best_coding_stop_bits(coder,quant,&current_code_size,
 969                                           initial_coding_parameter,natoms);
 970           Ptngc_coder_deinit(coder);
 971         }
 972       else if (*initial_coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE)
 973         {
 974           struct coder *coder=Ptngc_coder_init();
 975           int current_code_size=natoms*3;
 976           determine_best_coding_triple(coder,quant,&current_code_size,initial_coding_parameter,natoms);
 977           Ptngc_coder_deinit(coder);
 978         }
 979     }
 980 }
 981
 982 static void determine_best_vel_coding(int *quant, int *quant_inter, int natoms, int nframes, int speed,
 983                                       fix_t prec_hi, fix_t prec_lo,
 984                                       int *coding, int *coding_parameter)
 985 {
 986   if (*coding==-1)
 987     {
 988       /* Determine all parameters automatically */
 989       int best_coding;
 990       int best_coding_parameter;
 991       int best_code_size;
 992       int current_coding;
 993       int current_coding_parameter;
 994       int current_code_size;
 995       int initial_code_size;
 996       int initial_numbits=5;
 997       struct coder *coder;
 998       /* Use stopbits one-to-one coding for the initial coding. */
 999       compress_quantized_vel(quant,NULL,natoms,1,speed,
1000                              TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE,initial_numbits,
1001                              0,0,prec_hi,prec_lo,&initial_code_size,NULL);
1002
1003       /* Test stopbit one-to-one */
1004       current_coding=TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE;
1005       current_code_size=natoms*3*(nframes-1);
1006       current_coding_parameter=0;
1007       coder=Ptngc_coder_init();
1008       determine_best_coding_stop_bits(coder,quant+natoms*3,&current_code_size,
1009                                       &current_coding_parameter,natoms);
1010       Ptngc_coder_deinit(coder);
1011       best_coding=current_coding;
1012       best_code_size=current_code_size;
1013       best_coding_parameter=current_coding_parameter;
1014
1015       /* Test triplet interframe */
1016       current_coding=TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER;
1017       current_code_size=natoms*3*(nframes-1);
1018       current_coding_parameter=0;
1019       coder=Ptngc_coder_init();
1020       if (!determine_best_coding_triple(coder,quant_inter+natoms*3,&current_code_size,
1021                                         &current_coding_parameter,natoms))
1022         {
1023           if (current_code_size<best_code_size)
1024             {
1025               best_coding=current_coding;
1026               best_code_size=current_code_size;
1027               best_coding_parameter=current_coding_parameter;
1028             }
1029         }
1030       Ptngc_coder_deinit(coder);
1031
1032       /* Test triplet one-to-one */
1033       current_coding=TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE;
1034       current_code_size=natoms*3*(nframes-1);
1035       current_coding_parameter=0;
1036       coder=Ptngc_coder_init();
1037       if (!determine_best_coding_triple(coder,quant+natoms*3,&current_code_size,
1038                                         &current_coding_parameter,natoms))
1039         {
1040           if (current_code_size<best_code_size)
1041             {
1042               best_coding=current_coding;
1043               best_code_size=current_code_size;
1044               best_coding_parameter=current_coding_parameter;
1045             }
1046         }
1047       Ptngc_coder_deinit(coder);
1048
1049       /* Test stopbit interframe */
1050       current_coding=TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER;
1051       current_code_size=natoms*3*(nframes-1);
1052       current_coding_parameter=0;
1053       coder=Ptngc_coder_init();
1054       if (!determine_best_coding_stop_bits(coder,quant_inter+natoms*3,&current_code_size,
1055                                            &current_coding_parameter,natoms))
1056         {
1057           if (current_code_size<best_code_size)
1058             {
1059               best_coding=current_coding;
1060               best_code_size=current_code_size;
1061               best_coding_parameter=current_coding_parameter;
1062             }
1063         }
1064       Ptngc_coder_deinit(coder);
1065
1066       if (speed>=4)
1067         {
1068           /* Test BWLZH inter */
1069           current_coding=TNG_COMPRESS_ALGO_VEL_BWLZH_INTER;
1070           current_coding_parameter=0;
1071           compress_quantized_vel(quant,quant_inter,natoms,nframes,speed,
1072                                  TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE,initial_numbits,
1073                                  current_coding,current_coding_parameter,
1074                                  prec_hi,prec_lo,&current_code_size,NULL);
1075           current_code_size-=initial_code_size; /* Correct for the initial frame */
1076           if (current_code_size<best_code_size)
1077             {
1078               best_coding=current_coding;
1079               best_code_size=current_code_size;
1080               best_coding_parameter=current_coding_parameter;
1081             }
1082
1083           /* Test BWLZH one-to-one */
1084           current_coding=TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE;
1085           current_coding_parameter=0;
1086           compress_quantized_vel(quant,quant_inter,natoms,nframes,speed,
1087                                  TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE,initial_numbits,
1088                                  current_coding,current_coding_parameter,
1089                                  prec_hi,prec_lo,&current_code_size,NULL);
1090           current_code_size-=initial_code_size; /* Correct for the initial frame */
1091           if (current_code_size<best_code_size)
1092             {
1093               best_coding=current_coding;
1094               best_coding_parameter=current_coding_parameter;
1095             }
1096         }
1097       *coding=best_coding;
1098       *coding_parameter=best_coding_parameter;
1099     }
1100   else if (*coding_parameter==-1)
1101     {
1102       if ((*coding==TNG_COMPRESS_ALGO_VEL_BWLZH_INTER) ||
1103           (*coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE))
1104         *coding_parameter=0;
1105       else if (*coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE)
1106         {
1107           struct coder *coder=Ptngc_coder_init();
1108           int current_code_size=natoms*3*(nframes-1);
1109           determine_best_coding_stop_bits(coder,quant+natoms*3,&current_code_size,
1110                                           coding_parameter,natoms);
1111           Ptngc_coder_deinit(coder);
1112         }
1113       else if (*coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER)
1114         {
1115           struct coder *coder=Ptngc_coder_init();
1116           int current_code_size=natoms*3*(nframes-1);
1117           determine_best_coding_triple(coder,quant_inter+natoms*3,&current_code_size,
1118                                        coding_parameter,natoms);
1119           Ptngc_coder_deinit(coder);
1120         }
1121       else if (*coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE)
1122         {
1123           struct coder *coder=Ptngc_coder_init();
1124           int current_code_size=natoms*3*(nframes-1);
1125           determine_best_coding_triple(coder,quant+natoms*3,&current_code_size,
1126                                        coding_parameter,natoms);
1127           Ptngc_coder_deinit(coder);
1128         }
1129       else if (*coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER)
1130         {
1131           struct coder *coder=Ptngc_coder_init();
1132           int current_code_size=natoms*3*(nframes-1);
1133           determine_best_coding_stop_bits(coder,quant_inter+natoms*3,&current_code_size,
1134                                           coding_parameter,natoms);
1135           Ptngc_coder_deinit(coder);
1136         }
1137     }
1138 }
1139
1140 char DECLSPECDLLEXPORT *tng_compress_pos_int(int *pos, int natoms, int nframes,
1141                                              unsigned long prec_hi, unsigned long prec_lo,
1142                                              int speed,int *algo,
1143                                              int *nitems)
1144 {
1145   char *data=malloc(natoms*nframes*14+11*4); /* 12 bytes are required to store 4 32 bit integers
1146                                                This is 17% extra. The final 11*4 is to store information
1147                                                needed for decompression. */
1148   int *quant=pos; /* Already quantized positions. */
1149   int *quant_intra=malloc(natoms*nframes*3*sizeof *quant_intra);
1150   int *quant_inter=malloc(natoms*nframes*3*sizeof *quant_inter);
1151
1152   int initial_coding, initial_coding_parameter;
1153   int coding, coding_parameter;
1154   if (speed==0)
1155     speed=SPEED_DEFAULT; /* Set the default speed. */
1156   /* Boundaries of speed. */
1157   if (speed<1)
1158     speed=1;
1159   if (speed>6)
1160     speed=6;
1161   initial_coding=algo[0];
1162   initial_coding_parameter=algo[1];
1163   coding=algo[2];
1164   coding_parameter=algo[3];
1165
1166   quant_inter_differences(quant,natoms,nframes,quant_inter);
1167   quant_intra_differences(quant,natoms,nframes,quant_intra);
1168
1169   /* If any of the above codings / coding parameters are == -1, the optimal parameters must be found */
1170   if (initial_coding==-1)
1171     {
1172       initial_coding_parameter=-1;
1173       determine_best_pos_initial_coding(quant,quant_intra,natoms,speed,prec_hi,prec_lo,
1174                                         &initial_coding,&initial_coding_parameter);
1175     }
1176   else if (initial_coding_parameter==-1)
1177     {
1178       determine_best_pos_initial_coding(quant,quant_intra,natoms,speed,prec_hi,prec_lo,
1179                                         &initial_coding,&initial_coding_parameter);
1180     }
1181
1182   if (nframes==1)
1183     {
1184       coding=0;
1185       coding_parameter=0;
1186     }
1187
1188   if (nframes>1)
1189     {
1190       if (coding==-1)
1191         {
1192           coding_parameter=-1;
1193           determine_best_pos_coding(quant,quant_inter,quant_intra,natoms,nframes,speed,prec_hi,prec_lo,
1194                                     &coding,&coding_parameter);
1195         }
1196       else if (coding_parameter==-1)
1197         {
1198           determine_best_pos_coding(quant,quant_inter,quant_intra,natoms,nframes,speed,prec_hi,prec_lo,
1199                                     &coding,&coding_parameter);
1200         }
1201     }
1202
1203   compress_quantized_pos(quant,quant_inter,quant_intra,natoms,nframes,speed,
1204                          initial_coding,initial_coding_parameter,
1205                          coding,coding_parameter,
1206                          prec_hi,prec_lo,nitems,data);
1207   free(quant_inter);
1208   free(quant_intra);
1209   if (algo[0]==-1)
1210     algo[0]=initial_coding;
1211   if (algo[1]==-1)
1212     algo[1]=initial_coding_parameter;
1213   if (algo[2]==-1)
1214     algo[2]=coding;
1215   if (algo[3]==-1)
1216     algo[3]=coding_parameter;
1217   return data;
1218 }
1219
1220 char DECLSPECDLLEXPORT *tng_compress_pos(double *pos, int natoms, int nframes,
1221                                          double desired_precision,
1222                                          int speed,int *algo,
1223                                          int *nitems)
1224 {
1225   int *quant=malloc(natoms*nframes*3*sizeof *quant);
1226   char *data;
1227   fix_t prec_hi, prec_lo;
1228   Ptngc_d_to_i32x2(desired_precision,&prec_hi,&prec_lo);
1229
1230   if (quantize(pos,natoms,nframes,PRECISION(prec_hi,prec_lo),quant))
1231     data=NULL; /* Error occured. Too large input values. */
1232   else
1233     data=tng_compress_pos_int(quant,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1234   free(quant);
1235   return data;
1236 }
1237
1238 char DECLSPECDLLEXPORT *tng_compress_pos_float(float *pos, int natoms, int nframes,
1239                                                float desired_precision,
1240                                                int speed,int *algo,
1241                                                int *nitems)
1242 {
1243   int *quant=malloc(natoms*nframes*3*sizeof *quant);
1244   char *data;
1245   fix_t prec_hi, prec_lo;
1246   Ptngc_d_to_i32x2((double)desired_precision,&prec_hi,&prec_lo);
1247
1248   if (quantize_float(pos,natoms,nframes,(float)PRECISION(prec_hi,prec_lo),quant))
1249     data=NULL; /* Error occured. Too large input values. */
1250   else
1251     data=tng_compress_pos_int(quant,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1252   free(quant);
1253   return data;
1254 }
1255
1256 char DECLSPECDLLEXPORT *tng_compress_pos_find_algo(double *pos, int natoms, int nframes,
1257                                                    double desired_precision,
1258                                                    int speed,
1259                                                    int *algo,
1260                                                    int *nitems)
1261 {
1262   algo[0]=-1;
1263   algo[1]=-1;
1264   algo[2]=-1;
1265   algo[3]=-1;
1266   return tng_compress_pos(pos,natoms,nframes,desired_precision,speed,algo,nitems);
1267 }
1268
1269 char DECLSPECDLLEXPORT *tng_compress_pos_float_find_algo(float *pos, int natoms, int nframes,
1270                                                          float desired_precision,
1271                                                          int speed,
1272                                                          int *algo,
1273                                                          int *nitems)
1274 {
1275   algo[0]=-1;
1276   algo[1]=-1;
1277   algo[2]=-1;
1278   algo[3]=-1;
1279   return tng_compress_pos_float(pos,natoms,nframes,desired_precision,speed,algo,nitems);
1280 }
1281
1282 char DECLSPECDLLEXPORT *tng_compress_pos_int_find_algo(int *pos, int natoms, int nframes,
1283                                                        unsigned long prec_hi, unsigned long prec_lo,
1284                                                        int speed,int *algo,
1285                                                        int *nitems)
1286 {
1287   algo[0]=-1;
1288   algo[1]=-1;
1289   algo[2]=-1;
1290   algo[3]=-1;
1291   return tng_compress_pos_int(pos,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1292 }
1293
1294
1295
1296 int DECLSPECDLLEXPORT tng_compress_nalgo(void)
1297 {
1298   return 4; /* There are currently four parameters required:
1299
1300  1) The compression algorithm for the first frame (initial_coding).
1301  2) One parameter to the algorithm for the first frame (the initial coding parameter).
1302  3) The compression algorithm for the remaining frames (coding).
1303  4) One parameter to the algorithm for the remaining frames (the coding parameter). */
1304 }
1305
1306 char DECLSPECDLLEXPORT *tng_compress_vel_int(int *vel, int natoms, int nframes,
1307                                              unsigned long prec_hi, unsigned long prec_lo,
1308                                              int speed, int *algo,
1309                                              int *nitems)
1310 {
1311   char *data=malloc(natoms*nframes*14+11*4); /* 12 bytes are required to store 4 32 bit integers
1312                                                This is 17% extra. The final 11*4 is to store information
1313                                                needed for decompression. */
1314   int *quant=vel;
1315   int *quant_inter=malloc(natoms*nframes*3*sizeof *quant_inter);
1316
1317   int initial_coding, initial_coding_parameter;
1318   int coding, coding_parameter;
1319   if (speed==0)
1320     speed=SPEED_DEFAULT; /* Set the default speed. */
1321   /* Boundaries of speed. */
1322   if (speed<1)
1323     speed=1;
1324   if (speed>6)
1325     speed=6;
1326   initial_coding=algo[0];
1327   initial_coding_parameter=algo[1];
1328   coding=algo[2];
1329   coding_parameter=algo[3];
1330
1331   quant_inter_differences(quant,natoms,nframes,quant_inter);
1332
1333   /* If any of the above codings / coding parameters are == -1, the optimal parameters must be found */
1334   if (initial_coding==-1)
1335     {
1336       initial_coding_parameter=-1;
1337       determine_best_vel_initial_coding(quant,natoms,speed,prec_hi,prec_lo,
1338                                         &initial_coding,&initial_coding_parameter);
1339     }
1340   else if (initial_coding_parameter==-1)
1341     {
1342       determine_best_vel_initial_coding(quant,natoms,speed,prec_hi,prec_lo,
1343                                         &initial_coding,&initial_coding_parameter);
1344     }
1345
1346   if (nframes==1)
1347     {
1348       coding=0;
1349       coding_parameter=0;
1350     }
1351
1352   if (nframes>1)
1353     {
1354       if (coding==-1)
1355         {
1356           coding_parameter=-1;
1357           determine_best_vel_coding(quant,quant_inter,natoms,nframes,speed,prec_hi,prec_lo,
1358                                     &coding,&coding_parameter);
1359         }
1360       else if (coding_parameter==-1)
1361         {
1362           determine_best_vel_coding(quant,quant_inter,natoms,nframes,speed,prec_hi,prec_lo,
1363                                     &coding,&coding_parameter);
1364         }
1365     }
1366
1367   compress_quantized_vel(quant,quant_inter,natoms,nframes,speed,
1368                          initial_coding,initial_coding_parameter,
1369                          coding,coding_parameter,
1370                          prec_hi,prec_lo,nitems,data);
1371   free(quant_inter);
1372   if (algo[0]==-1)
1373     algo[0]=initial_coding;
1374   if (algo[1]==-1)
1375     algo[1]=initial_coding_parameter;
1376   if (algo[2]==-1)
1377     algo[2]=coding;
1378   if (algo[3]==-1)
1379     algo[3]=coding_parameter;
1380   return data;
1381 }
1382
1383 char DECLSPECDLLEXPORT *tng_compress_vel(double *vel, int natoms, int nframes,
1384                                          double desired_precision,
1385                                          int speed, int *algo,
1386                                          int *nitems)
1387 {
1388   int *quant=malloc(natoms*nframes*3*sizeof *quant);
1389   char *data;
1390   fix_t prec_hi, prec_lo;
1391   Ptngc_d_to_i32x2(desired_precision,&prec_hi,&prec_lo);
1392   if (quantize(vel,natoms,nframes,PRECISION(prec_hi,prec_lo),quant))
1393     data=NULL; /* Error occured. Too large input values. */
1394   else
1395     data=tng_compress_vel_int(quant,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1396   free(quant);
1397   return data;
1398 }
1399
1400 char DECLSPECDLLEXPORT *tng_compress_vel_float(float *vel, int natoms, int nframes,
1401                                                float desired_precision,
1402                                                int speed, int *algo,
1403                                                int *nitems)
1404 {
1405   int *quant=malloc(natoms*nframes*3*sizeof *quant);
1406   char *data;
1407   fix_t prec_hi, prec_lo;
1408   Ptngc_d_to_i32x2((double)desired_precision,&prec_hi,&prec_lo);
1409   if (quantize_float(vel,natoms,nframes,(float)PRECISION(prec_hi,prec_lo),quant))
1410     data=NULL; /* Error occured. Too large input values. */
1411   else
1412     data=tng_compress_vel_int(quant,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1413   free(quant);
1414   return data;
1415 }
1416
1417 char DECLSPECDLLEXPORT *tng_compress_vel_find_algo(double *vel, int natoms, int nframes,
1418                                                    double desired_precision,
1419                                                    int speed,
1420                                                    int *algo,
1421                                                    int *nitems)
1422 {
1423   algo[0]=-1;
1424   algo[1]=-1;
1425   algo[2]=-1;
1426   algo[3]=-1;
1427   return tng_compress_vel(vel,natoms,nframes,desired_precision,speed,algo,nitems);
1428 }
1429
1430 char DECLSPECDLLEXPORT *tng_compress_vel_float_find_algo(float *vel, int natoms, int nframes,
1431                                                          float desired_precision,
1432                                                          int speed,
1433                                                          int *algo,
1434                                                          int *nitems)
1435 {
1436   algo[0]=-1;
1437   algo[1]=-1;
1438   algo[2]=-1;
1439   algo[3]=-1;
1440   return tng_compress_vel_float(vel,natoms,nframes,desired_precision,speed,algo,nitems);
1441 }
1442
1443 char DECLSPECDLLEXPORT *tng_compress_vel_int_find_algo(int *vel, int natoms, int nframes,
1444                                                        unsigned long prec_hi, unsigned long prec_lo,
1445                                                        int speed,
1446                                                        int *algo,
1447                                                        int *nitems)
1448 {
1449   algo[0]=-1;
1450   algo[1]=-1;
1451   algo[2]=-1;
1452   algo[3]=-1;
1453   return tng_compress_vel_int(vel,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1454 }
1455
1456 int DECLSPECDLLEXPORT tng_compress_inquire(char *data,int *vel, int *natoms,
1457                                             int *nframes, double *precision,
1458                                             int *algo)
1459 {
1460   int bufloc=0;
1461   fix_t prec_hi, prec_lo;
1462   int initial_coding, initial_coding_parameter;
1463   int coding, coding_parameter;
1464   int magic_int;
1465   magic_int=(int)readbufferfix((unsigned char *)data+bufloc,4);
1466   bufloc+=4;
1467   if (magic_int==MAGIC_INT_POS)
1468     *vel=0;
1469   else if (magic_int==MAGIC_INT_VEL)
1470     *vel=1;
1471   else
1472     return 1;
1473   /* Number of atoms. */
1474   *natoms=(int)readbufferfix((unsigned char *)data+bufloc,4);
1475   bufloc+=4;
1476   /* Number of frames. */
1477   *nframes=(int)readbufferfix((unsigned char *)data+bufloc,4);
1478   bufloc+=4;
1479   /* Initial coding. */
1480   initial_coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1481   bufloc+=4;
1482   /* Initial coding parameter. */
1483   initial_coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1484   bufloc+=4;
1485   /* Coding. */
1486   coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1487   bufloc+=4;
1488   /* Coding parameter. */
1489   coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1490   bufloc+=4;
1491   /* Precision. */
1492   prec_lo=readbufferfix((unsigned char *)data+bufloc,4);
1493   bufloc+=4;
1494   prec_hi=readbufferfix((unsigned char *)data+bufloc,4);
1495   *precision=PRECISION(prec_hi, prec_lo);
1496   algo[0]=initial_coding;
1497   algo[1]=initial_coding_parameter;
1498   algo[2]=coding;
1499   algo[3]=coding_parameter;
1500   return 0;
1501 }
1502
1503 static int tng_compress_uncompress_pos_gen(char *data,double *posd,float *posf,int *posi,unsigned long *prec_hi, unsigned long *prec_lo)
1504 {
1505   int bufloc=0;
1506   int length;
1507   int natoms, nframes;
1508   int initial_coding, initial_coding_parameter;
1509   int coding, coding_parameter;
1510   int *quant=NULL;
1511   struct coder *coder=NULL;
1512   int rval=0;
1513   int magic_int;
1514   /* Magic integer for positions */
1515   magic_int=(int)readbufferfix((unsigned char *)data+bufloc,4);
1516   bufloc+=4;
1517   if (magic_int!=MAGIC_INT_POS)
1518     {
1519       rval=1;
1520       goto error;
1521     }
1522   /* Number of atoms. */
1523   natoms=(int)readbufferfix((unsigned char *)data+bufloc,4);
1524   bufloc+=4;
1525   /* Number of frames. */
1526   nframes=(int)readbufferfix((unsigned char *)data+bufloc,4);
1527   bufloc+=4;
1528   /* Initial coding. */
1529   initial_coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1530   bufloc+=4;
1531   /* Initial coding parameter. */
1532   initial_coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1533   bufloc+=4;
1534   /* Coding. */
1535   coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1536   bufloc+=4;
1537   /* Coding parameter. */
1538   coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1539   bufloc+=4;
1540   /* Precision. */
1541   *prec_lo=readbufferfix((unsigned char *)data+bufloc,4);
1542   bufloc+=4;
1543   *prec_hi=readbufferfix((unsigned char *)data+bufloc,4);
1544   bufloc+=4;
1545   /* Allocate the memory for the quantized positions */
1546   quant=malloc(natoms*nframes*3*sizeof *quant);
1547   /* The data block length. */
1548   length=(int)readbufferfix((unsigned char *)data+bufloc,4);
1549   bufloc+=4;
1550   /* The initial frame */
1551   coder=Ptngc_coder_init();
1552   rval=Ptngc_unpack_array(coder,(unsigned char*)data+bufloc,quant,natoms*3,
1553                          initial_coding,initial_coding_parameter,natoms);
1554   Ptngc_coder_deinit(coder);
1555   if (rval)
1556     goto error;
1557   /* Skip past the actual data block. */
1558   bufloc+=length;
1559   /* Obtain the actual positions for the initial block. */
1560   if ((initial_coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
1561       (initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE) ||
1562       (initial_coding==TNG_COMPRESS_ALGO_POS_XTC3))
1563     {
1564       if (posd)
1565         unquantize(posd,natoms,1,PRECISION(*prec_hi,*prec_lo),quant);
1566       else if (posf)
1567         unquantize_float(posf,natoms,1,(float)PRECISION(*prec_hi,*prec_lo),quant);
1568       else if (posi)
1569         memcpy(posi,quant,natoms*3*sizeof *posi);
1570     }
1571   else if ((initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) ||
1572            (initial_coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
1573     {
1574       if (posd)
1575         unquantize_intra_differences(posd,natoms,1,PRECISION(*prec_hi,*prec_lo),quant);
1576       else if (posf)
1577         unquantize_intra_differences_float(posf,natoms,1,(float)PRECISION(*prec_hi,*prec_lo),quant);
1578       else if (posi)
1579         unquantize_intra_differences_int(posi,natoms,1,quant);
1580       unquant_intra_differences_first_frame(quant,natoms);
1581     }
1582   /* The remaining frames. */
1583   if (nframes>1)
1584     {
1585       bufloc+=4;
1586       coder=Ptngc_coder_init();
1587       rval=Ptngc_unpack_array(coder,(unsigned char *)data+bufloc,quant+natoms*3,(nframes-1)*natoms*3,
1588                                   coding,coding_parameter,natoms);
1589       Ptngc_coder_deinit(coder);
1590       if (rval)
1591         goto error;
1592       if ((coding==TNG_COMPRESS_ALGO_POS_STOPBIT_INTER) ||
1593           (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTER) ||
1594           (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTER))
1595         {
1596           /* This requires that the first frame is already in one-to-one format, even if intra-frame
1597              compression was done there. Therefore the unquant_intra_differences_first_frame should be called
1598              before to convert it correctly. */
1599           if (posd)
1600             unquantize_inter_differences(posd,natoms,nframes,PRECISION(*prec_hi,*prec_lo),quant);
1601           else if (posf)
1602             unquantize_inter_differences_float(posf,natoms,nframes,(float)PRECISION(*prec_hi,*prec_lo),quant);
1603           else if (posi)
1604             unquantize_inter_differences_int(posi,natoms,nframes,quant);
1605         }
1606       else if ((coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
1607                (coding==TNG_COMPRESS_ALGO_POS_XTC3) ||
1608                (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE))
1609         {
1610           if (posd)
1611             unquantize(posd+natoms*3,natoms,nframes-1,PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1612           else if (posf)
1613             unquantize_float(posf+natoms*3,natoms,nframes-1,(float)PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1614           else if (posi)
1615             memcpy(posi+natoms*3,quant+natoms*3,natoms*3*(nframes-1)*sizeof *posi);
1616         }
1617       else if ((coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) ||
1618                (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
1619         {
1620           if (posd)
1621             unquantize_intra_differences(posd+natoms*3,natoms,nframes-1,PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1622           else if (posf)
1623             unquantize_intra_differences_float(posf+natoms*3,natoms,nframes-1,(float)PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1624           else if (posi)
1625             unquantize_intra_differences_int(posi+natoms*3,natoms,nframes-1,quant+natoms*3);
1626         }
1627     }
1628  error:
1629   free(quant);
1630   return rval;
1631 }
1632
1633 static int tng_compress_uncompress_pos(char *data,double *pos)
1634 {
1635   unsigned long prec_hi, prec_lo;
1636   return tng_compress_uncompress_pos_gen(data,pos,NULL,NULL,&prec_hi,&prec_lo);
1637 }
1638
1639 static int tng_compress_uncompress_pos_float(char *data,float *pos)
1640 {
1641   unsigned long prec_hi, prec_lo;
1642   return tng_compress_uncompress_pos_gen(data,NULL,pos,NULL,&prec_hi,&prec_lo);
1643 }
1644
1645 static int tng_compress_uncompress_pos_int(char *data,int *pos, unsigned long *prec_hi, unsigned long *prec_lo)
1646 {
1647   return tng_compress_uncompress_pos_gen(data,NULL,NULL,pos,prec_hi,prec_lo);
1648 }
1649
1650 static int tng_compress_uncompress_vel_gen(char *data,double *veld,float *velf,int *veli,unsigned long *prec_hi, unsigned long *prec_lo)
1651 {
1652   int bufloc=0;
1653   int length;
1654   int natoms, nframes;
1655   int initial_coding, initial_coding_parameter;
1656   int coding, coding_parameter;
1657   int *quant=NULL;
1658   struct coder *coder=NULL;
1659   int rval=0;
1660   int magic_int;
1661   /* Magic integer for velocities */
1662   magic_int=(int)readbufferfix((unsigned char *)data+bufloc,4);
1663   bufloc+=4;
1664   if (magic_int!=MAGIC_INT_VEL)
1665     {
1666       rval=1;
1667       goto error;
1668     }
1669   /* Number of atoms. */
1670   natoms=(int)readbufferfix((unsigned char *)data+bufloc,4);
1671   bufloc+=4;
1672   /* Number of frames. */
1673   nframes=(int)readbufferfix((unsigned char *)data+bufloc,4);
1674   bufloc+=4;
1675   /* Initial coding. */
1676   initial_coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1677   bufloc+=4;
1678   /* Initial coding parameter. */
1679   initial_coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1680   bufloc+=4;
1681   /* Coding. */
1682   coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1683   bufloc+=4;
1684   /* Coding parameter. */
1685   coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1686   bufloc+=4;
1687   /* Precision. */
1688   *prec_lo=readbufferfix((unsigned char *)data+bufloc,4);
1689   bufloc+=4;
1690   *prec_hi=readbufferfix((unsigned char *)data+bufloc,4);
1691   bufloc+=4;
1692   /* Allocate the memory for the quantized positions */
1693   quant=malloc(natoms*nframes*3*sizeof *quant);
1694   /* The data block length. */
1695   length=(int)readbufferfix((unsigned char *)data+bufloc,4);
1696   bufloc+=4;
1697   /* The initial frame */
1698   coder=Ptngc_coder_init();
1699   rval=Ptngc_unpack_array(coder,(unsigned char*)data+bufloc,quant,natoms*3,
1700                          initial_coding,initial_coding_parameter,natoms);
1701   Ptngc_coder_deinit(coder);
1702   if (rval)
1703     goto error;
1704   /* Skip past the actual data block. */
1705   bufloc+=length;
1706   /* Obtain the actual positions for the initial block. */
1707   if ((initial_coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) ||
1708       (initial_coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) ||
1709       (initial_coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE))
1710     {
1711       if (veld)
1712         unquantize(veld,natoms,1,PRECISION(*prec_hi,*prec_lo),quant);
1713       else if (velf)
1714         unquantize_float(velf,natoms,1,(float)PRECISION(*prec_hi,*prec_lo),quant);
1715       else if (veli)
1716         memcpy(veli,quant,natoms*3*sizeof *veli);
1717     }
1718   /* The remaining frames. */
1719   if (nframes>1)
1720     {
1721       bufloc+=4;
1722       coder=Ptngc_coder_init();
1723       rval=Ptngc_unpack_array(coder,(unsigned char *)data+bufloc,quant+natoms*3,(nframes-1)*natoms*3,
1724                                   coding,coding_parameter,natoms);
1725       Ptngc_coder_deinit(coder);
1726       if (rval)
1727         goto error;
1728       /* Inter-frame compression? */
1729       if ((coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER) ||
1730           (coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER) ||
1731           (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_INTER))
1732         {
1733           /* This requires that the first frame is already in one-to-one format. */
1734           if (veld)
1735             unquantize_inter_differences(veld,natoms,nframes,PRECISION(*prec_hi,*prec_lo),quant);
1736           else if (velf)
1737             unquantize_inter_differences_float(velf,natoms,nframes,(float)PRECISION(*prec_hi,*prec_lo),quant);
1738           else if (veli)
1739             unquantize_inter_differences_int(veli,natoms,nframes,quant);
1740         }
1741       /* One-to-one compression? */
1742       else if ((coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) ||
1743                (coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) ||
1744                (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE))
1745         {
1746           if (veld)
1747             unquantize(veld+natoms*3,natoms,nframes-1,PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1748           else if (velf)
1749             unquantize_float(velf+natoms*3,natoms,nframes-1,(float)PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1750           else if (veli)
1751             memcpy(veli+natoms*3,quant+natoms*3,natoms*3*(nframes-1)*sizeof *veli);
1752         }
1753     }
1754  error:
1755   free(quant);
1756   return rval;
1757 }
1758
1759 static int tng_compress_uncompress_vel(char *data,double *vel)
1760 {
1761   unsigned long prec_hi, prec_lo;
1762   return tng_compress_uncompress_vel_gen(data,vel,NULL,NULL,&prec_hi,&prec_lo);
1763 }
1764
1765 static int tng_compress_uncompress_vel_float(char *data,float *vel)
1766 {
1767   unsigned long prec_hi, prec_lo;
1768   return tng_compress_uncompress_vel_gen(data,NULL,vel,NULL,&prec_hi,&prec_lo);
1769 }
1770
1771 static int tng_compress_uncompress_vel_int(char *data,int *vel, unsigned long *prec_hi, unsigned long *prec_lo)
1772 {
1773   return tng_compress_uncompress_vel_gen(data,NULL,NULL,vel,prec_hi,prec_lo);
1774 }
1775
1776 /* Uncompresses any tng compress block, positions or velocities. It determines whether it is positions or velocities from the data buffer. The return value is 0 if ok, and 1 if not.
1777 */
1778 int DECLSPECDLLEXPORT tng_compress_uncompress(char *data,double *posvel)
1779 {
1780   int magic_int;
1781   magic_int=(int)readbufferfix((unsigned char *)data,4);
1782   if (magic_int==MAGIC_INT_POS)
1783     return tng_compress_uncompress_pos(data,posvel);
1784   else if (magic_int==MAGIC_INT_VEL)
1785     return tng_compress_uncompress_vel(data,posvel);
1786   else
1787     return 1;
1788 }
1789
1790 int DECLSPECDLLEXPORT tng_compress_uncompress_float(char *data,float *posvel)
1791 {
1792   int magic_int;
1793   magic_int=(int)readbufferfix((unsigned char *)data,4);
1794   if (magic_int==MAGIC_INT_POS)
1795     return tng_compress_uncompress_pos_float(data,posvel);
1796   else if (magic_int==MAGIC_INT_VEL)
1797     return tng_compress_uncompress_vel_float(data,posvel);
1798   else
1799     return 1;
1800 }
1801
1802 int DECLSPECDLLEXPORT tng_compress_uncompress_int(char *data,int *posvel, unsigned long *prec_hi, unsigned long *prec_lo)
1803 {
1804   int magic_int;
1805   magic_int=(int)readbufferfix((unsigned char *)data,4);
1806   if (magic_int==MAGIC_INT_POS)
1807     return tng_compress_uncompress_pos_int(data,posvel,prec_hi,prec_lo);
1808   else if (magic_int==MAGIC_INT_VEL)
1809     return tng_compress_uncompress_vel_int(data,posvel,prec_hi,prec_lo);
1810   else
1811     return 1;
1812 }
1813
1814 void DECLSPECDLLEXPORT tng_compress_int_to_double(int *posvel_int,unsigned long prec_hi, unsigned long prec_lo,
1815                                                   int natoms,int nframes,
1816                                                   double *posvel_double)
1817 {
1818   unquantize(posvel_double,natoms,nframes,PRECISION(prec_hi,prec_lo),posvel_int);
1819 }
1820
1821 void DECLSPECDLLEXPORT tng_compress_int_to_float(int *posvel_int,unsigned long prec_hi, unsigned long prec_lo,
1822                                                  int natoms,int nframes,
1823                                                  float *posvel_float)
1824 {
1825   unquantize_float(posvel_float,natoms,nframes,(float)PRECISION(prec_hi,prec_lo),posvel_int);
1826 }
1827
1828 static char *compress_algo_pos[TNG_COMPRESS_ALGO_MAX]={
1829   "Positions invalid algorithm",
1830   "Positions stopbits interframe",
1831   "Positions triplet interframe",
1832   "Positions triplet intraframe",
1833   "Positions invalid algorithm",
1834   "Positions XTC2",
1835   "Positions invalid algorithm",
1836   "Positions triplet one to one",
1837   "Positions BWLZH interframe",
1838   "Positions BWLZH intraframe",
1839   "Positions XTC3"
1840 };
1841
1842 static char *compress_algo_vel[TNG_COMPRESS_ALGO_MAX]={
1843   "Velocities invalid algorithm",
1844   "Velocities stopbits one to one",
1845   "Velocities triplet interframe",
1846   "Velocities triplet one to one",
1847   "Velocities invalid algorithm",
1848   "Velocities invalid algorithm",
1849   "Velocities stopbits interframe",
1850   "Velocities invalid algorithm",
1851   "Velocities BWLZH interframe",
1852   "Velocities BWLZH one to one",
1853   "Velocities invalid algorithm"
1854 };
1855
1856 char DECLSPECDLLEXPORT *tng_compress_initial_pos_algo(int *algo)
1857 {
1858   int i=algo[0];
1859   if (i<0)
1860     i=0;
1861   if (i>=TNG_COMPRESS_ALGO_MAX)
1862     i=0;
1863   return compress_algo_pos[i];
1864 }
1865
1866 char DECLSPECDLLEXPORT *tng_compress_pos_algo(int *algo)
1867 {
1868   int i=algo[2];
1869   if (i<0)
1870     i=0;
1871   if (i>=TNG_COMPRESS_ALGO_MAX)
1872     i=0;
1873   return compress_algo_pos[i];
1874 }
1875
1876 char DECLSPECDLLEXPORT *tng_compress_initial_vel_algo(int *algo)
1877 {
1878   int i=algo[0];
1879   if (i<0)
1880     i=0;
1881   if (i>=TNG_COMPRESS_ALGO_MAX)
1882     i=0;
1883   return compress_algo_vel[i];
1884 }
1885
1886 char DECLSPECDLLEXPORT *tng_compress_vel_algo(int *algo)
1887 {
1888   int i=algo[2];
1889   if (i<0)
1890     i=0;
1891   if (i>=TNG_COMPRESS_ALGO_MAX)
1892     i=0;
1893   return compress_algo_vel[i];
1894 }