src/external/tng_io/src/compression/xtc2.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 /* This code is heavily influenced by
  12    http://hpcv100.rc.rug.nl/xdrf.html
  13    Based on coordinate compression (c) by Frans van Hoesel.
  14    and GROMACS xtc files (http://www.gromacs.org)
  15    (c) Copyright (c) Erik Lindahl, David van der Spoel
  16 */
  17
  18 #include <stdio.h>
  19 #include <stdlib.h>
  20 #include <string.h>
  21 #include <math.h>
  22 #include "../../include/compression/coder.h"
  23 #include "../../include/compression/widemuldiv.h"
  24 #include "../../include/compression/warnmalloc.h"
  25
  26 /* Generated by gen_magic.py */
  27 #define MAX_MAGIC  92
  28
  29 static unsigned int magic[MAX_MAGIC]={
  30 2U,  3U,  4U,  5U,
  31 6U,  8U,  10U,  12U,
  32 16U,  20U,  25U,  32U,
  33 40U,  50U,  64U,  80U,
  34 101U,  128U,  161U,  203U,
  35 256U,  322U,  406U,  512U,
  36 645U,  812U,  1024U,  1290U,
  37 1625U,  2048U,  2580U,  3250U,
  38 4096U,  5160U,  6501U,  8192U,
  39 10321U,  13003U,  16384U,  20642U,
  40 26007U,  32768U,  41285U,  52015U,
  41 65536U,  82570U,  104031U,  131072U,
  42 165140U,  208063U,  262144U,  330280U,
  43 416127U,  524288U,  660561U,  832255U,
  44 1048576U,  1321122U,  1664510U,  2097152U,
  45 2642245U,  3329021U,  4194304U,  5284491U,
  46 6658042U,  8388608U,  10568983U,  13316085U,
  47 16777216U,  21137967U,  26632170U,  33554432U,
  48 42275935U,  53264340U,  67108864U,  84551870U,
  49 106528681U,  134217728U,  169103740U,  213057362U,
  50 268435456U,  338207481U,  426114725U,  536870912U,
  51 676414963U,  852229450U,  1073741824U,  1352829926U,
  52 1704458900U,  2147483648U,  2705659852U,  3408917801U,
  53 };
  54
  55 static unsigned int magic_bits[MAX_MAGIC][8]={
  56 { 3,  6,  9,  12,  15,  18,  21,  24,  },
  57 { 5,  10,  15,  20,  24,  29,  34,  39,  },
  58 { 6,  12,  18,  24,  30,  36,  42,  48,  },
  59 { 7,  14,  21,  28,  35,  42,  49,  56,  },
  60 { 8,  16,  24,  32,  39,  47,  55,  63,  },
  61 { 9,  18,  27,  36,  45,  54,  63,  72,  },
  62 { 10,  20,  30,  40,  50,  60,  70,  80,  },
  63 { 11,  22,  33,  44,  54,  65,  76,  87,  },
  64 { 12,  24,  36,  48,  60,  72,  84,  97,  },
  65 { 13,  26,  39,  52,  65,  78,  91,  104,  },
  66 { 14,  28,  42,  56,  70,  84,  98,  112,  },
  67 { 15,  30,  45,  60,  75,  90,  105,  120,  },
  68 { 16,  32,  48,  64,  80,  96,  112,  128,  },
  69 { 17,  34,  51,  68,  85,  102,  119,  136,  },
  70 { 18,  36,  54,  72,  90,  108,  127,  144,  },
  71 { 19,  38,  57,  76,  95,  114,  133,  152,  },
  72 { 20,  40,  60,  80,  100,  120,  140,  160,  },
  73 { 21,  42,  63,  84,  105,  127,  147,  168,  },
  74 { 22,  44,  66,  88,  110,  132,  154,  176,  },
  75 { 23,  46,  69,  92,  115,  138,  161,  184,  },
  76 { 24,  48,  72,  97,  120,  144,  168,  192,  },
  77 { 25,  50,  75,  100,  125,  150,  175,  200,  },
  78 { 26,  52,  78,  104,  130,  156,  182,  208,  },
  79 { 27,  54,  81,  108,  135,  162,  190,  216,  },
  80 { 28,  56,  84,  112,  140,  168,  196,  224,  },
  81 { 29,  58,  87,  116,  145,  174,  203,  232,  },
  82 { 30,  60,  90,  120,  150,  180,  210,  240,  },
  83 { 31,  62,  93,  124,  155,  186,  217,  248,  },
  84 { 32,  64,  96,  128,  160,  192,  224,  256,  },
  85 { 33,  66,  99,  132,  165,  198,  231,  264,  },
  86 { 34,  68,  102,  136,  170,  204,  238,  272,  },
  87 { 35,  70,  105,  140,  175,  210,  245,  280,  },
  88 { 36,  72,  108,  144,  180,  216,  252,  288,  },
  89 { 37,  74,  111,  148,  185,  222,  259,  296,  },
  90 { 38,  76,  114,  152,  190,  228,  266,  304,  },
  91 { 39,  78,  117,  157,  195,  234,  273,  312,  },
  92 { 40,  80,  120,  160,  200,  240,  280,  320,  },
  93 { 41,  82,  123,  164,  205,  246,  287,  328,  },
  94 { 42,  84,  127,  168,  210,  252,  294,  336,  },
  95 { 43,  86,  129,  172,  215,  258,  301,  344,  },
  96 { 44,  88,  132,  176,  220,  264,  308,  352,  },
  97 { 45,  90,  135,  180,  225,  270,  315,  360,  },
  98 { 46,  92,  138,  184,  230,  276,  322,  368,  },
  99 { 47,  94,  141,  188,  235,  282,  329,  376,  },
 100 { 48,  97,  144,  192,  240,  288,  336,  384,  },
 101 { 49,  98,  147,  196,  245,  294,  343,  392,  },
 102 { 50,  100,  150,  200,  250,  300,  350,  400,  },
 103 { 52,  102,  153,  204,  255,  306,  357,  408,  },
 104 { 52,  104,  156,  208,  260,  312,  364,  416,  },
 105 { 53,  106,  159,  212,  265,  318,  371,  424,  },
 106 { 54,  108,  162,  216,  270,  324,  378,  432,  },
 107 { 55,  110,  165,  220,  275,  330,  385,  440,  },
 108 { 56,  112,  168,  224,  280,  336,  392,  448,  },
 109 { 57,  114,  172,  228,  285,  342,  399,  456,  },
 110 { 58,  116,  174,  232,  290,  348,  406,  464,  },
 111 { 59,  118,  177,  236,  295,  354,  413,  472,  },
 112 { 60,  120,  180,  240,  300,  360,  420,  480,  },
 113 { 61,  122,  183,  244,  305,  366,  427,  488,  },
 114 { 62,  124,  186,  248,  310,  372,  434,  496,  },
 115 { 63,  127,  190,  252,  315,  378,  442,  505,  },
 116 { 64,  128,  192,  256,  320,  384,  448,  512,  },
 117 { 65,  130,  195,  260,  325,  390,  455,  520,  },
 118 { 66,  132,  198,  264,  330,  396,  462,  528,  },
 119 { 67,  134,  201,  268,  335,  402,  469,  536,  },
 120 { 68,  136,  204,  272,  340,  408,  476,  544,  },
 121 { 69,  138,  207,  276,  345,  414,  483,  552,  },
 122 { 70,  140,  210,  280,  350,  420,  490,  560,  },
 123 { 71,  142,  213,  284,  355,  426,  497,  568,  },
 124 { 72,  144,  216,  288,  360,  432,  505,  576,  },
 125 { 73,  146,  219,  292,  365,  438,  511,  584,  },
 126 { 74,  148,  222,  296,  370,  444,  518,  592,  },
 127 { 75,  150,  225,  300,  375,  451,  525,  600,  },
 128 { 76,  152,  228,  304,  380,  456,  532,  608,  },
 129 { 77,  154,  231,  308,  385,  462,  539,  616,  },
 130 { 78,  157,  234,  312,  390,  469,  546,  625,  },
 131 { 79,  158,  237,  316,  395,  474,  553,  632,  },
 132 { 80,  160,  240,  320,  400,  480,  560,  640,  },
 133 { 81,  162,  243,  324,  406,  486,  568,  648,  },
 134 { 82,  164,  246,  328,  410,  492,  574,  656,  },
 135 { 83,  166,  249,  332,  415,  498,  581,  664,  },
 136 { 84,  168,  252,  336,  420,  505,  588,  672,  },
 137 { 85,  170,  255,  340,  425,  510,  595,  680,  },
 138 { 86,  172,  258,  344,  430,  516,  602,  688,  },
 139 { 87,  174,  261,  348,  435,  522,  609,  696,  },
 140 { 88,  176,  264,  352,  440,  528,  616,  704,  },
 141 { 89,  178,  267,  356,  445,  534,  623,  712,  },
 142 { 90,  180,  270,  360,  451,  540,  631,  720,  },
 143 { 91,  182,  273,  364,  455,  546,  637,  728,  },
 144 { 92,  184,  276,  368,  460,  552,  644,  736,  },
 145 { 94,  187,  279,  373,  466,  558,  651,  745,  },
 146 { 94,  188,  282,  376,  470,  564,  658,  752,  },
 147 { 95,  190,  285,  380,  475,  570,  665,  760,  },
 148 };
 149
 150
 151 static const double iflipgaincheck=0.89089871814033927; /*  1./(2**(1./6)) */
 152
 153
 154 /* Difference in indices used for determining whether to store as large or small */
 155 #define QUITE_LARGE 3
 156 #define IS_LARGE 6
 157
 158 #if 0
 159 #define SHOWIT
 160 #endif
 161
 162 int Ptngc_magic(unsigned int i)
 163 {
 164   return magic[i];
 165 }
 166
 167 int Ptngc_find_magic_index(unsigned int maxval)
 168 {
 169   int i=0;
 170   while (magic[i]<=maxval)
 171     i++;
 172   return i;
 173 }
 174
 175 static unsigned int positive_int(int item)
 176 {
 177   int s=0;
 178   if (item>0)
 179     s=1+(item-1)*2;
 180   else if (item<0)
 181     s=2+(-item-1)*2;
 182   return s;
 183 }
 184
 185 static int unpositive_int(int val)
 186 {
 187   int s=(val+1)/2;
 188   if ((val%2)==0)
 189     s=-s;
 190   return s;
 191 }
 192
 193
 194 /* Sequence instructions */
 195 #define INSTR_DEFAULT 0
 196 #define INSTR_BASE_RUNLENGTH 1
 197 #define INSTR_ONLY_LARGE 2
 198 #define INSTR_ONLY_SMALL 3
 199 #define INSTR_LARGE_BASE_CHANGE 4
 200 #define INSTR_FLIP 5
 201 #define INSTR_LARGE_RLE 6
 202
 203 #define MAXINSTR 7
 204
 205 #ifdef SHOWIT
 206 static char *instrnames[MAXINSTR]={
 207   "large+small",
 208   "base+run",
 209   "large",
 210   "small",
 211   "large base change",
 212   "flip",
 213   "large rle",
 214 };
 215 #endif
 216
 217 /* Bit patterns in the compressed code stream: */
 218
 219 static const int seq_instr[MAXINSTR][2]=
 220   {
 221     { 1,1 }, /* 1 : one large atom + runlength encoded small integers. Use same settings as before. */
 222     { 0,2 }, /* 00 : set base and runlength in next four bits (x). base (increase/keep/decrease)=x%3-1. runlength=1+x/3.
 223                       The special value 1111 in the four bits means runlength=6 and base change=0 */
 224     { 4,4 }, /* 0100 : next only a large atom comes. */
 225     { 5,4 }, /* 0101 : next only runlength encoded small integers. Use same settings as before. */
 226     { 6,4 }, /* 0110 : Large change in base. Change is encoded in the
 227                 following 2 bits. change direction (sign) is the first
 228                 bit. The next bit +1 is the actual change. This
 229                 allows the change of up to +/- 2 indices. */
 230     { 14,5 }, /* 01110 : flip whether integers should be modified to compress water better */
 231     { 15,5 }, /* 01111 : Large rle. The next 4 bits encode how many
 232                large atoms are in the following sequence: 3-18. (2 is
 233                more efficiently coded with two large instructions. */
 234  };
 235
 236 static void write_instruction(struct coder *coder,int instr,unsigned char **output_ptr)
 237 {
 238   Ptngc_writebits(coder,seq_instr[instr][0],seq_instr[instr][1],output_ptr);
 239 #ifdef SHOWIT
 240   fprintf(stderr,"INSTR: %s (%d bits)\n",instrnames[instr],seq_instr[instr][1]);
 241 #endif
 242 }
 243
 244 static unsigned int readbits(unsigned char **ptr, int *bitptr, int nbits)
 245 {
 246   unsigned int val=0U;
 247   unsigned int extract_mask=0x80U>>*bitptr;
 248   unsigned char thisval=**ptr;
 249 #ifdef SHOWIT
 250   fprintf(stderr,"Read nbits=%d val=",nbits);
 251 #endif
 252   while (nbits--)
 253     {
 254       val<<=1;
 255       val|=((extract_mask & thisval)!=0);
 256       *bitptr=(*bitptr)+1;
 257       extract_mask>>=1;
 258       if (!extract_mask)
 259         {
 260           extract_mask=0x80U;
 261           *ptr=(*ptr)+1;
 262           *bitptr=0;
 263           if (nbits)
 264             thisval=**ptr;
 265         }
 266     }
 267 #ifdef SHOWIT
 268   fprintf(stderr,"%x\n",val);
 269 #endif
 270   return val;
 271 }
 272
 273 static void readmanybits(unsigned char **ptr, int *bitptr, int nbits, unsigned char *buffer)
 274 {
 275   while (nbits>=8)
 276     {
 277       *buffer++=readbits(ptr,bitptr,8);
 278       nbits-=8;
 279 #ifdef SHOWIT
 280       fprintf(stderr,"Read value %02x\n",buffer[-1]);
 281 #endif
 282     }
 283   if (nbits)
 284     {
 285       *buffer++=readbits(ptr,bitptr,nbits);
 286 #ifdef SHOWIT
 287       fprintf(stderr,"Read value %02x\n",buffer[-1]);
 288 #endif
 289     }
 290 }
 291
 292 static int read_instruction(unsigned char **ptr, int *bitptr)
 293 {
 294   int instr=-1;
 295   unsigned int bits=readbits(ptr,bitptr,1);
 296   if (bits)
 297     instr=INSTR_DEFAULT;
 298   else
 299     {
 300       bits=readbits(ptr,bitptr,1);
 301       if (!bits)
 302         instr=INSTR_BASE_RUNLENGTH;
 303       else
 304         {
 305           bits=readbits(ptr,bitptr,2);
 306           if (bits==0)
 307             instr=INSTR_ONLY_LARGE;
 308           else if (bits==1)
 309             instr=INSTR_ONLY_SMALL;
 310           else if (bits==2)
 311             instr=INSTR_LARGE_BASE_CHANGE;
 312           else if (bits==3)
 313             {
 314               bits=readbits(ptr,bitptr,1);
 315               if (bits==0)
 316                 instr=INSTR_FLIP;
 317               else
 318                 instr=INSTR_LARGE_RLE;
 319             }
 320         }
 321     }
 322   return instr;
 323 }
 324
 325 /* Modifies three integer values for better compression of water */
 326 static void swap_ints(int *in, int *out)
 327 {
 328   out[0]=in[0]+in[1];
 329   out[1]=-in[1];
 330   out[2]=in[1]+in[2];
 331 }
 332
 333 static void swap_is_better(int *input, int *minint, int *sum_normal, int *sum_swapped)
 334 {
 335   int normal_max=0;
 336   int swapped_max=0;
 337   int i,j;
 338   int normal[3];
 339   int swapped[3];
 340   for (i=0; i<3; i++)
 341     {
 342       normal[0]=input[i]-minint[i];
 343       normal[1]=input[3+i]-input[i]; /* minint[i]-minint[i] cancels out */
 344       normal[2]=input[6+i]-input[3+i]; /* minint[i]-minint[i] cancels out */
 345       swap_ints(normal,swapped);
 346       for (j=1; j<3; j++)
 347         {
 348           if (positive_int(normal[j])>(unsigned int)normal_max)
 349             normal_max=positive_int(normal[j]);
 350           if (positive_int(swapped[j])>(unsigned int)swapped_max)
 351             swapped_max=positive_int(swapped[j]);
 352         }
 353     }
 354   if (normal_max==0)
 355     normal_max=1;
 356   if (swapped_max==0)
 357     swapped_max=1;
 358   *sum_normal=normal_max;
 359   *sum_swapped=swapped_max;
 360 }
 361
 362 static void swapdecide(struct coder *coder, int *input,int *swapatoms, int *large_index, int *minint, unsigned char **output_ptr)
 363 {
 364   int didswap=0;
 365   int normal,swapped;
 366   (void)large_index;
 367   swap_is_better(input,minint,&normal,&swapped);
 368   /* We have to determine if it is worth to change the behaviour.
 369      If diff is positive it means that it is worth something to
 370      swap. But it costs 4 bits to do the change. If we assume that
 371      we gain 0.17 bit by the swap per value, and the runlength>2
 372      for four molecules in a row, we gain something. So check if we
 373      gain at least 0.17 bits to even attempt the swap.
 374   */
 375 #ifdef SHOWIT
 376   fprintf(stderr,"Trying Flip: %g %g\n",(double)swapped/normal, (double)normal/swapped);
 377 #endif
 378   if (((swapped<normal) && (fabs((double)swapped/normal)<iflipgaincheck)) ||
 379       ((normal<swapped) && (fabs((double)normal/swapped)<iflipgaincheck)))
 380     {
 381       if (swapped<normal)
 382         {
 383           if (!*swapatoms)
 384             {
 385               *swapatoms=1;
 386               didswap=1;
 387             }
 388         }
 389       else
 390         {
 391           if (*swapatoms)
 392             {
 393               *swapatoms=0;
 394               didswap=1;
 395             }
 396         }
 397     }
 398   if (didswap)
 399     {
 400 #ifdef SHOWIT
 401       fprintf(stderr,"Flip. Swapatoms is now %d\n",*swapatoms);
 402 #endif
 403       write_instruction(coder,INSTR_FLIP,output_ptr);
 404     }
 405 }
 406
 407 /* Compute number of bits required to store values using three different bases in the index array */
 408 static int compute_magic_bits(int *index)
 409 {
 410   unsigned int largeint[4];
 411   unsigned int largeint_tmp[4];
 412   int i,j,onebit;
 413   for (i=0; i<4; i++)
 414     largeint[i]=0U;
 415   for (i=0; i<3; i++)
 416     {
 417       if (i!=0)
 418         {
 419           /* We must do the multiplication of the largeint with the integer base */
 420           Ptngc_largeint_mul(magic[index[i]],largeint,largeint_tmp,4);
 421           for (j=0; j<4; j++)
 422             largeint[j]=largeint_tmp[j];
 423         }
 424       Ptngc_largeint_add(magic[index[i]]-1,largeint,4);
 425     }
 426   /* Find last bit. */
 427 #if 0
 428   printf("Largeint is %u %u %u\n",largeint[0],largeint[1],largeint[2]);
 429 #endif
 430   onebit=0;
 431   for (i=0; i<3; i++)
 432     for (j=0; j<32; j++)
 433       if (largeint[i]&(1U<<j))
 434         onebit=i*32+j+1;
 435   return onebit;
 436 }
 437
 438 /* Convert a sequence of (hopefully) small positive integers
 439    using the base pointed to by index (x base, y base and z base can be different).
 440    The largest number of integers supported is 18 (29 to handle/detect overflow) */
 441 static void trajcoder_base_compress(int *input, int n, int *index, unsigned char *result)
 442 {
 443   unsigned int largeint[19];
 444   unsigned int largeint_tmp[19];
 445   int i,j;
 446   for (i=0; i<19; i++)
 447     largeint[i]=0U;
 448
 449   for (i=0; i<n; i++)
 450     {
 451       if (i!=0)
 452         {
 453           /* We must do the multiplication of the largeint with the integer base */
 454           Ptngc_largeint_mul(magic[index[i%3]],largeint,largeint_tmp,19);
 455           for (j=0; j<19; j++)
 456             largeint[j]=largeint_tmp[j];
 457         }
 458       Ptngc_largeint_add(input[i],largeint,19);
 459     }
 460   if (largeint[18])
 461     {
 462       fprintf(stderr,"TRAJNG: BUG! Overflow in compression detected.\n");
 463       exit(EXIT_FAILURE);
 464     }
 465 #if 0
 466 #ifdef SHOWIT
 467   for (i=0; i<19; i++)
 468     fprintf(stderr,"Largeint[%d]=0x%x\n",i,largeint[i]);
 469 #endif
 470 #endif
 471   /* Convert the largeint to a sequence of bytes. */
 472   for (i=0; i<18; i++)
 473     {
 474       int shift=0;
 475       for (j=0; j<4; j++)
 476         {
 477           result[i*4+j]=(largeint[i]>>shift) & 0xFFU;
 478           shift+=8;
 479         }
 480     }
 481 }
 482
 483 /* The opposite of base_compress. */
 484 static void trajcoder_base_decompress(unsigned char *input, int n, int *index, int *output)
 485 {
 486   unsigned int largeint[19];
 487   unsigned int largeint_tmp[19];
 488   int i,j;
 489   /* Convert the sequence of bytes to a largeint. */
 490   for (i=0; i<18; i++)
 491     {
 492       int shift=0;
 493       largeint[i]=0U;
 494       for (j=0; j<4; j++)
 495         {
 496           largeint[i]|=((unsigned int)input[i*4+j])<<shift;
 497           shift+=8;
 498         }
 499     }
 500   largeint[18]=0U;
 501 #if 0
 502 #ifdef SHOWIT
 503   for (i=0; i<19; i++)
 504     fprintf(stderr,"Largeint[%d]=0x%x\n",i,largeint[i]);
 505 #endif
 506 #endif
 507   for (i=n-1; i>=0; i--)
 508     {
 509       unsigned int remainder=Ptngc_largeint_div(magic[index[i%3]],largeint,largeint_tmp,19);
 510 #if 0
 511 #ifdef SHOWIT
 512       fprintf(stderr,"Remainder: %u\n",remainder);
 513 #endif
 514 #endif
 515 #if 0
 516       for (j=0; j<19; j++)
 517         largeint[j]=largeint_tmp[j];
 518 #endif
 519       memcpy(largeint,largeint_tmp,19*sizeof *largeint);
 520       output[i]=remainder;
 521     }
 522 }
 523
 524 /* It is "large" if we have to increase the small index quite a
 525    bit. Not so much to be rejected by the not very large check
 526    later. */
 527 static int is_quite_large(int *input, int small_index, int max_large_index)
 528 {
 529   int is=0;
 530   int i;
 531   if (small_index+QUITE_LARGE>=max_large_index)
 532     is=1;
 533   else
 534     {
 535       for (i=0; i<3; i++)
 536         if (positive_int(input[i])>magic[small_index+QUITE_LARGE])
 537           {
 538             is=1;
 539             break;
 540           }
 541     }
 542   return is;
 543 }
 544
 545 #ifdef SHOWIT
 546 int nbits_sum;
 547 int nvalues_sum;
 548 #endif
 549
 550 static void write_three_large(struct coder *coder, int *encode_ints, int *large_index, int nbits, unsigned char *compress_buffer, unsigned char **output_ptr)
 551 {
 552   trajcoder_base_compress(encode_ints,3,large_index,compress_buffer);
 553   Ptngc_writemanybits(coder,compress_buffer,nbits,output_ptr);
 554 #ifdef SHOWIT
 555   fprintf(stderr,"nbits=%d (%g)\n",nbits,nbits/3.);
 556   nbits_sum+=nbits;
 557   nvalues_sum+=3;
 558   fprintf(stderr,"Large: %d %d %d\n",encode_ints[0],encode_ints[1],encode_ints[2]);
 559 #endif
 560 }
 561
 562 static void insert_batch(int *input_ptr, int ntriplets_left, int *prevcoord,int *minint, int *encode_ints, int startenc, int *nenc)
 563 {
 564   int nencode=startenc*3;
 565   int tmp_prevcoord[3];
 566
 567   tmp_prevcoord[0]=prevcoord[0];
 568   tmp_prevcoord[1]=prevcoord[1];
 569   tmp_prevcoord[2]=prevcoord[2];
 570
 571   if (startenc)
 572     {
 573       int i;
 574       for (i=0; i<startenc; i++)
 575         {
 576           tmp_prevcoord[0]+=encode_ints[i*3];
 577           tmp_prevcoord[1]+=encode_ints[i*3+1];
 578           tmp_prevcoord[2]+=encode_ints[i*3+2];
 579 #ifdef SHOWIT
 580           fprintf(stderr,"%6d: %6d %6d %6d\t\t%6d %6d %6d\t\t%6d %6d %6d\n",i*3,
 581                   tmp_prevcoord[0],tmp_prevcoord[1],tmp_prevcoord[2],
 582                   encode_ints[i*3],
 583                   encode_ints[i*3+1],
 584                   encode_ints[i*3+2],
 585                   positive_int(encode_ints[i*3]),
 586                   positive_int(encode_ints[i*3+1]),
 587                   positive_int(encode_ints[i*3+2]));
 588 #endif
 589         }
 590     }
 591
 592 #ifdef SHOWIT
 593   fprintf(stderr,"New batch\n");
 594 #endif
 595   while ((nencode<21) && (nencode<ntriplets_left*3))
 596     {
 597       encode_ints[nencode]=input_ptr[nencode]-minint[0]-tmp_prevcoord[0];
 598       encode_ints[nencode+1]=input_ptr[nencode+1]-minint[1]-tmp_prevcoord[1];
 599       encode_ints[nencode+2]=input_ptr[nencode+2]-minint[2]-tmp_prevcoord[2];
 600 #ifdef SHOWIT
 601       fprintf(stderr,"%6d: %6d %6d %6d\t\t%6d %6d %6d\t\t%6d %6d %6d\n",nencode,
 602               input_ptr[nencode]-minint[0],
 603               input_ptr[nencode+1]-minint[1],
 604               input_ptr[nencode+2]-minint[2],
 605               encode_ints[nencode],
 606               encode_ints[nencode+1],
 607               encode_ints[nencode+2],
 608               positive_int(encode_ints[nencode]),
 609               positive_int(encode_ints[nencode+1]),
 610               positive_int(encode_ints[nencode+2]));
 611 #endif
 612       tmp_prevcoord[0]=input_ptr[nencode]-minint[0];
 613       tmp_prevcoord[1]=input_ptr[nencode+1]-minint[1];
 614       tmp_prevcoord[2]=input_ptr[nencode+2]-minint[2];
 615       nencode+=3;
 616     }
 617   *nenc=nencode;
 618 }
 619
 620 static void flush_large(struct coder *coder, int *has_large, int *has_large_ints, int n,
 621                         int *large_index, int large_nbits, unsigned char *compress_buffer,
 622                         unsigned char **output_ptr)
 623 {
 624   int i;
 625   if (n<3)
 626     {
 627       for (i=0; i<n; i++)
 628         {
 629           write_instruction(coder,INSTR_ONLY_LARGE,output_ptr);
 630           write_three_large(coder,has_large_ints+i*3,large_index,large_nbits,compress_buffer,output_ptr);
 631         }
 632     }
 633   else
 634     {
 635 #if 0
 636       printf("CODING RLE: %d\n",n);
 637 #endif
 638       write_instruction(coder,INSTR_LARGE_RLE,output_ptr);
 639       Ptngc_writebits(coder,n-3,4,output_ptr); /* Number of large atoms in this sequence: 3 to 18 */
 640       for (i=0; i<n; i++)
 641         write_three_large(coder,has_large_ints+i*3,large_index,large_nbits,compress_buffer,output_ptr);
 642     }
 643   if (((*has_large)-n)!=0)
 644     {
 645       int j;
 646       for (i=0; i<(*has_large)-n; i++)
 647         for (j=0; j<3; j++)
 648           has_large_ints[i*3+j]=has_large_ints[(i+n)*3+j];
 649     }
 650   *has_large=(*has_large)-n; /* Number of remaining large atoms in buffer */
 651 }
 652
 653 static void buffer_large(struct coder *coder, int *has_large, int *has_large_ints, int *new_large_ints,
 654                         int *large_index, int large_nbits, unsigned char *compress_buffer,
 655                         unsigned char **output_ptr)
 656 {
 657   /* If it is full we must write them all. */
 658   if (*has_large==18)
 659     flush_large(coder,has_large,has_large_ints, *has_large,
 660                 large_index,large_nbits,compress_buffer,output_ptr); /* Flush them all. */
 661   has_large_ints[(*has_large)*3]=new_large_ints[0];
 662   has_large_ints[(*has_large)*3+1]=new_large_ints[1];
 663   has_large_ints[(*has_large)*3+2]=new_large_ints[2];
 664   *has_large=(*has_large)+1;
 665 }
 666
 667
 668 unsigned char *Ptngc_pack_array_xtc2(struct coder *coder,int *input, int *length)
 669 {
 670   unsigned char *output=NULL;
 671   unsigned char *output_ptr=NULL;
 672   int i,ienc,j;
 673   int output_length=0;
 674   /* Pack triplets. */
 675   int ntriplets=*length/3;
 676   int intmax;
 677   int max_small;
 678   int large_index[3];
 679   int max_large_index;
 680   int large_nbits;
 681   int small_index;
 682   int small_idx[3];
 683   int minint[3],maxint[3];
 684   int has_large=0;
 685   int has_large_ints[54]; /* Large cache. Up to 18 large atoms. */
 686   int prevcoord[3];
 687   int runlength=0; /* Initial runlength. "Stupidly" set to zero for
 688                       simplicity and explicity */
 689   int swapatoms=0; /* Initial guess is that we should not swap the
 690                       first two atoms in each large+small
 691                       transition */
 692   int didswap; /* Whether swapping was actually done. */
 693   int *input_ptr=input;
 694   int encode_ints[21]; /* Up to 3 large + 18 small ints can be encoded at once */
 695   int nencode;
 696   unsigned char compress_buffer[18*4]; /* Holds compressed result for 3 large ints or up to 18 small ints. */
 697   int ntriplets_left=ntriplets;
 698   int refused=0;
 699 #ifdef SHOWIT
 700   nbits_sum=0;
 701   nvalues_sum=0;
 702 #endif
 703   /* Allocate enough memory for output */
 704   output=warnmalloc(8* *length*sizeof *output);
 705   output_ptr=output;
 706
 707   maxint[0]=minint[0]=input[0];
 708   maxint[1]=minint[1]=input[1];
 709   maxint[2]=minint[2]=input[2];
 710
 711   for (i=1; i<ntriplets; i++)
 712     for (j=0; j<3; j++)
 713       {
 714         if (input[i*3+j]>maxint[j])
 715           maxint[j]=input[i*3+j];
 716         if (input[i*3+j]<minint[j])
 717           minint[j]=input[i*3+j];
 718       }
 719
 720   large_index[0]=Ptngc_find_magic_index(maxint[0]-minint[0]+1);
 721   large_index[1]=Ptngc_find_magic_index(maxint[1]-minint[1]+1);
 722   large_index[2]=Ptngc_find_magic_index(maxint[2]-minint[2]+1);
 723   large_nbits=compute_magic_bits(large_index);
 724   max_large_index=large_index[0];
 725   if (large_index[1]>max_large_index)
 726     max_large_index=large_index[1];
 727   if (large_index[2]>max_large_index)
 728     max_large_index=large_index[2];
 729
 730 #ifdef SHOWIT
 731   for (j=0; j<3; j++)
 732     fprintf(stderr,"minint[%d]=%d. maxint[%d]=%d large_index[%d]=%d value=%d\n",j,minint[j],j,maxint[j],
 733             j,large_index[j],magic[large_index[j]]);
 734   fprintf(stderr,"large_nbits=%d\n",large_nbits);
 735 #endif
 736
 737
 738   /* Guess initial small index */
 739   small_index=max_large_index/2;
 740
 741   /* Find the largest value that is not large. Not large is half index of
 742      large. */
 743   max_small=magic[small_index];
 744   intmax=0;
 745   for (i=0; i<*length; i++)
 746     {
 747       int item=input[i];
 748       int s=positive_int(item);
 749       if (s>intmax)
 750         if (s<max_small)
 751           intmax=s;
 752     }
 753   /* This value is not critical, since if I guess wrong, the code will
 754      just insert instructions to increase this value. */
 755   small_index=Ptngc_find_magic_index(intmax);
 756 #ifdef SHOWIT
 757   fprintf(stderr,"initial_small intmax=%d. small_index=%d value=%d\n",intmax,small_index,magic[small_index]);
 758 #endif
 759
 760   /* Store min integers */
 761   coder->pack_temporary_bits=32;
 762   coder->pack_temporary=positive_int(minint[0]);
 763   Ptngc_out8bits(coder,&output_ptr);
 764   coder->pack_temporary_bits=32;
 765   coder->pack_temporary=positive_int(minint[1]);
 766   Ptngc_out8bits(coder,&output_ptr);
 767   coder->pack_temporary_bits=32;
 768   coder->pack_temporary=positive_int(minint[2]);
 769   Ptngc_out8bits(coder,&output_ptr);
 770   /* Store max indices */
 771   coder->pack_temporary_bits=8;
 772   coder->pack_temporary=large_index[0];
 773   Ptngc_out8bits(coder,&output_ptr);
 774   coder->pack_temporary_bits=8;
 775   coder->pack_temporary=large_index[1];
 776   Ptngc_out8bits(coder,&output_ptr);
 777   coder->pack_temporary_bits=8;
 778   coder->pack_temporary=large_index[2];
 779   Ptngc_out8bits(coder,&output_ptr);
 780   /* Store initial small index */
 781   coder->pack_temporary_bits=8;
 782   coder->pack_temporary=small_index;
 783   Ptngc_out8bits(coder,&output_ptr);
 784
 785 #if 0
 786 #ifdef SHOWIT
 787   for (i=0; i<ntriplets_left; i++)
 788     fprintf(stderr,"VALUE:%d %6d %6d %6d\n",
 789             i,
 790             input_ptr[i*3],
 791             input_ptr[i*3+1],
 792             input_ptr[i*3+2]);
 793 #endif
 794 #endif
 795
 796   /* Initial prevcoord is the minimum integers. */
 797   prevcoord[0]=minint[0];
 798   prevcoord[1]=minint[1];
 799   prevcoord[2]=minint[2];
 800
 801   while (ntriplets_left)
 802     {
 803       if (ntriplets_left<0)
 804         {
 805           fprintf(stderr,"TRAJNG: BUG! ntriplets_left<0!\n");
 806           exit(EXIT_FAILURE);
 807         }
 808       /* If only less than three atoms left we just write them all as large integers. Here no swapping is done! */
 809       if (ntriplets_left<3)
 810         {
 811           for (ienc=0; ienc<ntriplets_left; ienc++)
 812             {
 813               int jenc;
 814               for (jenc=0; jenc<3; jenc++)
 815                 encode_ints[jenc]=input_ptr[ienc*3+jenc]-minint[jenc];
 816               buffer_large(coder,&has_large,has_large_ints,encode_ints,large_index,large_nbits,compress_buffer,&output_ptr);
 817               input_ptr+=3;
 818               ntriplets_left--;
 819             }
 820           flush_large(coder,&has_large,has_large_ints,has_large,large_index,large_nbits,compress_buffer,&output_ptr);
 821         }
 822       else
 823         {
 824           int min_runlength=0;
 825           int largest_required_base;
 826           int largest_runlength_base;
 827           int largest_required_index;
 828           int largest_runlength_index;
 829           int new_runlength;
 830           int new_small_index;
 831           int iter_runlength;
 832           int iter_small_index;
 833           int rle_index_dep;
 834           didswap=0;
 835           /* Insert the next batch of integers to be encoded into the buffer */
 836 #ifdef SHOWIT
 837           fprintf(stderr,"Initial batch\n");
 838 #endif
 839           insert_batch(input_ptr,ntriplets_left,prevcoord,minint,encode_ints,0,&nencode);
 840
 841           /* First we must decide if the next value is large (does not reasonably fit in current small encoding)
 842              Also, if we have not written any values yet, we must begin by writing a large atom. */
 843           if ((input_ptr==input) || (is_quite_large(encode_ints,small_index,max_large_index)) || (refused))
 844             {
 845               /* If any of the next two atoms are large we should probably write them as large and not swap them */
 846               int no_swap=0;
 847               if ((is_quite_large(encode_ints+3,small_index,max_large_index)) || (is_quite_large(encode_ints+6,small_index,max_large_index)))
 848                 no_swap=1;
 849               if (!no_swap)
 850                 {
 851                   /* Next we must decide if we should swap the first
 852                      two values. */
 853 #if 1
 854                   swapdecide(coder,input_ptr,&swapatoms,large_index,minint,&output_ptr);
 855 #else
 856                   swapatoms=0;
 857 #endif
 858                   /* If we should do the integer swapping manipulation we should do it now */
 859                   if (swapatoms)
 860                     {
 861                       didswap=1;
 862                       for (i=0; i<3; i++)
 863                         {
 864                           int in[3], out[3];
 865                           in[0]=input_ptr[i]-minint[i];
 866                           in[1]=input_ptr[3+i]-input_ptr[i]; /* minint[i]-minint[i] cancels out */
 867                           in[2]=input_ptr[6+i]-input_ptr[3+i]; /* minint[i]-minint[i] cancels out */
 868                           swap_ints(in,out);
 869                           encode_ints[i]=out[0];
 870                           encode_ints[3+i]=out[1];
 871                           encode_ints[6+i]=out[2];
 872                         }
 873                       /* We have swapped atoms, so the minimum run-length is 2 */
 874 #ifdef SHOWIT
 875                       fprintf(stderr,"Swap atoms results in:\n");
 876                       for (i=0; i<3; i++)
 877                         fprintf(stderr,"%d: %6d %6d %6d\t\t%6d %6d %6d\n",i*3,
 878                                 encode_ints[i*3],
 879                                 encode_ints[i*3+1],
 880                                 encode_ints[i*3+2],
 881                                 positive_int(encode_ints[i*3]),
 882                                 positive_int(encode_ints[i*3+1]),
 883                                 positive_int(encode_ints[i*3+2]));
 884
 885 #endif
 886                       min_runlength=2;
 887                     }
 888                 }
 889               if ((swapatoms) && (didswap))
 890                 {
 891                   for (ienc=0; ienc<3; ienc++)
 892                     prevcoord[ienc]=encode_ints[ienc];
 893                 }
 894               else
 895                 {
 896                   for (ienc=0; ienc<3; ienc++)
 897                     prevcoord[ienc]=input_ptr[ienc]-minint[ienc];
 898                 }
 899               /* Cache large value for later possible combination with
 900                  a sequence of small integers. */
 901               buffer_large(coder,&has_large,has_large_ints,prevcoord,large_index,large_nbits,compress_buffer,&output_ptr);
 902 #ifdef SHOWIT
 903               fprintf(stderr,"Prevcoord after packing of large: %d %d %d\n",
 904                       prevcoord[0],prevcoord[1],prevcoord[2]);
 905 #endif
 906
 907               /* We have written a large integer so we have one less atoms to worry about */
 908               input_ptr+=3;
 909               ntriplets_left--;
 910
 911               refused=0;
 912
 913               /* Insert the next batch of integers to be encoded into the buffer */
 914 #ifdef SHOWIT
 915               fprintf(stderr,"Update batch due to large int.\n");
 916 #endif
 917               if ((swapatoms) && (didswap))
 918                 {
 919                   /* Keep swapped values. */
 920                   for (i=0; i<2; i++)
 921                     for (ienc=0; ienc<3; ienc++)
 922                       encode_ints[i*3+ienc]=encode_ints[(i+1)*3+ienc];
 923                 }
 924               insert_batch(input_ptr,ntriplets_left,prevcoord,minint,encode_ints,min_runlength,&nencode);
 925             }
 926           /* Here we should only have differences for the atom coordinates. */
 927           /* Convert the ints to positive ints */
 928           for (ienc=0; ienc<nencode; ienc++)
 929             {
 930               int pint=positive_int(encode_ints[ienc]);
 931               encode_ints[ienc]=pint;
 932             }
 933           /* Now we must decide what base and runlength to do. If we have swapped atoms it will be at least 2.
 934              If even the next atom is large, we will not do anything. */
 935           largest_required_base=0;
 936           /* Determine required base */
 937           for (ienc=0; ienc<min_runlength*3; ienc++)
 938             if (encode_ints[ienc]>largest_required_base)
 939               largest_required_base=encode_ints[ienc];
 940           /* Also compute what the largest base is for the current runlength setting! */
 941           largest_runlength_base=0;
 942           for (ienc=0; (ienc<runlength*3) && (ienc<nencode); ienc++)
 943             if (encode_ints[ienc]>largest_runlength_base)
 944               largest_runlength_base=encode_ints[ienc];
 945
 946           largest_required_index=Ptngc_find_magic_index(largest_required_base);
 947           largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base);
 948
 949           if (largest_required_index<largest_runlength_index)
 950             {
 951               new_runlength=min_runlength;
 952               new_small_index=largest_required_index;
 953             }
 954           else
 955             {
 956               new_runlength=runlength;
 957               new_small_index=largest_runlength_index;
 958             }
 959
 960           /* Only allow increase of runlength wrt min_runlength */
 961           if (new_runlength<min_runlength)
 962             new_runlength=min_runlength;
 963
 964           /* If the current runlength is longer than the number of
 965              triplets left stop it from being so. */
 966           if (new_runlength>ntriplets_left)
 967             new_runlength=ntriplets_left;
 968
 969           /* We must at least try to get some small integers going. */
 970           if (new_runlength==0)
 971             {
 972               new_runlength=1;
 973               new_small_index=small_index;
 974             }
 975
 976           iter_runlength=new_runlength;
 977           iter_small_index=new_small_index;
 978
 979           /* Iterate to find optimal encoding and runlength */
 980 #ifdef SHOWIT
 981           fprintf(stderr,"Entering iterative loop.\n");
 982           fflush(stderr);
 983 #endif
 984
 985           do {
 986             new_runlength=iter_runlength;
 987             new_small_index=iter_small_index;
 988
 989 #ifdef SHOWIT
 990             fprintf(stderr,"Test new_small_index=%d Base=%d\n",new_small_index,magic[new_small_index]);
 991 #endif
 992             /* What is the largest runlength
 993                we can do with the currently
 994                selected encoding? Also the max supported runlength is 6 triplets! */
 995             for (ienc=0; ienc<nencode && ienc<18; ienc++)
 996               {
 997                 int test_index=Ptngc_find_magic_index(encode_ints[ienc]);
 998                 if (test_index>new_small_index)
 999                   break;
1000               }
1001             if (ienc/3>new_runlength)
1002               {
1003                 iter_runlength=ienc/3;
1004 #ifdef SHOWIT
1005                 fprintf(stderr,"I found a new possible runlength: %d\n",iter_runlength);
1006 #endif
1007               }
1008
1009             /* How large encoding do we have to use? */
1010             largest_runlength_base=0;
1011             for (ienc=0; ienc<iter_runlength*3; ienc++)
1012               if (encode_ints[ienc]>largest_runlength_base)
1013                 largest_runlength_base=encode_ints[ienc];
1014             largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base);
1015             if (largest_runlength_index!=new_small_index)
1016               {
1017                 iter_small_index=largest_runlength_index;
1018 #ifdef SHOWIT
1019                 fprintf(stderr,"I found a new possible small index: %d Base=%d\n",iter_small_index,magic[iter_small_index]);
1020 #endif
1021               }
1022           } while ((new_runlength!=iter_runlength) ||
1023                    (new_small_index!=iter_small_index));
1024
1025 #ifdef SHOWIT
1026           fprintf(stderr,"Exit iterative loop.\n");
1027           fflush(stderr);
1028 #endif
1029
1030           /* Verify that we got something good. We may have caught a
1031              substantially larger atom. If so we should just bail
1032              out and let the loop get on another lap. We may have a
1033              minimum runlength though and then we have to fulfill
1034              the request to write out these atoms! */
1035           rle_index_dep=0;
1036           if (new_runlength<3)
1037             rle_index_dep=IS_LARGE;
1038           else if (new_runlength<6)
1039             rle_index_dep=QUITE_LARGE;
1040           if ((min_runlength)
1041               || ((new_small_index<small_index+IS_LARGE) && (new_small_index+rle_index_dep<max_large_index))
1042 #if 1
1043               || (new_small_index+IS_LARGE<max_large_index)
1044 #endif
1045 )
1046             {
1047               int nbits;
1048               if ((new_runlength!=runlength) || (new_small_index!=small_index))
1049                 {
1050                   int change=new_small_index-small_index;
1051
1052                   if (new_small_index<=0)
1053                     change=0;
1054
1055                   if (change<0)
1056                     {
1057                       int ixx;
1058                       for (ixx=0; ixx<new_runlength; ixx++)
1059                         {
1060                           int rejected;
1061                           do {
1062                             int ixyz;
1063                             double isum=0.; /* ints can be almost 32 bit so multiplication will overflow. So do doubles. */
1064                             for (ixyz=0; ixyz<3; ixyz++)
1065                               {
1066                                 /* encode_ints is already positive (and multiplied by 2 versus the original, just as magic ints) */
1067                                 double id=encode_ints[ixx*3+ixyz];
1068                                 isum+=id*id;
1069                               }
1070                             rejected=0;
1071 #ifdef SHOWIT
1072                             fprintf(stderr,"Tested decrease %d of index: %g>=%g?\n",change,isum,(double)magic[small_index+change]*(double)magic[small_index+change]);
1073 #endif
1074                             if (isum>(double)magic[small_index+change]*(double)magic[small_index+change])
1075                               {
1076 #ifdef SHOWIT
1077                                 fprintf(stderr,"Rejected decrease %d of index due to length of vector: %g>=%g\n",change,isum,(double)magic[small_index+change]*(double)magic[small_index+change]);
1078 #endif
1079                                 rejected=1;
1080                                 change++;
1081                               }
1082                           } while ((change<0) && (rejected));
1083                           if (change==0)
1084                             break;
1085                         }
1086                     }
1087
1088                   /* If the only thing to do is to change the base by
1089                      only one -1 it is probably not worth it. */
1090                   if (!((change==-1) && (runlength==new_runlength)))
1091                     {
1092                       /* If we have a very short runlength we do not
1093                          want to do large base changes.  It costs 6
1094                          extra bits to do -2. We gain 2/3
1095                          bits per value to decrease the index by -2,
1096                          ie 2 bits, so to any changes down we must
1097                          have a runlength of 3 or more to do it for
1098                          one molecule! If we have several molecules we
1099                          will gain of course, so not be so strict. */
1100                       if ((change==-2) && (new_runlength<3))
1101                         {
1102                           if (runlength==new_runlength)
1103                             change=0;
1104                           else
1105                             change=-1;
1106 #ifdef SHOWIT
1107                           fprintf(stderr,"Rejected change by -2 due to too short runlenght. Change set to %d\n",change);
1108 #endif
1109                         }
1110
1111                       /* First adjust base using large base change instruction (if necessary) */
1112                       while ((change>1) || (change<-1) || ((new_runlength==6) && (change)))
1113                         {
1114                           unsigned int code=0U;
1115                           int this_change=change;
1116                           if (this_change>2)
1117                             this_change=2;
1118                           if (this_change<-2)
1119                             this_change=-2;
1120                           change-=this_change;
1121 #ifdef SHOWIT
1122                           fprintf(stderr,"Large base change: %d.\n",this_change);
1123 #endif
1124                           small_index+=this_change;
1125                           if (this_change<0)
1126                             {
1127                               code|=2U;
1128                               this_change=-this_change;
1129                             }
1130                           code|=(unsigned int)(this_change-1);
1131                           write_instruction(coder,INSTR_LARGE_BASE_CHANGE,&output_ptr);
1132                           Ptngc_writebits(coder,code,2,&output_ptr);
1133                         }
1134                       /* If there is still base change or runlength changes to do, we do them now. */
1135                       if ((new_runlength!=runlength) || (change))
1136                         {
1137                           unsigned int ichange=(unsigned int)(change+1);
1138                           unsigned int code=0U;
1139                           unsigned int irun=(unsigned int)((new_runlength-1)*3);
1140                           if (new_runlength==6)
1141                             ichange=0; /* Means no change. The change has been taken care of explicitly by a large
1142                                           base change instruction above. */
1143                           code=ichange+irun;
1144 #ifdef SHOWIT
1145                           fprintf(stderr,"Small base change: %d Runlength change: %d\n",change,new_runlength);
1146 #endif
1147                           small_index+=change;
1148                           write_instruction(coder,INSTR_BASE_RUNLENGTH,&output_ptr);
1149                           Ptngc_writebits(coder,code,4,&output_ptr);
1150                           runlength=new_runlength;
1151                         }
1152                     }
1153 #ifdef SHOWIT
1154                   else
1155                     fprintf(stderr,"Rejected base change due to only change==-1\n");
1156 #endif
1157 #ifdef SHOWIT
1158                   fprintf(stderr,"Current small index: %d Base=%d\n",small_index,magic[small_index]);
1159 #endif
1160                 }
1161               /* If we have a large previous integer we can combine it with a sequence of small ints. */
1162               if (has_large)
1163                 {
1164                   /* If swapatoms is set to 1 but we did actually not
1165                      do any swapping, we must first write out the
1166                      large atom and then the small. If swapatoms is 1
1167                      and we did swapping we can use the efficient
1168                      encoding. */
1169                   if ((swapatoms) && (!didswap))
1170                     {
1171 #ifdef SHOWIT
1172                       fprintf(stderr,"Swapatoms was set to 1 but we did not do swapping!\n");
1173                       fprintf(stderr,"Only one large integer.\n");
1174 #endif
1175                       /* Flush all large atoms. */
1176                       flush_large(coder,&has_large,has_large_ints,has_large,large_index,large_nbits,compress_buffer,&output_ptr);
1177 #ifdef SHOWIT
1178                       fprintf(stderr,"Sequence of only small integers.\n");
1179 #endif
1180                       write_instruction(coder,INSTR_ONLY_SMALL,&output_ptr);
1181                     }
1182                   else
1183                     {
1184
1185 #ifdef SHOWIT
1186                       fprintf(stderr,"Sequence of one large and small integers (good compression).\n");
1187 #endif
1188                       /* Flush all large atoms but one! */
1189                       if (has_large>1)
1190                         flush_large(coder,&has_large,has_large_ints,has_large-1,large_index,large_nbits,compress_buffer,&output_ptr);
1191                       write_instruction(coder,INSTR_DEFAULT,&output_ptr);
1192                       write_three_large(coder,has_large_ints,large_index,large_nbits,compress_buffer,&output_ptr);
1193                       has_large=0;
1194                     }
1195                 }
1196               else
1197                 {
1198 #ifdef SHOWIT
1199                   fprintf(stderr,"Sequence of only small integers.\n");
1200 #endif
1201                   write_instruction(coder,INSTR_ONLY_SMALL,&output_ptr);
1202                 }
1203               /* Base compress small integers using the current parameters. */
1204               nbits=magic_bits[small_index][runlength-1];
1205               /* The same base is used for the small changes. */
1206               small_idx[0]=small_index;
1207               small_idx[1]=small_index;
1208               small_idx[2]=small_index;
1209               trajcoder_base_compress(encode_ints,runlength*3,small_idx,compress_buffer);
1210 #ifdef SHOWIT
1211               fprintf(stderr,"nbits=%d (%g)\n",nbits,nbits/(runlength*3.));
1212               nbits_sum+=nbits;
1213               nvalues_sum+=runlength*3;
1214               fprintf(stderr,"Runlength encoded small integers. runlength=%d\n",runlength);
1215 #endif
1216               /* write out base compressed small integers */
1217               Ptngc_writemanybits(coder,compress_buffer,nbits,&output_ptr);
1218 #ifdef SHOWIT
1219               for (ienc=0; ienc<runlength; ienc++)
1220                 fprintf(stderr,"Small: %d %d %d\n",
1221                         encode_ints[ienc*3],
1222                         encode_ints[ienc*3+1],
1223                         encode_ints[ienc*3+2]);
1224 #endif
1225               /* Update prevcoord. */
1226               for (ienc=0; ienc<runlength; ienc++)
1227                 {
1228 #ifdef SHOWIT
1229                   fprintf(stderr,"Prevcoord in packing: %d %d %d\n",
1230                           prevcoord[0],prevcoord[1],prevcoord[2]);
1231 #endif
1232                   prevcoord[0]+=unpositive_int(encode_ints[ienc*3]);
1233                   prevcoord[1]+=unpositive_int(encode_ints[ienc*3+1]);
1234                   prevcoord[2]+=unpositive_int(encode_ints[ienc*3+2]);
1235                 }
1236 #ifdef SHOWIT
1237               fprintf(stderr,"Prevcoord in packing: %d %d %d\n",
1238                       prevcoord[0],prevcoord[1],prevcoord[2]);
1239 #endif
1240
1241               input_ptr+=3*runlength;
1242               ntriplets_left-=runlength;
1243             }
1244           else
1245             {
1246 #ifdef SHOWIT
1247               fprintf(stderr,"Refused value: %d old is %d max is %d\n",new_small_index,small_index,max_large_index);
1248               fflush(stderr);
1249 #endif
1250               refused=1;
1251             }
1252         }
1253 #ifdef SHOWIT
1254       fprintf(stderr,"Number of triplets left is %d\n",ntriplets_left);
1255 #endif
1256     }
1257   /* If we have large previous integers we must flush them now. */
1258   if (has_large)
1259     flush_large(coder,&has_large,has_large_ints,has_large,large_index,large_nbits,compress_buffer,&output_ptr);
1260   Ptngc_pack_flush(coder,&output_ptr);
1261   output_length=(int)(output_ptr-output);
1262 #ifdef SHOWIT
1263   fprintf(stderr,"Done block: nbits=%d nvalues=%d (%g)\n",nbits_sum,nvalues_sum,(double)nbits_sum/nvalues_sum);
1264 #endif
1265   *length=output_length;
1266   return output;
1267 }
1268
1269
1270 int Ptngc_unpack_array_xtc2(struct coder *coder,unsigned char *packed,int *output, int length)
1271 {
1272   unsigned char *ptr=packed;
1273   int bitptr=0;
1274   int minint[3];
1275   int large_index[3];
1276   int small_index;
1277   int prevcoord[3];
1278   int ntriplets_left=length/3;
1279   int swapatoms=0;
1280   int runlength=0;
1281   int large_nbits;
1282   unsigned char compress_buffer[18*4]; /* Holds compressed result for 3 large ints or up to 18 small ints. */
1283   int encode_ints[21]; /* Up to 3 large + 18 small ints can be encoded at once */
1284   (void)coder;
1285
1286   /* Read min integers. */
1287   minint[0]=unpositive_int(readbits(&ptr,&bitptr,32));
1288   minint[1]=unpositive_int(readbits(&ptr,&bitptr,32));
1289   minint[2]=unpositive_int(readbits(&ptr,&bitptr,32));
1290   /* Read large indices */
1291   large_index[0]=readbits(&ptr,&bitptr,8);
1292   large_index[1]=readbits(&ptr,&bitptr,8);
1293   large_index[2]=readbits(&ptr,&bitptr,8);
1294   /* Read small index */
1295   small_index=readbits(&ptr,&bitptr,8);
1296
1297   large_nbits=compute_magic_bits(large_index);
1298
1299 #ifdef SHOWIT
1300   fprintf(stderr,"Minimum integers: %d %d %d\n",minint[0],minint[1],minint[2]);
1301   fprintf(stderr,"Large indices: %d %d %d\n",large_index[0],large_index[1],large_index[2]);
1302   fprintf(stderr,"Small index: %d\n",small_index);
1303   fprintf(stderr,"large_nbits=%d\n",large_nbits);
1304 #endif
1305
1306   /* Initial prevcoord is the minimum integers. */
1307   prevcoord[0]=minint[0];
1308   prevcoord[1]=minint[1];
1309   prevcoord[2]=minint[2];
1310
1311   while (ntriplets_left)
1312     {
1313       int instr=read_instruction(&ptr,&bitptr);
1314 #ifdef SHOWIT
1315       if ((instr>=0) && (instr<MAXINSTR))
1316         fprintf(stderr,"Decoded instruction %s\n",instrnames[instr]);
1317 #endif
1318       if ((instr==INSTR_DEFAULT) /* large+small */
1319           || (instr==INSTR_ONLY_LARGE) /* only large */
1320           || (instr==INSTR_ONLY_SMALL)) /* only small */
1321         {
1322           int large_ints[3]={0,0,0};
1323           if (instr!=INSTR_ONLY_SMALL)
1324             {
1325               /* Clear the compress buffer. */
1326               int i;
1327               for (i=0; i<18*4; i++)
1328                 compress_buffer[i]=0;
1329               /* Get the large value. */
1330               readmanybits(&ptr,&bitptr,large_nbits,compress_buffer);
1331               trajcoder_base_decompress(compress_buffer,3,large_index,encode_ints);
1332               large_ints[0]=encode_ints[0];
1333               large_ints[1]=encode_ints[1];
1334               large_ints[2]=encode_ints[2];
1335 #ifdef SHOWIT
1336               fprintf(stderr,"large ints: %d %d %d\n",large_ints[0],large_ints[1],large_ints[2]);
1337 #endif
1338             }
1339           if (instr!=INSTR_ONLY_LARGE)
1340             {
1341               int small_idx[3];
1342               int i;
1343               /* The same base is used for the small changes. */
1344               small_idx[0]=small_index;
1345               small_idx[1]=small_index;
1346               small_idx[2]=small_index;
1347               /* Clear the compress buffer. */
1348               for (i=0; i<18*4; i++)
1349                 compress_buffer[i]=0;
1350               /* Get the small values. */
1351               readmanybits(&ptr,&bitptr,magic_bits[small_index][runlength-1],compress_buffer);
1352               trajcoder_base_decompress(compress_buffer,3*runlength,small_idx,encode_ints);
1353 #ifdef SHOWIT
1354               for (i=0; i<runlength; i++)
1355                 fprintf(stderr,"small ints: %d %d %d\n",encode_ints[i*3+0],encode_ints[i*3+1],encode_ints[i*3+2]);
1356 #endif
1357             }
1358           if (instr==INSTR_DEFAULT)
1359             {
1360               /* Check for swapped atoms */
1361               if (swapatoms)
1362                 {
1363                   /* Unswap the atoms. */
1364                   int i;
1365                   for (i=0; i<3; i++)
1366                     {
1367                       int in[3], out[3];
1368                       in[0]=large_ints[i];
1369                       in[1]=unpositive_int(encode_ints[i]);
1370                       in[2]=unpositive_int(encode_ints[3+i]);
1371                       swap_ints(in,out);
1372                       large_ints[i]=out[0];
1373                       encode_ints[i]=positive_int(out[1]);
1374                       encode_ints[3+i]=positive_int(out[2]);
1375                     }
1376                 }
1377             }
1378           /* Output result. */
1379           if (instr!=INSTR_ONLY_SMALL)
1380             {
1381               /* Output large value */
1382               *output++=large_ints[0]+minint[0];
1383               *output++=large_ints[1]+minint[1];
1384               *output++=large_ints[2]+minint[2];
1385               prevcoord[0]=large_ints[0];
1386               prevcoord[1]=large_ints[1];
1387               prevcoord[2]=large_ints[2];
1388 #ifdef SHOWIT
1389               fprintf(stderr,"Prevcoord after unpacking of large: %d %d %d\n",
1390                       prevcoord[0],prevcoord[1],prevcoord[2]);
1391               fprintf(stderr,"VALUE:%d %6d %6d %6d\n",
1392                       length/3-ntriplets_left,
1393                       prevcoord[0]+minint[0],
1394                       prevcoord[1]+minint[1],
1395                       prevcoord[2]+minint[2]);
1396 #endif
1397               ntriplets_left--;
1398             }
1399           if (instr!=INSTR_ONLY_LARGE)
1400             {
1401               /* Output small values */
1402               int i;
1403 #ifdef SHOWIT
1404               fprintf(stderr,"Prevcoord before unpacking of small: %d %d %d\n",
1405                           prevcoord[0],prevcoord[1],prevcoord[2]);
1406 #endif
1407               for (i=0; i<runlength; i++)
1408                 {
1409                   int v[3];
1410                   v[0]=unpositive_int(encode_ints[i*3]);
1411                   v[1]=unpositive_int(encode_ints[i*3+1]);
1412                   v[2]=unpositive_int(encode_ints[i*3+2]);
1413                   prevcoord[0]+=v[0];
1414                   prevcoord[1]+=v[1];
1415                   prevcoord[2]+=v[2];
1416 #ifdef SHOWIT
1417                   fprintf(stderr,"Prevcoord after unpacking of small: %d %d %d\n",
1418                           prevcoord[0],prevcoord[1],prevcoord[2]);
1419                   fprintf(stderr,"Unpacked small values: %6d %6d %6d\t\t%6d %6d %6d\n",v[0],v[1],v[2],prevcoord[0],prevcoord[1],prevcoord[2]);
1420                   fprintf(stderr,"VALUE:%d %6d %6d %6d\n",
1421                           length/3-(ntriplets_left-i),
1422                           prevcoord[0]+minint[0],
1423                           prevcoord[1]+minint[1],
1424                           prevcoord[2]+minint[2]);
1425 #endif
1426                   *output++=prevcoord[0]+minint[0];
1427                   *output++=prevcoord[1]+minint[1];
1428                   *output++=prevcoord[2]+minint[2];
1429                 }
1430               ntriplets_left-=runlength;
1431             }
1432         }
1433       else if (instr==INSTR_LARGE_RLE)
1434         {
1435           int i,j;
1436           int large_ints[3];
1437           /* How many large atoms in this sequence? */
1438           int n=(int)readbits(&ptr,&bitptr,4)+3; /* 3-18 large atoms */
1439           for (i=0; i<n; i++)
1440             {
1441               /* Clear the compress buffer. */
1442               for (j=0; j<18*4; j++)
1443                 compress_buffer[j]=0;
1444               /* Get the large value. */
1445               readmanybits(&ptr,&bitptr,large_nbits,compress_buffer);
1446               trajcoder_base_decompress(compress_buffer,3,large_index,encode_ints);
1447               large_ints[0]=encode_ints[0];
1448               large_ints[1]=encode_ints[1];
1449               large_ints[2]=encode_ints[2];
1450               /* Output large value */
1451               *output++=large_ints[0]+minint[0];
1452               *output++=large_ints[1]+minint[1];
1453               *output++=large_ints[2]+minint[2];
1454               prevcoord[0]=large_ints[0];
1455               prevcoord[1]=large_ints[1];
1456               prevcoord[2]=large_ints[2];
1457             }
1458           ntriplets_left-=n;
1459         }
1460       else if (instr==INSTR_BASE_RUNLENGTH)
1461         {
1462           unsigned int code=readbits(&ptr,&bitptr,4);
1463           int change;
1464           if (code==15)
1465             {
1466               change=0;
1467               runlength=6;
1468             }
1469           else
1470             {
1471               int ichange=code%3;
1472               runlength=code/3+1;
1473               change=ichange-1;
1474             }
1475           small_index+=change;
1476         }
1477       else if (instr==INSTR_FLIP)
1478         {
1479           swapatoms=1-swapatoms;
1480         }
1481       else if (instr==INSTR_LARGE_BASE_CHANGE)
1482         {
1483           unsigned int ichange=readbits(&ptr,&bitptr,2);
1484           int change=(int)(ichange&0x1U)+1;
1485           if (ichange&0x2U)
1486             change=-change;
1487           small_index+=change;
1488         }
1489       else
1490         {
1491           fprintf(stderr,"TRAJNG: BUG! Encoded unknown instruction.\n");
1492           exit(EXIT_FAILURE);
1493         }
1494 #ifdef SHOWIT
1495       fprintf(stderr,"Number of triplets left is %d\n",ntriplets_left);
1496 #endif
1497     }
1498   return 0;
1499 }