src/external/tng_io/src/compression/xtc3.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 /* The cost estimates are ripped right out of xtc2.c, so take these
  19    with a grain (truckload) of salt. */
  20
  21 #include <limits.h>
  22 #include <stdio.h>
  23 #include <stdlib.h>
  24 #include <string.h>
  25 #include <math.h>
  26 #include "../../include/compression/warnmalloc.h"
  27 #include "../../include/compression/widemuldiv.h"
  28 #include "../../include/compression/bwlzh.h"
  29
  30 static const double iflipgaincheck=0.89089871814033927; /*  1./(2**(1./6)) */
  31
  32 #define MAX_LARGE_RLE 1024 /* Maximum number of large atoms for large RLE. */
  33 #define MAX_SMALL_RLE 12 /* Maximum number of small atoms in one group. */
  34
  35 #define TRESHOLD_INTRA_INTER_DIRECT 1.5 /* How much larger can the direct
  36                                            frame deltas for the small
  37                                            triplets be and be accepted anyway
  38                                            as better than the intra/inter frame
  39                                            deltas. For better instructions/RLEs. */
  40
  41 #define TRESHOLD_INTER_INTRA 5.0 /* How much larger can the intra
  42                                     frame deltas for the small
  43                                     triplets be and be accepted anyway
  44                                     as better than the inter frame
  45                                     deltas. */
  46
  47 /* Difference in indices used for determining whether to store as
  48    large or small. A fun detail in this compression algorithm is that
  49    if everything works fine, large can often be smaller than small, or
  50    at least not as large as is large in magic.c. This is a key idea of
  51    xtc3. */
  52 #define QUITE_LARGE 3
  53 #define IS_LARGE 6
  54
  55 #if 0
  56 #define SHOWIT
  57 #endif
  58
  59 #if 0
  60 #define SHOWIT_LIGHT
  61 #endif
  62
  63 /* These routines are in xtc2.c */
  64 int Ptngc_magic(unsigned int i);
  65 int Ptngc_find_magic_index(unsigned int maxval);
  66
  67 static unsigned int positive_int(int item)
  68 {
  69   int s=0;
  70   if (item>0)
  71     s=1+(item-1)*2;
  72   else if (item<0)
  73     s=2+(-item-1)*2;
  74   return s;
  75 }
  76
  77 static int unpositive_int(int val)
  78 {
  79   int s=(val+1)/2;
  80   if ((val%2)==0)
  81     s=-s;
  82   return s;
  83 }
  84
  85
  86 /* Sequence instructions */
  87 #define INSTR_DEFAULT 0U
  88 #define INSTR_SMALL_RUNLENGTH 1U
  89 #define INSTR_ONLY_LARGE 2U
  90 #define INSTR_ONLY_SMALL 3U
  91 #define INSTR_FLIP 4U
  92 #define INSTR_LARGE_RLE 5U
  93 #define INSTR_LARGE_DIRECT 6U
  94 #define INSTR_LARGE_INTRA_DELTA 7U
  95 #define INSTR_LARGE_INTER_DELTA 8U
  96
  97 #define MAXINSTR 9
  98
  99 struct xtc3_context
 100 {
 101   unsigned int *instructions;
 102   int ninstr, ninstr_alloc;
 103   unsigned int *rle;
 104   int nrle, nrle_alloc;
 105   unsigned int *large_direct;
 106   int nlargedir, nlargedir_alloc;
 107   unsigned int *large_intra_delta;
 108   int nlargeintra, nlargeintra_alloc;
 109   unsigned int *large_inter_delta;
 110   int nlargeinter, nlargeinter_alloc;
 111   unsigned int *smallintra;
 112   int nsmallintra, nsmallintra_alloc;
 113   int minint[3],maxint[3];
 114   int has_large;
 115   int has_large_ints[MAX_LARGE_RLE*3]; /* Large cache. */
 116   int has_large_type[MAX_LARGE_RLE]; /* What kind of type this large
 117                                         int is. */
 118   int current_large_type;
 119 };
 120
 121 static void init_xtc3_context(struct xtc3_context *xtc3_context)
 122 {
 123   xtc3_context->instructions=NULL;
 124   xtc3_context->ninstr=0;
 125   xtc3_context->ninstr_alloc=0;
 126   xtc3_context->rle=NULL;
 127   xtc3_context->nrle=0;
 128   xtc3_context->nrle_alloc=0;
 129   xtc3_context->large_direct=NULL;
 130   xtc3_context->nlargedir=0;
 131   xtc3_context->nlargedir_alloc=0;
 132   xtc3_context->large_intra_delta=NULL;
 133   xtc3_context->nlargeintra=0;
 134   xtc3_context->nlargeintra_alloc=0;
 135   xtc3_context->large_inter_delta=NULL;
 136   xtc3_context->nlargeinter=0;
 137   xtc3_context->nlargeinter_alloc=0;
 138   xtc3_context->smallintra=NULL;
 139   xtc3_context->nsmallintra=0;
 140   xtc3_context->nsmallintra_alloc=0;
 141   xtc3_context->has_large=0;
 142   xtc3_context->current_large_type=0;
 143 }
 144
 145 static void free_xtc3_context(struct xtc3_context *xtc3_context)
 146 {
 147   free(xtc3_context->instructions);
 148   free(xtc3_context->rle);
 149   free(xtc3_context->large_direct);
 150   free(xtc3_context->large_intra_delta);
 151   free(xtc3_context->large_inter_delta);
 152   free(xtc3_context->smallintra);
 153 }
 154
 155 /* Modifies three integer values for better compression of water */
 156 static void swap_ints(int *in, int *out)
 157 {
 158   out[0]=in[0]+in[1];
 159   out[1]=-in[1];
 160   out[2]=in[1]+in[2];
 161 }
 162
 163 static void swap_is_better(int *input, int *minint, int *sum_normal, int *sum_swapped)
 164 {
 165   int normal_max=0;
 166   int swapped_max=0;
 167   int i,j;
 168   int normal[3];
 169   int swapped[3];
 170   for (i=0; i<3; i++)
 171     {
 172       normal[0]=input[i]-minint[i];
 173       normal[1]=input[3+i]-input[i]; /* minint[i]-minint[i] cancels out */
 174       normal[2]=input[6+i]-input[3+i]; /* minint[i]-minint[i] cancels out */
 175       swap_ints(normal,swapped);
 176       for (j=1; j<3; j++)
 177         {
 178           if (positive_int(normal[j])>(unsigned int)normal_max)
 179             normal_max=positive_int(normal[j]);
 180           if (positive_int(swapped[j])>(unsigned int)swapped_max)
 181             swapped_max=positive_int(swapped[j]);
 182         }
 183     }
 184   if (normal_max==0)
 185     normal_max=1;
 186   if (swapped_max==0)
 187     swapped_max=1;
 188   *sum_normal=normal_max;
 189   *sum_swapped=swapped_max;
 190 }
 191
 192 static void allocate_enough_memory(unsigned int **ptr, int *nele, int *nele_alloc)
 193 {
 194   (*nele)++;
 195   if (*nele>*nele_alloc)
 196     {
 197       *nele_alloc=*nele + *nele/2;
 198       *ptr=warnrealloc(*ptr,*nele_alloc*sizeof **ptr);
 199     }
 200 }
 201
 202 static void insert_value_in_array(unsigned int **ptr, int *nele, int *nele_alloc,
 203                                   unsigned int value,
 204                                   char *arrayname)
 205 {
 206 #ifndef SHOWIT
 207   (void)arrayname;
 208 #endif
 209   allocate_enough_memory(ptr,nele,nele_alloc);
 210 #ifdef SHOWIT
 211   fprintf(stderr,"Inserting value %u into array %s @ %d\n",value,arrayname,(*nele)-1);
 212 #endif
 213   (*ptr)[(*nele)-1]=value;
 214 }
 215
 216
 217
 218 static void swapdecide(struct xtc3_context *xtc3_context, int *input,int *swapatoms, int *large_index, int *minint)
 219 {
 220   int didswap=0;
 221   int normal,swapped;
 222   (void)large_index;
 223   swap_is_better(input,minint,&normal,&swapped);
 224   /* We have to determine if it is worth to change the behaviour.
 225      If diff is positive it means that it is worth something to
 226      swap. But it costs 4 bits to do the change. If we assume that
 227      we gain 0.17 bit by the swap per value, and the runlength>2
 228      for four molecules in a row, we gain something. So check if we
 229      gain at least 0.17 bits to even attempt the swap.
 230   */
 231 #ifdef SHOWIT
 232   fprintf(stderr,"Trying Flip: %g %g\n",(double)swapped/normal, (double)normal/swapped);
 233 #endif
 234   if (((swapped<normal) && (fabs((double)swapped/normal)<iflipgaincheck)) ||
 235       ((normal<swapped) && (fabs((double)normal/swapped)<iflipgaincheck)))
 236     {
 237       if (swapped<normal)
 238         {
 239           if (!*swapatoms)
 240             {
 241               *swapatoms=1;
 242               didswap=1;
 243             }
 244         }
 245       else
 246         {
 247           if (*swapatoms)
 248             {
 249               *swapatoms=0;
 250               didswap=1;
 251             }
 252         }
 253     }
 254   if (didswap)
 255     {
 256 #ifdef SHOWIT
 257       fprintf(stderr,"Flip. Swapatoms is now %d\n",*swapatoms);
 258 #endif
 259       insert_value_in_array(&xtc3_context->instructions,
 260                             &xtc3_context->ninstr,
 261                             &xtc3_context->ninstr_alloc,
 262                             INSTR_FLIP,"instr");
 263     }
 264 }
 265
 266 /* It is "large" if we have to increase the small index quite a
 267    bit. Not so much to be rejected by the not very large check
 268    later. */
 269 static int is_quite_large(int *input, int small_index, int max_large_index)
 270 {
 271   int is=0;
 272   int i;
 273   if (small_index+QUITE_LARGE>=max_large_index)
 274     is=1;
 275   else
 276     {
 277       for (i=0; i<3; i++)
 278         if (positive_int(input[i])>(unsigned int)Ptngc_magic(small_index+QUITE_LARGE))
 279           {
 280             is=1;
 281             break;
 282           }
 283     }
 284   return is;
 285 }
 286
 287 #ifdef SHOWIT
 288 int nbits_sum;
 289 int nvalues_sum;
 290 #endif
 291
 292 static void insert_batch(int *input_ptr, int ntriplets_left, int *prevcoord, int *encode_ints, int startenc, int *nenc)
 293 {
 294   int nencode=startenc*3;
 295   int tmp_prevcoord[3];
 296
 297   tmp_prevcoord[0]=prevcoord[0];
 298   tmp_prevcoord[1]=prevcoord[1];
 299   tmp_prevcoord[2]=prevcoord[2];
 300
 301   if (startenc)
 302     {
 303       int i;
 304       for (i=0; i<startenc; i++)
 305         {
 306           tmp_prevcoord[0]+=encode_ints[i*3];
 307           tmp_prevcoord[1]+=encode_ints[i*3+1];
 308           tmp_prevcoord[2]+=encode_ints[i*3+2];
 309 #ifdef SHOWIT
 310           fprintf(stderr,"%6d: %6d %6d %6d\t\t%6d %6d %6d\t\t%6d %6d %6d\n",i*3,
 311                   tmp_prevcoord[0],tmp_prevcoord[1],tmp_prevcoord[2],
 312                   encode_ints[i*3],
 313                   encode_ints[i*3+1],
 314                   encode_ints[i*3+2],
 315                   positive_int(encode_ints[i*3]),
 316                   positive_int(encode_ints[i*3+1]),
 317                   positive_int(encode_ints[i*3+2]));
 318 #endif
 319         }
 320     }
 321
 322 #ifdef SHOWIT
 323   fprintf(stderr,"New batch\n");
 324 #endif
 325   while ((nencode<3+MAX_SMALL_RLE*3) && (nencode<ntriplets_left*3))
 326     {
 327       encode_ints[nencode]=input_ptr[nencode]-tmp_prevcoord[0];
 328       encode_ints[nencode+1]=input_ptr[nencode+1]-tmp_prevcoord[1];
 329       encode_ints[nencode+2]=input_ptr[nencode+2]-tmp_prevcoord[2];
 330 #ifdef SHOWIT
 331       fprintf(stderr,"%6d: %6d %6d %6d\t\t%6d %6d %6d\t\t%6d %6d %6d\n",nencode,
 332               input_ptr[nencode],
 333               input_ptr[nencode+1],
 334               input_ptr[nencode+2],
 335               encode_ints[nencode],
 336               encode_ints[nencode+1],
 337               encode_ints[nencode+2],
 338               positive_int(encode_ints[nencode]),
 339               positive_int(encode_ints[nencode+1]),
 340               positive_int(encode_ints[nencode+2]));
 341 #endif
 342       tmp_prevcoord[0]=input_ptr[nencode];
 343       tmp_prevcoord[1]=input_ptr[nencode+1];
 344       tmp_prevcoord[2]=input_ptr[nencode+2];
 345       nencode+=3;
 346     }
 347   *nenc=nencode;
 348 }
 349
 350 static void large_instruction_change(struct xtc3_context *xtc3_context, int i)
 351 {
 352   /* If the first large is of a different kind than the currently used we must
 353      emit an "instruction" to change the large type. */
 354   if (xtc3_context->has_large_type[i]!=xtc3_context->current_large_type)
 355     {
 356       unsigned int instr;
 357       xtc3_context->current_large_type=xtc3_context->has_large_type[i];
 358       if (xtc3_context->current_large_type==0)
 359         instr=INSTR_LARGE_DIRECT;
 360       else if (xtc3_context->current_large_type==1)
 361         instr=INSTR_LARGE_INTRA_DELTA;
 362       else
 363         instr=INSTR_LARGE_INTER_DELTA;
 364       insert_value_in_array(&xtc3_context->instructions,
 365                             &xtc3_context->ninstr,
 366                             &xtc3_context->ninstr_alloc,
 367                             instr,"instr");
 368     }
 369 }
 370
 371 static void write_three_large(struct xtc3_context *xtc3_context,
 372                               int i)
 373 {
 374   int m;
 375   if (xtc3_context->current_large_type==0)
 376     {
 377       for (m=0; m<3; m++)
 378         insert_value_in_array(&xtc3_context->large_direct,
 379                               &xtc3_context->nlargedir,
 380                               &xtc3_context->nlargedir_alloc,
 381                               xtc3_context->has_large_ints[i*3+m],"large direct");
 382     }
 383   else if (xtc3_context->current_large_type==1)
 384     {
 385       for (m=0; m<3; m++)
 386         insert_value_in_array(&xtc3_context->large_intra_delta,
 387                               &xtc3_context->nlargeintra,
 388                               &xtc3_context->nlargeintra_alloc,
 389                               xtc3_context->has_large_ints[i*3+m],"large intra");
 390     }
 391   else
 392     {
 393       for (m=0; m<3; m++)
 394         insert_value_in_array(&xtc3_context->large_inter_delta,
 395                               &xtc3_context->nlargeinter,
 396                               &xtc3_context->nlargeinter_alloc,
 397                               xtc3_context->has_large_ints[i*3+m],"large inter");
 398     }
 399 }
 400
 401 static void flush_large(struct xtc3_context *xtc3_context,
 402                         int n) /* How many to flush. */
 403 {
 404   int i;
 405   i=0;
 406   while (i<n)
 407     {
 408       int j,k;
 409       /* If the first large is of a different kind than the currently used we must
 410          emit an "instruction" to change the large type. */
 411       large_instruction_change(xtc3_context,i);
 412       /* How many large of the same kind in a row? */
 413       for (j=0;
 414            (i+j<n) &&
 415              (xtc3_context->has_large_type[i+j]==xtc3_context->has_large_type[i]);
 416            j++);
 417       if (j<3)
 418         {
 419           for (k=0; k<j; k++)
 420             {
 421               insert_value_in_array(&xtc3_context->instructions,
 422                                     &xtc3_context->ninstr,
 423                                     &xtc3_context->ninstr_alloc,
 424                                     INSTR_ONLY_LARGE,"instr");
 425               write_three_large(xtc3_context,i+k);
 426             }
 427         }
 428       else
 429         {
 430           insert_value_in_array(&xtc3_context->instructions,
 431                                 &xtc3_context->ninstr,
 432                                 &xtc3_context->ninstr_alloc,
 433                                 INSTR_LARGE_RLE,"instr");
 434           insert_value_in_array(&xtc3_context->rle,
 435                                 &xtc3_context->nrle,
 436                                 &xtc3_context->nrle_alloc,
 437                                 (unsigned int)j,"rle (large)");
 438           for (k=0; k<j; k++)
 439             write_three_large(xtc3_context,i+k);
 440         }
 441       i+=j;
 442     }
 443   if ((xtc3_context->has_large-n)!=0)
 444     {
 445       int j;
 446       for (i=0; i<xtc3_context->has_large-n; i++)
 447         {
 448           xtc3_context->has_large_type[i]=xtc3_context->has_large_type[i+n];
 449           for (j=0; j<3; j++)
 450             xtc3_context->has_large_ints[i*3+j]=xtc3_context->has_large_ints[(i+n)*3+j];
 451         }
 452     }
 453   xtc3_context->has_large-=n; /* Number of remaining large atoms in buffer */
 454 }
 455
 456 static double compute_intlen(unsigned int *ints)
 457 {
 458   /* The largest value. */
 459   unsigned int m=ints[0];
 460   if (ints[1]>m)
 461     m=ints[1];
 462   if (ints[2]>m)
 463     m=ints[2];
 464   return (double)m;
 465 }
 466
 467 static void buffer_large(struct xtc3_context *xtc3_context, int *input, int inpdata,
 468                          int natoms, int intradelta_ok)
 469 {
 470   unsigned int direct[3], intradelta[3]={0,}, interdelta[3]={0,};
 471   double minlen;
 472   int best_type;
 473   int frame=inpdata/(natoms*3);
 474   int atomframe=inpdata%(natoms*3);
 475   /* If it is full we must write them all. */
 476   if (xtc3_context->has_large==MAX_LARGE_RLE)
 477     flush_large(xtc3_context,xtc3_context->has_large); /* Flush all. */
 478   /* Find out which is the best choice for the large integer. Direct coding, or some
 479      kind of delta coding? */
 480   /* First create direct coding. */
 481   direct[0]=(unsigned int)(input[inpdata]-xtc3_context->minint[0]);
 482   direct[1]=(unsigned int)(input[inpdata+1]-xtc3_context->minint[1]);
 483   direct[2]=(unsigned int)(input[inpdata+2]-xtc3_context->minint[2]);
 484   minlen=compute_intlen(direct);
 485   best_type=0; /* Direct. */
 486 #if 1
 487   /* Then try intra coding if we can. */
 488   if ((intradelta_ok) && (atomframe>=3))
 489     {
 490       double thislen;
 491       intradelta[0]=positive_int(input[inpdata]-input[inpdata-3]);
 492       intradelta[1]=positive_int(input[inpdata+1]-input[inpdata-2]);
 493       intradelta[2]=positive_int(input[inpdata+2]-input[inpdata-1]);
 494       thislen=compute_intlen(intradelta);
 495       if (thislen*TRESHOLD_INTRA_INTER_DIRECT<minlen)
 496         {
 497           minlen=thislen;
 498           best_type=1; /* Intra delta */
 499         }
 500     }
 501 #endif
 502 #if 1
 503   /* Then try inter coding if we can. */
 504   if (frame>0)
 505     {
 506       double thislen;
 507       interdelta[0]=positive_int(input[inpdata]-input[inpdata-natoms*3]);
 508       interdelta[1]=positive_int(input[inpdata+1]-input[inpdata-natoms*3+1]);
 509       interdelta[2]=positive_int(input[inpdata+2]-input[inpdata-natoms*3+2]);
 510       thislen=compute_intlen(interdelta);
 511       if (thislen*TRESHOLD_INTRA_INTER_DIRECT<minlen)
 512         {
 513           best_type=2; /* Inter delta */
 514         }
 515     }
 516 #endif
 517   xtc3_context->has_large_type[xtc3_context->has_large]=best_type;
 518   if (best_type==0)
 519     {
 520       xtc3_context->has_large_ints[xtc3_context->has_large*3]=direct[0];
 521       xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=direct[1];
 522       xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=direct[2];
 523     }
 524   else if (best_type==1)
 525     {
 526       xtc3_context->has_large_ints[xtc3_context->has_large*3]=intradelta[0];
 527       xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=intradelta[1];
 528       xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=intradelta[2];
 529     }
 530   else if (best_type==2)
 531     {
 532       xtc3_context->has_large_ints[xtc3_context->has_large*3]=interdelta[0];
 533       xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=interdelta[1];
 534       xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=interdelta[2];
 535     }
 536   xtc3_context->has_large++;
 537 }
 538
 539 static void output_int(unsigned char *output,int *outdata, unsigned int n)
 540 {
 541   output[(*outdata)++]=((unsigned int)n)&0xFFU;
 542   output[(*outdata)++]=(((unsigned int)n)>>8)&0xFFU;
 543   output[(*outdata)++]=(((unsigned int)n)>>16)&0xFFU;
 544   output[(*outdata)++]=(((unsigned int)n)>>24)&0xFFU;
 545 }
 546
 547 #if 0
 548 static void printarray(unsigned int *a, int n, char *name)
 549 {
 550   int i;
 551   for (i=0; i<n; i++)
 552     fprintf(stderr,"%u %s\n",a[i],name);
 553 }
 554 #endif
 555
 556 /* The base_compress routine first compresses all x coordinates, then
 557    y and finally z. The bases used for each can be different. The
 558    MAXBASEVALS value determines how many coordinates are compressed
 559    into a single number. Only resulting whole bytes are dealt with for
 560    simplicity. MAXMAXBASEVALS is the insanely large value to accept
 561    files written with that value. BASEINTERVAL determines how often a
 562    new base is actually computed and stored in the output
 563    file. MAXBASEVALS*BASEINTERVAL values are stored using the same
 564    base in BASEINTERVAL different integers. Note that the primarily
 565    the decompression using a large MAXBASEVALS becomes very slow. */
 566 #define MAXMAXBASEVALS 16384U
 567 #define MAXBASEVALS 24U
 568 #define BASEINTERVAL 8
 569
 570 /* How many bytes are needed to store n values in base base */
 571 static int base_bytes(unsigned int base, int n)
 572 {
 573   int i,j;
 574   unsigned int largeint[MAXMAXBASEVALS+1];
 575   unsigned int largeint_tmp[MAXMAXBASEVALS+1];
 576   int numbytes=0;
 577   for (i=0; i<n+1; i++)
 578     largeint[i]=0U;
 579   for (i=0; i<n; i++)
 580     {
 581       if (i!=0)
 582         {
 583           Ptngc_largeint_mul(base,largeint,largeint_tmp,n+1);
 584           for (j=0; j<n+1; j++)
 585             largeint[j]=largeint_tmp[j];
 586         }
 587       Ptngc_largeint_add(base-1U,largeint,n+1);
 588     }
 589   for (i=0; i<n; i++)
 590     if (largeint[i])
 591       for (j=0; j<4; j++)
 592         if ((largeint[i]>>(j*8))&0xFFU)
 593           numbytes=i*4+j+1;
 594   return numbytes;
 595 }
 596
 597 static void base_compress(unsigned int *data, int len, unsigned char *output, int *outlen)
 598 {
 599   unsigned int largeint[MAXBASEVALS+1];
 600   unsigned int largeint_tmp[MAXBASEVALS+1];
 601   int ixyz, i;
 602   unsigned int j;
 603   int nwrittenout=0;
 604   unsigned int numbytes=0;
 605   /* Store the MAXBASEVALS value in the output. */
 606   output[nwrittenout++]=(unsigned char)(MAXBASEVALS&0xFFU);
 607   output[nwrittenout++]=(unsigned char)((MAXBASEVALS>>8)&0xFFU);
 608   /* Store the BASEINTERVAL value in the output. */
 609   output[nwrittenout++]=(unsigned char)(BASEINTERVAL&0xFFU);
 610   for (ixyz=0; ixyz<3; ixyz++)
 611     {
 612       unsigned int base=0U;
 613       int nvals=0;
 614       int basegiven=0;
 615       for (j=0; j<MAXBASEVALS+1; j++)
 616         largeint[j]=0U;
 617       for (i=ixyz; i<len; i+=3)
 618         {
 619          if (nvals==0)
 620            {
 621              int basecheckvals=0;
 622              int k;
 623              if (basegiven==0)
 624                {
 625                  base=0U;
 626                  /* Find the largest value for this particular coordinate. */
 627                  for (k=i; k<len; k+=3)
 628                    {
 629                      if (data[k]>base)
 630                        base=data[k];
 631                      basecheckvals++;
 632                      if (basecheckvals==MAXBASEVALS*BASEINTERVAL)
 633                        break;
 634                    }
 635                  /* The base is one larger than the largest values. */
 636                  base++;
 637                  if (base<2)
 638                    base=2;
 639                  /* Store the base in the output. */
 640                  output[nwrittenout++]=(unsigned char)(base&0xFFU);
 641                  output[nwrittenout++]=(unsigned char)((base>>8)&0xFFU);
 642                  output[nwrittenout++]=(unsigned char)((base>>16)&0xFFU);
 643                  output[nwrittenout++]=(unsigned char)((base>>24)&0xFFU);
 644                  basegiven=BASEINTERVAL;
 645                  /* How many bytes is needed to store MAXBASEVALS values using this base? */
 646                  numbytes=base_bytes(base,MAXBASEVALS);
 647                }
 648              basegiven--;
 649 #ifdef SHOWIT
 650              fprintf(stderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,MAXBASEVALS);
 651 #endif
 652            }
 653           if (nvals!=0)
 654             {
 655               Ptngc_largeint_mul(base,largeint,largeint_tmp,MAXBASEVALS+1);
 656               for (j=0; j<MAXBASEVALS+1; j++)
 657                 largeint[j]=largeint_tmp[j];
 658             }
 659           Ptngc_largeint_add(data[i],largeint,MAXBASEVALS+1);
 660 #ifdef SHOWIT
 661           fprintf(stderr,"outputting value %u\n",data[i]);
 662 #endif
 663           nvals++;
 664           if (nvals==MAXBASEVALS)
 665             {
 666 #ifdef SHOWIT
 667               fprintf(stderr,"Writing largeint: ");
 668 #endif
 669               for (j=0; j<numbytes; j++)
 670                 {
 671                   int ilarge=j/4;
 672                   int ibyte=j%4;
 673                   output[nwrittenout++]=(unsigned char)((largeint[ilarge]>>(ibyte*8))&(0xFFU));
 674 #ifdef SHOWIT
 675                   fprintf(stderr,"%02x",(unsigned int)output[nwrittenout-1]);
 676 #endif
 677                 }
 678 #ifdef SHOWIT
 679               fprintf(stderr,"\n");
 680 #endif
 681               nvals=0;
 682               for (j=0; j<MAXBASEVALS+1; j++)
 683                 largeint[j]=0U;
 684             }
 685         }
 686       if (nvals)
 687         {
 688           numbytes=base_bytes(base,nvals);
 689 #ifdef SHOWIT
 690           fprintf(stderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,nvals);
 691 #endif
 692           for (j=0; j<numbytes; j++)
 693             {
 694               int ilarge=j/4;
 695               int ibyte=j%4;
 696               output[nwrittenout++]=(unsigned char)((largeint[ilarge]>>(ibyte*8))&(0xFFU));
 697             }
 698         }
 699     }
 700   *outlen=nwrittenout;
 701 }
 702
 703 static void base_decompress(unsigned char *input, int len, unsigned int *output)
 704 {
 705   unsigned int largeint[MAXMAXBASEVALS+1];
 706   unsigned int largeint_tmp[MAXMAXBASEVALS+1];
 707   int ixyz, i, j;
 708   int maxbasevals=(int)((unsigned int)(input[0])|(((unsigned int)(input[1]))<<8));
 709   int baseinterval=(int)input[2];
 710   if (maxbasevals>(int)MAXMAXBASEVALS)
 711     {
 712       fprintf(stderr,"Read a larger maxbasevals value from the file than I can handle. Fix"
 713               " by increasing MAXMAXBASEVALS to at least %d. Although, this is"
 714               " probably a bug in TRAJNG, since MAXMAXBASEVALS should already be insanely large enough.\n",maxbasevals);
 715       exit(EXIT_FAILURE);
 716     }
 717   input+=3;
 718   for (ixyz=0; ixyz<3; ixyz++)
 719     {
 720       int numbytes=0;
 721       int nvals_left=len/3;
 722       int outvals=ixyz;
 723       int basegiven=0;
 724       unsigned int base=0U;
 725 #ifdef SHOWIT
 726       fprintf(stderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,maxbasevals);
 727 #endif
 728       while (nvals_left)
 729         {
 730           int n;
 731           if (basegiven==0)
 732             {
 733               base=(unsigned int)(input[0])|
 734                 (((unsigned int)(input[1]))<<8)|
 735                 (((unsigned int)(input[2]))<<16)|
 736                 (((unsigned int)(input[3]))<<24);
 737               input+=4;
 738               basegiven=baseinterval;
 739               /* How many bytes is needed to store maxbasevals values using this base? */
 740               numbytes=base_bytes(base,maxbasevals);
 741             }
 742           basegiven--;
 743           if (nvals_left<maxbasevals)
 744             {
 745               numbytes=base_bytes(base,nvals_left);
 746 #ifdef SHOWIT
 747               fprintf(stderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,nvals_left);
 748 #endif
 749             }
 750           for (j=0; j<maxbasevals+1; j++)
 751             largeint[j]=0U;
 752 #ifdef SHOWIT
 753           fprintf(stderr,"Reading largeint: ");
 754 #endif
 755           if (numbytes/4 < maxbasevals+1)
 756             {
 757               for (j=0; j<numbytes; j++)
 758                 {
 759                   int ilarge=j/4;
 760                   int ibyte=j%4;
 761                   largeint[ilarge]|=((unsigned int)input[j])<<(ibyte*8);
 762 #ifdef SHOWIT
 763                   fprintf(stderr,"%02x",(unsigned int)input[j]);
 764 #endif
 765                 }
 766             }
 767 #ifdef SHOWIT
 768           fprintf(stderr,"\n");
 769 #endif
 770           input+=numbytes;
 771           /* Do the long division required to get the output values. */
 772           n=maxbasevals;
 773           if (n>nvals_left)
 774             n=nvals_left;
 775           for (i=n-1; i>=0; i--)
 776             {
 777               output[outvals+i*3]=Ptngc_largeint_div(base,largeint,largeint_tmp,maxbasevals+1);
 778               for (j=0; j<maxbasevals+1; j++)
 779                 largeint[j]=largeint_tmp[j];
 780             }
 781 #ifdef SHOWIT
 782           for (i=0; i<n; i++)
 783             fprintf(stderr,"outputting value %u\n",output[outvals+i*3]);
 784 #endif
 785           outvals+=n*3;
 786           nvals_left-=n;
 787         }
 788     }
 789 }
 790
 791 /* If a large proportion of the integers are large (More than 10\% are >14 bits) we return 0, otherwise 1 */
 792 static int heuristic_bwlzh(unsigned int *ints, int nints)
 793 {
 794   int i,num;
 795   num=0;
 796   for (i=0; i<nints; i++)
 797     if (ints[i]>=16384)
 798       num++;
 799   if (num>nints/10)
 800     return 0;
 801   else
 802     return 1;
 803 }
 804
 805 /* Speed selects how careful to try to find the most efficient compression. The BWLZH algo is expensive!
 806    Speed <=2 always avoids BWLZH everywhere it is possible.
 807    Speed 3 and 4 and 5 use heuristics (check proportion of large value). This should mostly be safe.
 808    Speed 5 enables the LZ77 component of BWLZH.
 809    Speed 6 always tests if BWLZH is better and if it is uses it. This can be very slow.
 810  */
 811 unsigned char *Ptngc_pack_array_xtc3(int *input, int *length, int natoms, int speed)
 812 {
 813   unsigned char *output=NULL;
 814   int i,ienc,j;
 815   int outdata=0;
 816   /* Pack triplets. */
 817   int ntriplets=*length/3;
 818   int intmax;
 819   int max_small;
 820   int small_index;
 821   int max_large_index;
 822   int large_index[3];
 823   int prevcoord[3];
 824   int runlength=0; /* Initial runlength. "Stupidly" set to zero for
 825                       simplicity and explicity */
 826   int swapatoms=0; /* Initial guess is that we should not swap the
 827                       first two atoms in each large+small
 828                       transition */
 829   int didswap; /* Whether swapping was actually done. */
 830   int inpdata=0;
 831   int encode_ints[3+MAX_SMALL_RLE*3]; /* Up to 3 large + 24 small ints can be encoded at once */
 832   int nencode;
 833   int ntriplets_left=ntriplets;
 834   int refused=0;
 835   unsigned char *bwlzh_buf=NULL;
 836   int bwlzh_buf_len;
 837   unsigned char *base_buf=NULL;
 838   int base_buf_len;
 839
 840   struct xtc3_context xtc3_context;
 841   init_xtc3_context(&xtc3_context);
 842
 843   xtc3_context.maxint[0]=xtc3_context.minint[0]=input[0];
 844   xtc3_context.maxint[1]=xtc3_context.minint[1]=input[1];
 845   xtc3_context.maxint[2]=xtc3_context.minint[2]=input[2];
 846
 847   /* Values of speed should be sane. */
 848   if (speed<1)
 849     speed=1;
 850   if (speed>6)
 851     speed=6;
 852
 853 #ifdef SHOWIT
 854   nbits_sum=0;
 855   nvalues_sum=0;
 856 #endif
 857   /* Allocate enough memory for output */
 858   if (*length < 48)
 859     output=warnmalloc(8*48*sizeof *output);
 860   else
 861     output=warnmalloc(8* *length*sizeof *output);
 862
 863
 864   for (i=1; i<ntriplets; i++)
 865     for (j=0; j<3; j++)
 866       {
 867         if (input[i*3+j]>xtc3_context.maxint[j])
 868           xtc3_context.maxint[j]=input[i*3+j];
 869         if (input[i*3+j]<xtc3_context.minint[j])
 870           xtc3_context.minint[j]=input[i*3+j];
 871       }
 872
 873   large_index[0]=Ptngc_find_magic_index(xtc3_context.maxint[0]-xtc3_context.minint[0]+1);
 874   large_index[1]=Ptngc_find_magic_index(xtc3_context.maxint[1]-xtc3_context.minint[1]+1);
 875   large_index[2]=Ptngc_find_magic_index(xtc3_context.maxint[2]-xtc3_context.minint[2]+1);
 876   max_large_index=large_index[0];
 877   if (large_index[1]>max_large_index)
 878     max_large_index=large_index[1];
 879   if (large_index[2]>max_large_index)
 880     max_large_index=large_index[2];
 881
 882 #ifdef SHOWIT
 883   for (j=0; j<3; j++)
 884     fprintf(stderr,"minint[%d]=%d. maxint[%d]=%d large_index[%d]=%d value=%d\n",j,xtc3_context.minint[j],j,xtc3_context.maxint[j],
 885             j,large_index[j],Ptngc_magic(large_index[j]));
 886 #endif
 887
 888   /* Guess initial small index */
 889   small_index=max_large_index/2;
 890
 891   /* Find the largest value that is not large. Not large is half index of
 892      large. */
 893   max_small=Ptngc_magic(small_index);
 894   intmax=0;
 895   for (i=0; i<*length; i++)
 896     {
 897       int item=input[i];
 898       int s=positive_int(item);
 899       if (s>intmax)
 900         if (s<max_small)
 901           intmax=s;
 902     }
 903   /* This value is not critical, since if I guess wrong, the code will
 904      just insert instructions to increase this value. */
 905   small_index=Ptngc_find_magic_index(intmax);
 906 #ifdef SHOWIT
 907   fprintf(stderr,"initial small_index=%d value=%d\n",small_index,Ptngc_magic(small_index));
 908 #endif
 909
 910   output_int(output,&outdata,positive_int(xtc3_context.minint[0]));
 911   output_int(output,&outdata,positive_int(xtc3_context.minint[1]));
 912   output_int(output,&outdata,positive_int(xtc3_context.minint[2]));
 913
 914 #if 0
 915 #ifdef SHOWIT
 916   for (i=0; i<ntriplets_left; i++)
 917     fprintf(stderr,"VALUE:%d %6d %6d %6d\n",
 918             i,
 919             input[inpdata+i*3],
 920             input[inpdata+i*3+1],
 921             input[inpdata+i*3+2]);
 922 #endif
 923 #endif
 924
 925   /* Initial prevcoord is the minimum integers. */
 926   prevcoord[0]=xtc3_context.minint[0];
 927   prevcoord[1]=xtc3_context.minint[1];
 928   prevcoord[2]=xtc3_context.minint[2];
 929
 930   while (ntriplets_left)
 931     {
 932       if (ntriplets_left<0)
 933         {
 934           fprintf(stderr,"TRAJNG: BUG! ntriplets_left<0!\n");
 935           exit(EXIT_FAILURE);
 936         }
 937       /* If only less than three atoms left we just write them all as large integers. Here no swapping is done! */
 938       if (ntriplets_left<3)
 939         {
 940           for (ienc=0; ienc<ntriplets_left; ienc++)
 941             {
 942               buffer_large(&xtc3_context,input,inpdata,natoms,1);
 943               inpdata+=3;
 944               ntriplets_left--;
 945             }
 946           flush_large(&xtc3_context,xtc3_context.has_large); /* Flush all */
 947         }
 948       else
 949         {
 950           int min_runlength=0;
 951           int largest_required_base;
 952           int largest_runlength_base;
 953           int largest_required_index;
 954           int largest_runlength_index;
 955           int new_runlength;
 956           int new_small_index;
 957           int iter_runlength;
 958           int iter_small_index;
 959           int rle_index_dep;
 960           didswap=0;
 961           /* Insert the next batch of integers to be encoded into the buffer */
 962 #ifdef SHOWIT
 963           fprintf(stderr,"Initial batch\n");
 964 #endif
 965           insert_batch(input+inpdata,ntriplets_left,prevcoord,encode_ints,0,&nencode);
 966
 967           /* First we must decide if the next value is large (does not reasonably fit in current small encoding)
 968              Also, if we have not written any values yet, we must begin by writing a large atom. */
 969           if ((inpdata==0) || (is_quite_large(encode_ints,small_index,max_large_index)) || (refused))
 970             {
 971               /* If any of the next two atoms are large we should probably write them as large and not swap them */
 972               int no_swap=0;
 973               if ((is_quite_large(encode_ints+3,small_index,max_large_index)) || (is_quite_large(encode_ints+6,small_index,max_large_index)))
 974                 no_swap=1;
 975 #if 1
 976               if (!no_swap)
 977                 {
 978                   /* If doing inter-frame coding results in smaller values we should not do any swapping either. */
 979                   int frame=inpdata/(natoms*3);
 980                   if (frame>0)
 981                     {
 982                       unsigned int delta[3], delta2[3];
 983                       delta[0]=positive_int(input[inpdata+3]-input[inpdata-natoms*3+3]);
 984                       delta[1]=positive_int(input[inpdata+4]-input[inpdata-natoms*3+4]);
 985                       delta[2]=positive_int(input[inpdata+5]-input[inpdata-natoms*3+5]);
 986                       delta2[0]=positive_int(encode_ints[3]);
 987                       delta2[1]=positive_int(encode_ints[4]);
 988                       delta2[2]=positive_int(encode_ints[5]);
 989 #ifdef SHOWIT
 990                       fprintf(stderr,"A1: inter delta: %u %u %u. intra delta=%u %u %u\n",
 991                               delta[0],delta[1],delta[2],
 992                               delta2[0],delta2[1],delta2[2]);
 993 #endif
 994                       if (compute_intlen(delta)*TRESHOLD_INTER_INTRA<compute_intlen(delta2))
 995                         {
 996                           delta[0]=positive_int(input[inpdata+6]-input[inpdata-natoms*3+6]);
 997                           delta[1]=positive_int(input[inpdata+7]-input[inpdata-natoms*3+7]);
 998                           delta[2]=positive_int(input[inpdata+8]-input[inpdata-natoms*3+8]);
 999                           delta2[0]=positive_int(encode_ints[6]);
1000                           delta2[1]=positive_int(encode_ints[7]);
1001                           delta2[2]=positive_int(encode_ints[8]);
1002 #ifdef SHOWIT
1003                           fprintf(stderr,"A2: inter delta: %u %u %u. intra delta=%u %u %u\n",
1004                                   delta[0],delta[1],delta[2],
1005                                   delta2[0],delta2[1],delta2[2]);
1006 #endif
1007                           if (compute_intlen(delta)*TRESHOLD_INTER_INTRA<compute_intlen(delta2))
1008                             {
1009                               no_swap=1;
1010 #ifdef SHOWIT
1011                               fprintf(stderr,"SETTING NO SWAP!\n");
1012 #endif
1013                             }
1014                         }
1015                     }
1016                 }
1017 #endif
1018               if (!no_swap)
1019                 {
1020                   /* Next we must decide if we should swap the first
1021                      two values. */
1022 #if 1
1023                   swapdecide(&xtc3_context,input+inpdata,&swapatoms,large_index,xtc3_context.minint);
1024 #else
1025                   swapatoms=0;
1026 #endif
1027                   /* If we should do the integer swapping manipulation we should do it now */
1028                   if (swapatoms)
1029                     {
1030                       didswap=1;
1031                       for (i=0; i<3; i++)
1032                         {
1033                           int in[3], out[3];
1034                           in[0]=input[inpdata+i];
1035                           in[1]=input[inpdata+3+i]-input[inpdata+i];
1036                           in[2]=input[inpdata+6+i]-input[inpdata+3+i];
1037                           swap_ints(in,out);
1038                           encode_ints[i]=out[0];
1039                           encode_ints[3+i]=out[1];
1040                           encode_ints[6+i]=out[2];
1041                         }
1042                       /* We have swapped atoms, so the minimum run-length is 2 */
1043 #ifdef SHOWIT
1044                       fprintf(stderr,"Swap atoms results in:\n");
1045                       for (i=0; i<3; i++)
1046                         fprintf(stderr,"%d: %6d %6d %6d\t\t%6d %6d %6d\n",i*3,
1047                                 encode_ints[i*3],
1048                                 encode_ints[i*3+1],
1049                                 encode_ints[i*3+2],
1050                                 positive_int(encode_ints[i*3]),
1051                                 positive_int(encode_ints[i*3+1]),
1052                                 positive_int(encode_ints[i*3+2]));
1053
1054 #endif
1055                       min_runlength=2;
1056                     }
1057                 }
1058               /* Cache large value for later possible combination with
1059                  a sequence of small integers. */
1060               if ((swapatoms) && (didswap))
1061                 {
1062                   buffer_large(&xtc3_context,input,inpdata+3,natoms,0); /* This is a swapped integer, so inpdata is one atom later and
1063                                                                            intra coding is not ok. */
1064                   for (ienc=0; ienc<3; ienc++)
1065                     prevcoord[ienc]=input[inpdata+3+ienc];
1066                 }
1067               else
1068                 {
1069                   buffer_large(&xtc3_context,input,inpdata,natoms,1);
1070                   for (ienc=0; ienc<3; ienc++)
1071                     prevcoord[ienc]=input[inpdata+ienc];
1072                 }
1073
1074
1075 #ifdef SHOWIT
1076               fprintf(stderr,"Prevcoord after packing of large: %d %d %d\n",
1077                       prevcoord[0],prevcoord[1],prevcoord[2]);
1078 #endif
1079
1080               /* We have written a large integer so we have one less atoms to worry about */
1081               inpdata+=3;
1082               ntriplets_left--;
1083
1084               refused=0;
1085
1086               /* Insert the next batch of integers to be encoded into the buffer */
1087 #ifdef SHOWIT
1088               fprintf(stderr,"Update batch due to large int.\n");
1089 #endif
1090               if ((swapatoms) && (didswap))
1091                 {
1092                   /* Keep swapped values. */
1093                   for (i=0; i<2; i++)
1094                     for (ienc=0; ienc<3; ienc++)
1095                       encode_ints[i*3+ienc]=encode_ints[(i+1)*3+ienc];
1096                 }
1097               insert_batch(input+inpdata,ntriplets_left,prevcoord,encode_ints,min_runlength,&nencode);
1098             }
1099           /* Here we should only have differences for the atom coordinates. */
1100           /* Convert the ints to positive ints */
1101           for (ienc=0; ienc<nencode; ienc++)
1102             {
1103               int pint=positive_int(encode_ints[ienc]);
1104               encode_ints[ienc]=pint;
1105             }
1106           /* Now we must decide what base and runlength to do. If we have swapped atoms it will be at least 2.
1107              If even the next atom is large, we will not do anything. */
1108           largest_required_base=0;
1109           /* Determine required base */
1110           for (ienc=0; ienc<min_runlength*3; ienc++)
1111             if (encode_ints[ienc]>largest_required_base)
1112               largest_required_base=encode_ints[ienc];
1113           /* Also compute what the largest base is for the current runlength setting! */
1114           largest_runlength_base=0;
1115           for (ienc=0; (ienc<runlength*3) && (ienc<nencode); ienc++)
1116             if (encode_ints[ienc]>largest_runlength_base)
1117               largest_runlength_base=encode_ints[ienc];
1118
1119           largest_required_index=Ptngc_find_magic_index(largest_required_base);
1120           largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base);
1121
1122           if (largest_required_index<largest_runlength_index)
1123             {
1124               new_runlength=min_runlength;
1125               new_small_index=largest_required_index;
1126             }
1127           else
1128             {
1129               new_runlength=runlength;
1130               new_small_index=largest_runlength_index;
1131             }
1132
1133           /* Only allow increase of runlength wrt min_runlength */
1134           if (new_runlength<min_runlength)
1135             new_runlength=min_runlength;
1136
1137           /* If the current runlength is longer than the number of
1138              triplets left stop it from being so. */
1139           if (new_runlength>ntriplets_left)
1140             new_runlength=ntriplets_left;
1141
1142           /* We must at least try to get some small integers going. */
1143           if (new_runlength==0)
1144             {
1145               new_runlength=1;
1146               new_small_index=small_index;
1147             }
1148
1149           iter_runlength=new_runlength;
1150           iter_small_index=new_small_index;
1151
1152           /* Iterate to find optimal encoding and runlength */
1153 #ifdef SHOWIT
1154           fprintf(stderr,"Entering iterative loop.\n");
1155           fflush(stderr);
1156 #endif
1157
1158           do {
1159             new_runlength=iter_runlength;
1160             new_small_index=iter_small_index;
1161
1162 #ifdef SHOWIT
1163             fprintf(stderr,"Test new_small_index=%d Base=%d\n",new_small_index,Ptngc_magic(new_small_index));
1164 #endif
1165             /* What is the largest runlength
1166                we can do with the currently
1167                selected encoding? Also the max supported runlength is MAX_SMALL_RLE triplets! */
1168             for (ienc=0; ienc<nencode && ienc<MAX_SMALL_RLE*3; ienc++)
1169               {
1170                 int test_index=Ptngc_find_magic_index(encode_ints[ienc]);
1171                 if (test_index>new_small_index)
1172                   break;
1173               }
1174             if (ienc/3>new_runlength)
1175               {
1176                 iter_runlength=ienc/3;
1177 #ifdef SHOWIT
1178                 fprintf(stderr,"I found a new possible runlength: %d\n",iter_runlength);
1179 #endif
1180               }
1181
1182             /* How large encoding do we have to use? */
1183             largest_runlength_base=0;
1184             for (ienc=0; ienc<iter_runlength*3; ienc++)
1185               if (encode_ints[ienc]>largest_runlength_base)
1186                 largest_runlength_base=encode_ints[ienc];
1187             largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base);
1188             if (largest_runlength_index!=new_small_index)
1189               {
1190                 iter_small_index=largest_runlength_index;
1191 #ifdef SHOWIT
1192                 fprintf(stderr,"I found a new possible small index: %d Base=%d\n",iter_small_index,Ptngc_magic(iter_small_index));
1193 #endif
1194               }
1195           } while ((new_runlength!=iter_runlength) ||
1196                    (new_small_index!=iter_small_index));
1197
1198 #ifdef SHOWIT
1199           fprintf(stderr,"Exit iterative loop.\n");
1200           fflush(stderr);
1201 #endif
1202
1203           /* Verify that we got something good. We may have caught a
1204              substantially larger atom. If so we should just bail
1205              out and let the loop get on another lap. We may have a
1206              minimum runlength though and then we have to fulfill
1207              the request to write out these atoms! */
1208           rle_index_dep=0;
1209           if (new_runlength<3)
1210             rle_index_dep=IS_LARGE;
1211           else if (new_runlength<6)
1212             rle_index_dep=QUITE_LARGE;
1213           if ((min_runlength)
1214               || ((new_small_index<small_index+IS_LARGE) && (new_small_index+rle_index_dep<max_large_index))
1215 #if 1
1216               || (new_small_index+IS_LARGE<max_large_index)
1217 #endif
1218 )
1219             {
1220               /* If doing inter-frame coding of large integers results
1221                  in smaller values than the small value we should not
1222                  produce a sequence of small values here. */
1223               int frame=inpdata/(natoms*3);
1224               int numsmaller=0;
1225 #if 1
1226               if ((!swapatoms) && (frame>0))
1227                 {
1228                   for (i=0; i<new_runlength; i++)
1229                     {
1230                       unsigned int delta[3];
1231                       unsigned int delta2[3];
1232                       delta[0]=positive_int(input[inpdata+i*3]-input[inpdata-natoms*3+i*3]);
1233                       delta[1]=positive_int(input[inpdata+i*3+1]-input[inpdata-natoms*3+i*3+1]);
1234                       delta[2]=positive_int(input[inpdata+i*3+2]-input[inpdata-natoms*3+i*3+2]);
1235                       delta2[0]=positive_int(encode_ints[i*3]);
1236                       delta2[1]=positive_int(encode_ints[i*3+1]);
1237                       delta2[2]=positive_int(encode_ints[i*3+2]);
1238                       if (compute_intlen(delta)*TRESHOLD_INTER_INTRA<compute_intlen(delta2))
1239                         numsmaller++;
1240                     }
1241                 }
1242 #endif
1243               /* Most of the values should become smaller, otherwise
1244                  we should encode them with intra coding. */
1245               if ((!swapatoms) && (numsmaller>=2*new_runlength/3))
1246                 {
1247                   /* Put all the values in large arrays, instead of the small array */
1248                   if (new_runlength)
1249                     {
1250                       for (i=0; i<new_runlength; i++)
1251                         buffer_large(&xtc3_context,input,inpdata+i*3,natoms,1);
1252                       for (i=0; i<3; i++)
1253                         prevcoord[i]=input[inpdata+(new_runlength-1)*3+i];
1254                       inpdata+=3*new_runlength;
1255                       ntriplets_left-=new_runlength;
1256                     }
1257                 }
1258               else
1259                 {
1260                   if ((new_runlength!=runlength) || (new_small_index!=small_index))
1261                     {
1262                       int change=new_small_index-small_index;
1263
1264                       if (new_small_index<=0)
1265                         change=0;
1266
1267                       if (change<0)
1268                         {
1269                           int ixx;
1270                           for (ixx=0; ixx<new_runlength; ixx++)
1271                             {
1272                               int rejected;
1273                               do {
1274                                 int ixyz;
1275                                 double isum=0.; /* ints can be almost 32 bit so multiplication will overflow. So do doubles. */
1276                                 for (ixyz=0; ixyz<3; ixyz++)
1277                                   {
1278                                     /* encode_ints is already positive (and multiplied by 2 versus the original, just as magic ints) */
1279                                     double id=encode_ints[ixx*3+ixyz];
1280                                     isum+=id*id;
1281                                   }
1282                                 rejected=0;
1283 #ifdef SHOWIT
1284                                 fprintf(stderr,"Tested decrease %d of index: %g>=%g?\n",change,isum,(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change));
1285 #endif
1286                                 if (isum>(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change))
1287                                   {
1288 #ifdef SHOWIT
1289                                     fprintf(stderr,"Rejected decrease %d of index due to length of vector: %g>=%g\n",change,isum,(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change));
1290 #endif
1291                                     rejected=1;
1292                                     change++;
1293                                   }
1294                               } while ((change<0) && (rejected));
1295                               if (change==0)
1296                                 break;
1297                             }
1298                         }
1299
1300                       /* Always accept the new small indices here. */
1301                       small_index=new_small_index;
1302                       /* If we have a new runlength emit it */
1303                       if (runlength!=new_runlength)
1304                         {
1305                           runlength=new_runlength;
1306                           insert_value_in_array(&xtc3_context.instructions,
1307                                                 &xtc3_context.ninstr,
1308                                                 &xtc3_context.ninstr_alloc,
1309                                                 INSTR_SMALL_RUNLENGTH,"instr");
1310                           insert_value_in_array(&xtc3_context.rle,
1311                                                 &xtc3_context.nrle,
1312                                                 &xtc3_context.nrle_alloc,
1313                                                 (unsigned int)runlength,"rle (small)");
1314                         }
1315 #ifdef SHOWIT
1316                       fprintf(stderr,"Current small index: %d Base=%d\n",small_index,Ptngc_magic(small_index));
1317 #endif
1318                     }
1319                   /* If we have a large previous integer we can combine it with a sequence of small ints. */
1320                   if (xtc3_context.has_large)
1321                     {
1322                       /* If swapatoms is set to 1 but we did actually not
1323                          do any swapping, we must first write out the
1324                          large atom and then the small. If swapatoms is 1
1325                          and we did swapping we can use the efficient
1326                          encoding. */
1327                       if ((swapatoms) && (!didswap))
1328                         {
1329 #ifdef SHOWIT
1330                           fprintf(stderr,"Swapatoms was set to 1 but we did not do swapping!\n");
1331                           fprintf(stderr,"Only one large integer.\n");
1332 #endif
1333                           /* Flush all large atoms. */
1334                           flush_large(&xtc3_context,xtc3_context.has_large);
1335 #ifdef SHOWIT
1336                           fprintf(stderr,"Sequence of only small integers.\n");
1337 #endif
1338                           insert_value_in_array(&xtc3_context.instructions,
1339                                                 &xtc3_context.ninstr,
1340                                                 &xtc3_context.ninstr_alloc,
1341                                                 INSTR_ONLY_SMALL,"instr");
1342                         }
1343                       else
1344                         {
1345
1346 #ifdef SHOWIT
1347                           fprintf(stderr,"Sequence of one large and small integers (good compression).\n");
1348 #endif
1349                           /* Flush all large atoms but one! */
1350                           if (xtc3_context.has_large>1)
1351                             flush_large(&xtc3_context,xtc3_context.has_large-1);
1352
1353                           /* Here we must check if we should emit a large
1354                              type change instruction. */
1355                           large_instruction_change(&xtc3_context,0);
1356
1357                           insert_value_in_array(&xtc3_context.instructions,
1358                                                 &xtc3_context.ninstr,
1359                                                 &xtc3_context.ninstr_alloc,
1360                                                 INSTR_DEFAULT,"instr");
1361
1362                           write_three_large(&xtc3_context,0);
1363                           xtc3_context.has_large=0;
1364                         }
1365                     }
1366                   else
1367                     {
1368 #ifdef SHOWIT
1369                       fprintf(stderr,"Sequence of only small integers.\n");
1370 #endif
1371                       insert_value_in_array(&xtc3_context.instructions,
1372                                             &xtc3_context.ninstr,
1373                                             &xtc3_context.ninstr_alloc,
1374                                             INSTR_ONLY_SMALL,"instr");
1375                     }
1376                   /* Insert the small integers into the small integer array. */
1377                   for (ienc=0; ienc<runlength*3; ienc++)
1378                     insert_value_in_array(&xtc3_context.smallintra,
1379                                           &xtc3_context.nsmallintra,
1380                                           &xtc3_context.nsmallintra_alloc,
1381                                           (unsigned int)encode_ints[ienc],"smallintra");
1382
1383 #ifdef SHOWIT
1384                   for (ienc=0; ienc<runlength; ienc++)
1385                     fprintf(stderr,"Small: %d %d %d\n",
1386                             encode_ints[ienc*3],
1387                             encode_ints[ienc*3+1],
1388                             encode_ints[ienc*3+2]);
1389 #endif
1390                   /* Update prevcoord. */
1391                   for (ienc=0; ienc<runlength; ienc++)
1392                     {
1393 #ifdef SHOWIT
1394                       fprintf(stderr,"Prevcoord in packing: %d %d %d\n",
1395                               prevcoord[0],prevcoord[1],prevcoord[2]);
1396 #endif
1397                       prevcoord[0]+=unpositive_int(encode_ints[ienc*3]);
1398                       prevcoord[1]+=unpositive_int(encode_ints[ienc*3+1]);
1399                       prevcoord[2]+=unpositive_int(encode_ints[ienc*3+2]);
1400                     }
1401 #ifdef SHOWIT
1402                   fprintf(stderr,"Prevcoord in packing: %d %d %d\n",
1403                           prevcoord[0],prevcoord[1],prevcoord[2]);
1404 #endif
1405
1406                   inpdata+=3*runlength;
1407                   ntriplets_left-=runlength;
1408 #if 1
1409                 }
1410 #endif
1411             }
1412           else
1413             {
1414 #ifdef SHOWIT
1415               fprintf(stderr,"Refused value: %d old is %d max is %d\n",new_small_index,small_index,max_large_index);
1416               fflush(stderr);
1417 #endif
1418               refused=1;
1419             }
1420         }
1421 #ifdef SHOWIT
1422       fprintf(stderr,"Number of triplets left is %d\n",ntriplets_left);
1423 #endif
1424     }
1425
1426   /* If we have large previous integers we must flush them now. */
1427   if (xtc3_context.has_large)
1428     flush_large(&xtc3_context,xtc3_context.has_large);
1429
1430   /* Now it is time to compress all the data in the buffers with the bwlzh or base algo. */
1431
1432 #if 0
1433   /* Inspect the data. */
1434   printarray(xtc3_context.instructions,xtc3_context.ninstr,"A instr");
1435   printarray(xtc3_context.rle,xtc3_context.nrle,"A rle");
1436   printarray(xtc3_context.large_direct,xtc3_context.nlargedir,"A largedir");
1437   printarray(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,"A largeintra");
1438   printarray(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,"A largeinter");
1439   printarray(xtc3_context.smallintra,xtc3_context.nsmallintra,"A smallintra");
1440   exit(0);
1441 #endif
1442
1443 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1444   fprintf(stderr,"instructions: %d\n",xtc3_context.ninstr);
1445 #endif
1446
1447 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1448 #define bwlzh_compress bwlzh_compress_verbose
1449 #define bwlzh_compress_no_lz77 bwlzh_compress_no_lz77_verbose
1450 #endif
1451
1452   output_int(output,&outdata,(unsigned int)xtc3_context.ninstr);
1453   if (xtc3_context.ninstr)
1454     {
1455       bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.ninstr));
1456       if (speed>=5)
1457         bwlzh_compress(xtc3_context.instructions,xtc3_context.ninstr,bwlzh_buf,&bwlzh_buf_len);
1458       else
1459         bwlzh_compress_no_lz77(xtc3_context.instructions,xtc3_context.ninstr,bwlzh_buf,&bwlzh_buf_len);
1460       output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1461       memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1462       outdata+=bwlzh_buf_len;
1463       free(bwlzh_buf);
1464     }
1465
1466 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1467   fprintf(stderr,"rle: %d\n",xtc3_context.nrle);
1468 #endif
1469
1470   output_int(output,&outdata,(unsigned int)xtc3_context.nrle);
1471   if (xtc3_context.nrle)
1472     {
1473       bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nrle));
1474       if (speed>=5)
1475         bwlzh_compress(xtc3_context.rle,xtc3_context.nrle,bwlzh_buf,&bwlzh_buf_len);
1476       else
1477         bwlzh_compress_no_lz77(xtc3_context.rle,xtc3_context.nrle,bwlzh_buf,&bwlzh_buf_len);
1478       output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1479       memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1480       outdata+=bwlzh_buf_len;
1481       free(bwlzh_buf);
1482     }
1483
1484 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1485   fprintf(stderr,"large direct: %d\n",xtc3_context.nlargedir);
1486 #endif
1487
1488   output_int(output,&outdata,(unsigned int)xtc3_context.nlargedir);
1489   if (xtc3_context.nlargedir)
1490     {
1491       if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.large_direct,xtc3_context.nlargedir))))
1492         {
1493           bwlzh_buf=NULL;
1494           bwlzh_buf_len=INT_MAX;
1495         }
1496       else
1497         {
1498           bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nlargedir));
1499           if (speed>=5)
1500             bwlzh_compress(xtc3_context.large_direct,xtc3_context.nlargedir,bwlzh_buf,&bwlzh_buf_len);
1501           else
1502             bwlzh_compress_no_lz77(xtc3_context.large_direct,xtc3_context.nlargedir,bwlzh_buf,&bwlzh_buf_len);
1503         }
1504       /* If this can be written smaller using base compression we should do that. */
1505       base_buf=warnmalloc((xtc3_context.nlargedir+3)*sizeof(int));
1506       base_compress(xtc3_context.large_direct,xtc3_context.nlargedir,base_buf,&base_buf_len);
1507 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1508       fprintf(stderr,"Large direct: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len);
1509 #endif
1510       if (base_buf_len<bwlzh_buf_len)
1511         {
1512           output[outdata++]=0U;
1513           output_int(output,&outdata,(unsigned int)base_buf_len);
1514           memcpy(output+outdata,base_buf,base_buf_len);
1515           outdata+=base_buf_len;
1516         }
1517       else
1518         {
1519           output[outdata++]=1U;
1520           output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1521           memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1522           outdata+=bwlzh_buf_len;
1523         }
1524       free(bwlzh_buf);
1525       free(base_buf);
1526     }
1527
1528 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1529   fprintf(stderr,"large intra: %d\n",xtc3_context.nlargeintra);
1530 #endif
1531
1532   output_int(output,&outdata,(unsigned int)xtc3_context.nlargeintra);
1533   if (xtc3_context.nlargeintra)
1534     {
1535       if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.large_intra_delta,xtc3_context.nlargeintra))))
1536         {
1537           bwlzh_buf=NULL;
1538           bwlzh_buf_len=INT_MAX;
1539         }
1540       else
1541         {
1542           bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nlargeintra));
1543           if (speed>=5)
1544             bwlzh_compress(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,bwlzh_buf,&bwlzh_buf_len);
1545           else
1546             bwlzh_compress_no_lz77(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,bwlzh_buf,&bwlzh_buf_len);
1547         }
1548       /* If this can be written smaller using base compression we should do that. */
1549       base_buf=warnmalloc((xtc3_context.nlargeintra+3)*sizeof(int));
1550       base_compress(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,base_buf,&base_buf_len);
1551 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1552       fprintf(stderr,"Large intra: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len);
1553 #endif
1554       if (base_buf_len<bwlzh_buf_len)
1555         {
1556           output[outdata++]=0U;
1557           output_int(output,&outdata,(unsigned int)base_buf_len);
1558           memcpy(output+outdata,base_buf,base_buf_len);
1559           outdata+=base_buf_len;
1560         }
1561       else
1562         {
1563           output[outdata++]=1U;
1564           output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1565           memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1566           outdata+=bwlzh_buf_len;
1567         }
1568       free(bwlzh_buf);
1569       free(base_buf);
1570     }
1571
1572 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1573   fprintf(stderr,"large inter: %d\n",xtc3_context.nlargeinter);
1574 #endif
1575
1576   output_int(output,&outdata,(unsigned int)xtc3_context.nlargeinter);
1577   if (xtc3_context.nlargeinter)
1578     {
1579       if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.large_inter_delta,xtc3_context.nlargeinter))))
1580         {
1581           bwlzh_buf=NULL;
1582           bwlzh_buf_len=INT_MAX;
1583         }
1584       else
1585         {
1586           bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nlargeinter));
1587           if (speed>=5)
1588             bwlzh_compress(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,bwlzh_buf,&bwlzh_buf_len);
1589           else
1590             bwlzh_compress_no_lz77(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,bwlzh_buf,&bwlzh_buf_len);
1591         }
1592       /* If this can be written smaller using base compression we should do that. */
1593       base_buf=warnmalloc((xtc3_context.nlargeinter+3)*sizeof(int));
1594       base_compress(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,base_buf,&base_buf_len);
1595 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1596       fprintf(stderr,"Large inter: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len);
1597 #endif
1598       if (base_buf_len<bwlzh_buf_len)
1599         {
1600           output[outdata++]=0U;
1601           output_int(output,&outdata,(unsigned int)base_buf_len);
1602           memcpy(output+outdata,base_buf,base_buf_len);
1603           outdata+=base_buf_len;
1604         }
1605       else
1606         {
1607           output[outdata++]=1U;
1608           output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1609           memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1610           outdata+=bwlzh_buf_len;
1611         }
1612       free(bwlzh_buf);
1613       free(base_buf);
1614     }
1615
1616 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1617   fprintf(stderr,"small intra: %d\n",xtc3_context.nsmallintra);
1618 #endif
1619
1620   output_int(output,&outdata,(unsigned int)xtc3_context.nsmallintra);
1621   if (xtc3_context.nsmallintra)
1622     {
1623       if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.smallintra,xtc3_context.nsmallintra))))
1624         {
1625           bwlzh_buf=NULL;
1626           bwlzh_buf_len=INT_MAX;
1627         }
1628       else
1629         {
1630           bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nsmallintra));
1631           if (speed>=5)
1632             bwlzh_compress(xtc3_context.smallintra,xtc3_context.nsmallintra,bwlzh_buf,&bwlzh_buf_len);
1633           else
1634             bwlzh_compress_no_lz77(xtc3_context.smallintra,xtc3_context.nsmallintra,bwlzh_buf,&bwlzh_buf_len);
1635         }
1636       /* If this can be written smaller using base compression we should do that. */
1637       base_buf=warnmalloc((xtc3_context.nsmallintra+3)*sizeof(int));
1638       base_compress(xtc3_context.smallintra,xtc3_context.nsmallintra,base_buf,&base_buf_len);
1639 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1640       fprintf(stderr,"Small intra: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len);
1641 #endif
1642       if (base_buf_len<bwlzh_buf_len)
1643         {
1644           output[outdata++]=0U;
1645           output_int(output,&outdata,(unsigned int)base_buf_len);
1646           memcpy(output+outdata,base_buf,base_buf_len);
1647           outdata+=base_buf_len;
1648         }
1649       else
1650         {
1651           output[outdata++]=1U;
1652           output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1653           memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1654           outdata+=bwlzh_buf_len;
1655         }
1656       free(bwlzh_buf);
1657       free(base_buf);
1658     }
1659   *length=outdata;
1660
1661   free_xtc3_context(&xtc3_context);
1662   return output;
1663 }
1664
1665 static void decompress_bwlzh_block(unsigned char **ptr,
1666                                    int nvals,
1667                                    unsigned int **vals)
1668 {
1669   int bwlzh_buf_len=(int)(((unsigned int)(*ptr)[0]) |
1670                           (((unsigned int)(*ptr)[1])<<8) |
1671                           (((unsigned int)(*ptr)[2])<<16) |
1672                           (((unsigned int)(*ptr)[3])<<24));
1673   (*ptr)+=4;
1674   *vals=warnmalloc(nvals*sizeof (**vals));
1675   bwlzh_decompress(*ptr,nvals,*vals);
1676   (*ptr)+=bwlzh_buf_len;
1677 }
1678
1679 static void decompress_base_block(unsigned char **ptr,
1680                                   int nvals,
1681                                   unsigned int **vals)
1682 {
1683   int base_buf_len=(int)(((unsigned int)(*ptr)[0]) |
1684                          (((unsigned int)(*ptr)[1])<<8) |
1685                          (((unsigned int)(*ptr)[2])<<16) |
1686                          (((unsigned int)(*ptr)[3])<<24));
1687   (*ptr)+=4;
1688   *vals=warnmalloc(nvals*sizeof (**vals));
1689   base_decompress(*ptr,nvals,*vals);
1690   (*ptr)+=base_buf_len;
1691 }
1692
1693 static void unpack_one_large(struct xtc3_context *xtc3_context,
1694                              int *ilargedir, int *ilargeintra,
1695                              int *ilargeinter, int *prevcoord,
1696                              int *minint, int *output,
1697                              int outdata, int didswap,
1698                              int natoms, int current_large_type)
1699 {
1700   int large_ints[3]={0,0,0};
1701   if (current_large_type==0 && xtc3_context->large_direct)
1702     {
1703       large_ints[0]=(int)xtc3_context->large_direct[(*ilargedir)]+minint[0];
1704       large_ints[1]=(int)xtc3_context->large_direct[(*ilargedir)+1]+minint[1];
1705       large_ints[2]=(int)xtc3_context->large_direct[(*ilargedir)+2]+minint[2];
1706       (*ilargedir)+=3;
1707     }
1708   else if (current_large_type==1 && xtc3_context->large_intra_delta)
1709     {
1710       large_ints[0]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)])+prevcoord[0];
1711       large_ints[1]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)+1])+prevcoord[1];
1712       large_ints[2]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)+2])+prevcoord[2];
1713       (*ilargeintra)+=3;
1714     }
1715   else if (xtc3_context->large_inter_delta)
1716     {
1717       large_ints[0]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)])
1718         +output[outdata-natoms*3+didswap*3];
1719       large_ints[1]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)+1])
1720         +output[outdata-natoms*3+1+didswap*3];
1721       large_ints[2]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)+2])
1722         +output[outdata-natoms*3+2+didswap*3];
1723       (*ilargeinter)+=3;
1724     }
1725   prevcoord[0]=large_ints[0];
1726   prevcoord[1]=large_ints[1];
1727   prevcoord[2]=large_ints[2];
1728   output[outdata]=large_ints[0];
1729   output[outdata+1]=large_ints[1];
1730   output[outdata+2]=large_ints[2];
1731 #ifdef SHOWIT
1732   fprintf(stderr,"Unpack one large: %d %d %d\n",prevcoord[0],prevcoord[1],prevcoord[2]);
1733 #endif
1734 }
1735
1736
1737 int Ptngc_unpack_array_xtc3(unsigned char *packed,int *output, int length, int natoms)
1738 {
1739   int i;
1740   int minint[3];
1741   unsigned char *ptr=packed;
1742   int prevcoord[3];
1743   int outdata=0;
1744   int ntriplets_left=length/3;
1745   int swapatoms=0;
1746   int runlength=0;
1747   int current_large_type=0;
1748   int iinstr=0;
1749   int irle=0;
1750   int ilargedir=0;
1751   int ilargeintra=0;
1752   int ilargeinter=0;
1753   int ismallintra=0;
1754
1755   struct xtc3_context xtc3_context;
1756   init_xtc3_context(&xtc3_context);
1757
1758   for (i=0; i<3; i++)
1759     {
1760       minint[i]=unpositive_int((int)(((unsigned int)ptr[0]) |
1761                                      (((unsigned int)ptr[1])<<8) |
1762                                      (((unsigned int)ptr[2])<<16) |
1763                                      (((unsigned int)ptr[3])<<24)));
1764       ptr+=4;
1765     }
1766
1767   xtc3_context.ninstr=(int)(((unsigned int)ptr[0]) |
1768                             (((unsigned int)ptr[1])<<8) |
1769                             (((unsigned int)ptr[2])<<16) |
1770                             (((unsigned int)ptr[3])<<24));
1771   ptr+=4;
1772   if (xtc3_context.ninstr)
1773     decompress_bwlzh_block(&ptr,xtc3_context.ninstr,&xtc3_context.instructions);
1774
1775   xtc3_context.nrle=(int)(((unsigned int)ptr[0]) |
1776                           (((unsigned int)ptr[1])<<8) |
1777                           (((unsigned int)ptr[2])<<16) |
1778                           (((unsigned int)ptr[3])<<24));
1779   ptr+=4;
1780   if (xtc3_context.nrle)
1781     decompress_bwlzh_block(&ptr,xtc3_context.nrle,&xtc3_context.rle);
1782
1783   xtc3_context.nlargedir=(int)(((unsigned int)ptr[0]) |
1784                                (((unsigned int)ptr[1])<<8) |
1785                                (((unsigned int)ptr[2])<<16) |
1786                                (((unsigned int)ptr[3])<<24));
1787   ptr+=4;
1788   if (xtc3_context.nlargedir)
1789     {
1790       if (*ptr++==1)
1791         decompress_bwlzh_block(&ptr,xtc3_context.nlargedir,&xtc3_context.large_direct);
1792       else
1793         decompress_base_block(&ptr,xtc3_context.nlargedir,&xtc3_context.large_direct);
1794     }
1795
1796   xtc3_context.nlargeintra=(int)(((unsigned int)ptr[0]) |
1797                                  (((unsigned int)ptr[1])<<8) |
1798                                  (((unsigned int)ptr[2])<<16) |
1799                                  (((unsigned int)ptr[3])<<24));
1800   ptr+=4;
1801   if (xtc3_context.nlargeintra)
1802     {
1803       if (*ptr++==1)
1804         decompress_bwlzh_block(&ptr,xtc3_context.nlargeintra,&xtc3_context.large_intra_delta);
1805       else
1806         decompress_base_block(&ptr,xtc3_context.nlargeintra,&xtc3_context.large_intra_delta);
1807     }
1808
1809   xtc3_context.nlargeinter=(int)(((unsigned int)ptr[0]) |
1810                                  (((unsigned int)ptr[1])<<8) |
1811                                  (((unsigned int)ptr[2])<<16) |
1812                                  (((unsigned int)ptr[3])<<24));
1813   ptr+=4;
1814   if (xtc3_context.nlargeinter)
1815     {
1816       if (*ptr++==1)
1817         decompress_bwlzh_block(&ptr,xtc3_context.nlargeinter,&xtc3_context.large_inter_delta);
1818       else
1819         decompress_base_block(&ptr,xtc3_context.nlargeinter,&xtc3_context.large_inter_delta);
1820     }
1821
1822   xtc3_context.nsmallintra=(int)(((unsigned int)ptr[0]) |
1823                                  (((unsigned int)ptr[1])<<8) |
1824                                  (((unsigned int)ptr[2])<<16) |
1825                                  (((unsigned int)ptr[3])<<24));
1826   ptr+=4;
1827   if (xtc3_context.nsmallintra)
1828     {
1829       if (*ptr++==1)
1830         decompress_bwlzh_block(&ptr,xtc3_context.nsmallintra,&xtc3_context.smallintra);
1831       else
1832         decompress_base_block(&ptr,xtc3_context.nsmallintra,&xtc3_context.smallintra);
1833     }
1834
1835   /* Initial prevcoord is the minimum integers. */
1836   prevcoord[0]=minint[0];
1837   prevcoord[1]=minint[1];
1838   prevcoord[2]=minint[2];
1839
1840   while (ntriplets_left>0 && iinstr<xtc3_context.ninstr)
1841     {
1842       int instr=xtc3_context.instructions[iinstr++];
1843 #ifdef SHOWIT
1844       fprintf(stderr,"instr=%d @ %d\n",instr,iinstr-1);
1845 #endif
1846 #ifdef SHOWIT
1847       fprintf(stderr,"ntriplets left=%d\n",ntriplets_left);
1848 #endif
1849       if ((instr==INSTR_DEFAULT) /* large+small */
1850           || (instr==INSTR_ONLY_LARGE) /* only large */
1851           || (instr==INSTR_ONLY_SMALL)) /* only small */
1852         {
1853           if (instr!=INSTR_ONLY_SMALL)
1854             {
1855               int didswap=0;
1856               if ((instr==INSTR_DEFAULT) && (swapatoms))
1857                 didswap=1;
1858               unpack_one_large(&xtc3_context,&ilargedir, &ilargeintra, &ilargeinter,
1859                                prevcoord, minint, output, outdata, didswap,
1860                                natoms, current_large_type);
1861               ntriplets_left--;
1862               outdata+=3;
1863             }
1864           if (instr!=INSTR_ONLY_LARGE)
1865             {
1866               for (i=0; i<runlength; i++)
1867                 {
1868                   prevcoord[0]+=unpositive_int(xtc3_context.smallintra[ismallintra]);
1869                   prevcoord[1]+=unpositive_int(xtc3_context.smallintra[ismallintra+1]);
1870                   prevcoord[2]+=unpositive_int(xtc3_context.smallintra[ismallintra+2]);
1871                   ismallintra+=3;
1872                   output[outdata+i*3]=prevcoord[0];
1873                   output[outdata+i*3+1]=prevcoord[1];
1874                   output[outdata+i*3+2]=prevcoord[2];
1875 #ifdef SHOWIT
1876                   fprintf(stderr,"Unpack small: %d %d %d\n",prevcoord[0],prevcoord[1],prevcoord[2]);
1877 #endif
1878                 }
1879               if ((instr==INSTR_DEFAULT) && (swapatoms))
1880                 {
1881                   for (i=0; i<3; i++)
1882                     {
1883                       int tmp=output[outdata-3+i];
1884                       output[outdata-3+i]=output[outdata+i];
1885                       output[outdata+i]=tmp;
1886                     }
1887 #ifdef SHOWIT
1888                   fprintf(stderr,"Unswap results in\n");
1889                   for (i=0; i<3; i++)
1890                     fprintf(stderr,"%d %d %d\n",output[outdata-3+i*3],output[outdata-2+i*3],output[outdata-1+i*3]);
1891 #endif
1892                 }
1893               ntriplets_left-=runlength;
1894               outdata+=runlength*3;
1895             }
1896         }
1897       else if (instr==INSTR_LARGE_RLE && irle<xtc3_context.nrle)
1898         {
1899           int large_rle=xtc3_context.rle[irle++];
1900 #ifdef SHOWIT
1901           fprintf(stderr,"large_rle=%d @ %d\n",large_rle,irle-1);
1902 #endif
1903           for (i=0; i<large_rle; i++)
1904             {
1905               unpack_one_large(&xtc3_context,&ilargedir, &ilargeintra, &ilargeinter,
1906                                prevcoord, minint, output, outdata, 0,
1907                                natoms, current_large_type);
1908               ntriplets_left--;
1909               outdata+=3;
1910             }
1911         }
1912       else if (instr==INSTR_SMALL_RUNLENGTH && irle<xtc3_context.nrle)
1913         {
1914           runlength=xtc3_context.rle[irle++];
1915 #ifdef SHOWIT
1916           fprintf(stderr,"small_rle=%d @ %d\n",runlength,irle-1);
1917 #endif
1918         }
1919       else if (instr==INSTR_FLIP)
1920         {
1921           swapatoms=1-swapatoms;
1922 #ifdef SHOWIT
1923           fprintf(stderr,"new flip=%d\n",swapatoms);
1924 #endif
1925         }
1926       else if (instr==INSTR_LARGE_DIRECT)
1927         {
1928           current_large_type=0;
1929 #ifdef SHOWIT
1930           fprintf(stderr,"large direct\n");
1931 #endif
1932         }
1933       else if (instr==INSTR_LARGE_INTRA_DELTA)
1934         {
1935           current_large_type=1;
1936 #ifdef SHOWIT
1937           fprintf(stderr,"large intra delta\n");
1938 #endif
1939         }
1940       else if (instr==INSTR_LARGE_INTER_DELTA)
1941         {
1942           current_large_type=2;
1943 #ifdef SHOWIT
1944           fprintf(stderr,"large inter delta\n");
1945 #endif
1946         }
1947     }
1948   if (ntriplets_left<0)
1949     {
1950       fprintf(stderr,"TRAJNG XTC3: A bug has been found. At end ntriplets_left<0\n");
1951       exit(EXIT_FAILURE);
1952     }
1953   free_xtc3_context(&xtc3_context);
1954   return 0;
1955 }