src/external/tng_io/src/compression/xtc2.c

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