Added tng_compress trajectory compression algorithms

1 files changed, 805 insertions, 0 deletions

diff --git a/src/compression/bwlzh.c b/src/compression/bwlzh.c
new file mode 100644
index 0000000..b5590e3
--- /dev/null
+++ b/src/compression/bwlzh.c
@@ -0,0 +1,805 @@
+/* This code is part of the tng compression routines.
+ *
+ * Written by Daniel Spangberg
+ * Copyright (c) 2010, 2013, The GROMACS development team.
+ *
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ */
+
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "warnmalloc.h"
+#include "tng_compress.h"
+#include "bwlzh.h"
+#include "dict.h"
+#include "vals16.h"
+#include "rle.h"
+#include "mtf.h"
+#include "bwt.h"
+#include "lz77.h"
+
+#if 0
+#define SHOWIT
+#endif
+
+#if 0
+#define SHOWTEST
+#endif
+
+#define MAX_VALS_PER_BLOCK 200000
+
+#if 1
+#define PARTIAL_MTF3
+#endif
+
+#if 0
+#define PARTIAL_MTF
+#endif
+
+int bwlzh_get_buflen(int nvals)
+{
+  return 132000+nvals*8+12*((nvals+MAX_VALS_PER_BLOCK)/MAX_VALS_PER_BLOCK);
+}
+
+#ifdef SHOWIT
+static void printvals(char *name, unsigned int *vals, int nvals)
+{
+  int i;
+  int nvalsmax=nvals;
+  if (nvalsmax>99)
+    nvalsmax=99;
+#if 0
+  for (i=0; i<nvalsmax; i++)
+    fprintf(stderr,"%s %d %u\n",name,i,vals[i]);
+#else
+  fprintf(stderr,"%s\n",name);
+  {
+    unsigned char *x=(unsigned char*)vals;
+    for (i=0; i<nvalsmax*4; i++)
+      fprintf(stderr,"%02x",(unsigned int)x[i]);
+    fprintf(stderr,"\n");
+  }
+#endif
+  
+}
+#endif
+
+
+static void bwlzh_compress_gen(unsigned int *vals, int nvals,
+			       unsigned char *output, int *output_len,
+			       int enable_lz77,
+			       int verbose)
+{
+  unsigned int *vals16;
+  int nvals16;
+  int huffdatalen;
+  int nhufflen[N_HUFFMAN_ALGO];
+  int huffalgo;
+  int ndict;
+  int bwt_index;
+  unsigned int *bwt=NULL;
+#ifdef PARTIAL_MTF3
+  unsigned char *mtf3=NULL;
+  int imtfinner;
+#endif
+  unsigned int *mtf=NULL;
+  unsigned int *rle=NULL;
+  unsigned int *offsets=NULL;
+  unsigned int *lens=NULL;
+  unsigned int *dict=warnmalloc(0x20004*sizeof *dict);
+  unsigned int *hist=warnmalloc(0x20004*sizeof *hist);
+  int nrle;
+  int noffsets;
+  int nlens;
+  unsigned char *bwlzhhuff=NULL;
+  int bwlzhhufflen;
+  int max_vals_per_block=MAX_VALS_PER_BLOCK;
+  int valsleft;
+  int thisvals;
+  int valstart;
+  int outdata=0;
+
+  unsigned int *tmpmem=warnmalloc(max_vals_per_block*18*sizeof *tmpmem);
+
+#if 0
+  verbose=1;
+#endif
+
+  bwlzhhuff=warnmalloc(Ptngc_comp_huff_buflen(3*nvals));
+  vals16=tmpmem;
+  bwt=tmpmem+max_vals_per_block*3;
+  mtf=tmpmem+max_vals_per_block*6;
+  rle=tmpmem+max_vals_per_block*9;
+  offsets=tmpmem+max_vals_per_block*12;
+  lens=tmpmem+max_vals_per_block*15;
+#ifdef PARTIAL_MTF3
+  mtf3=warnmalloc(max_vals_per_block*3*3*sizeof *mtf3); /* 3 due to expansion of 32 bit to 16 bit, 3 due to up to 3 bytes
+							   per 16 value. */
+#endif
+  if (verbose)
+    {
+      fprintf(stderr,"Number of input values: %d\n",nvals);
+    }
+
+  /* Store the number of real values in the whole block. */
+  output[outdata++]=((unsigned int)nvals)&0xFFU;
+  output[outdata++]=(((unsigned int)nvals)>>8)&0xFFU;
+  output[outdata++]=(((unsigned int)nvals)>>16)&0xFFU;
+  output[outdata++]=(((unsigned int)nvals)>>24)&0xFFU;
+  
+  valsleft=nvals;
+  valstart=0;
+  while (valsleft)
+    {
+      int reducealgo=1; /* Reduce algo is LZ77. */
+      if (!enable_lz77)
+	reducealgo=0;
+      thisvals=valsleft;
+      if (thisvals>max_vals_per_block)
+	thisvals=max_vals_per_block;
+      valsleft-=thisvals;
+      if (verbose)
+	fprintf(stderr,"Creating vals16 block from %d values.\n",thisvals);
+
+#ifdef SHOWIT
+      printvals("vals",vals+valstart,thisvals);
+#endif
+
+      Ptngc_comp_conv_to_vals16(vals+valstart,thisvals,vals16,&nvals16);
+      valstart+=thisvals;
+
+#ifdef SHOWTEST
+  nvals16=99;
+#endif
+
+#ifdef SHOWIT
+      printvals("vals16",vals16,nvals16);
+#endif
+      
+      if (verbose)
+	{
+	  fprintf(stderr,"Resulting vals16 values: %d\n",nvals16);
+	}
+      if (verbose)
+	{
+	  fprintf(stderr,"BWT\n");
+	}
+      Ptngc_comp_to_bwt(vals16,nvals16,bwt,&bwt_index);
+
+#ifdef SHOWIT
+      printvals("bwt",bwt,nvals16);
+      fprintf(stderr,"BWT INDEX is %d\n",bwt_index);
+#endif
+
+      /* Store the number of real values in this block. */
+      output[outdata++]=((unsigned int)thisvals)&0xFFU;
+      output[outdata++]=(((unsigned int)thisvals)>>8)&0xFFU;
+      output[outdata++]=(((unsigned int)thisvals)>>16)&0xFFU;
+      output[outdata++]=(((unsigned int)thisvals)>>24)&0xFFU;
+      
+      /* Store the number of nvals16 values in this block. */
+      output[outdata++]=((unsigned int)nvals16)&0xFFU;
+      output[outdata++]=(((unsigned int)nvals16)>>8)&0xFFU;
+      output[outdata++]=(((unsigned int)nvals16)>>16)&0xFFU;
+      output[outdata++]=(((unsigned int)nvals16)>>24)&0xFFU;
+
+      /* Store the BWT index. */
+      output[outdata++]=((unsigned int)bwt_index)&0xFFU;
+      output[outdata++]=(((unsigned int)bwt_index)>>8)&0xFFU;
+      output[outdata++]=(((unsigned int)bwt_index)>>16)&0xFFU;
+      output[outdata++]=(((unsigned int)bwt_index)>>24)&0xFFU;
+      
+      if (verbose)
+	fprintf(stderr,"MTF\n");
+#ifdef PARTIAL_MTF3
+      Ptngc_comp_conv_to_mtf_partial3(bwt,nvals16,
+				mtf3);
+      for (imtfinner=0; imtfinner<3; imtfinner++)
+	{
+	  int i;
+	  if (verbose)
+	    fprintf(stderr,"Doing partial MTF: %d\n",imtfinner);
+	  for (i=0; i<nvals16; i++)
+	    mtf[i]=(unsigned int)mtf3[imtfinner*nvals16+i];
+#else
+#ifdef PARTIAL_MTF
+	  Ptngc_comp_conv_to_mtf_partial(bwt,nvals16,mtf);
+#else
+	  Ptngc_comp_canonical_dict(dict,&ndict);
+	  Ptngc_comp_conv_to_mtf(bwt,nvals16,
+			   dict,ndict,mtf);
+#endif
+
+#ifdef SHOWIT
+	  printvals("mtf",mtf,nvals16);
+#endif
+#endif
+      
+
+	  if (reducealgo==1)
+	    {
+	      if (verbose)
+		fprintf(stderr,"LZ77\n");
+	      reducealgo=1;
+	      Ptngc_comp_to_lz77(mtf,nvals16,rle,&nrle,lens,&nlens,offsets,&noffsets);
+      
+	      if (verbose)
+		{
+		  fprintf(stderr,"Resulting LZ77 values: %d\n",nrle);
+		  fprintf(stderr,"Resulting LZ77 lens: %d\n",nlens);
+		  fprintf(stderr,"Resulting LZ77 offsets: %d\n",noffsets);
+		}
+#ifdef SHOWIT
+	      printvals("lz77 table",rle,nrle);
+	      printvals("lz77 lengths",lens,nlens);
+	      printvals("lz77 offsets",offsets,noffsets);
+#endif
+
+#if 0
+	      if (noffsets)
+		{
+		  unsigned int thist[0x20004];
+		  unsigned int coarse[17]={0,};
+		  int jj;
+		  Ptngc_comp_make_dict_hist(lens,nlens,dict,&ndict,thist);
+		  for (jj=0; jj<ndict; jj++)
+		    fprintf(stderr,"%d %u %u L\n",jj,dict[jj],thist[jj]);
+
+		  Ptngc_comp_make_dict_hist(offsets,noffsets,dict,&ndict,thist);
+		  for (jj=0; jj<ndict; jj++)
+		    {
+		      unsigned int v=dict[jj];
+		      int numbits=0;
+		      while (v)
+			{
+			  numbits++;
+			  v>>=1;
+			}
+		      coarse[numbits-1]+=thist[jj];
+		    }
+#if 1
+		  for (jj=0; jj<ndict; jj++)
+		    fprintf(stderr,"%d %u %u O\n",jj,dict[jj],thist[jj]);
+#else
+		  for (jj=0; jj<17; jj++)
+		    fprintf(stderr,"%d %u\n",jj+1,coarse[jj]);
+#endif
+
+		}
+	      exit(0);
+#endif
+
+	      if (nlens<2)
+		reducealgo=0;
+
+#ifdef SHOWTEST
+	      reducealgo=1;
+#endif
+	    }
+	  if (reducealgo==0)
+	    {
+	      if (verbose)
+		fprintf(stderr,"RLE\n");
+	      /* Do RLE. For any repetetitive characters. */
+	      Ptngc_comp_conv_to_rle(mtf,nvals16,rle,&nrle,1);
+	  
+#ifdef SHOWIT
+	      printvals("rle",rle,nrle);
+#endif
+	      if (verbose)
+		fprintf(stderr,"Resulting RLE values: %d\n",nrle);
+	    }
+
+	  /* reducealgo: RLE == 0, LZ77 == 1 */
+	  output[outdata++]=reducealgo;
+
+	  if (verbose)
+	    fprintf(stderr,"Huffman\n");
+      
+	  huffalgo=-1;
+	  Ptngc_comp_huff_compress_verbose(rle,nrle,bwlzhhuff,&bwlzhhufflen,&huffdatalen,nhufflen,&huffalgo,1);
+#ifdef SHOWTEST
+	  {
+	    int i;
+	    fprintf(stderr,"Huffman\n");
+	    for (i=0; i<bwlzhhufflen; i++)
+	      fprintf(stderr,"%02x",(unsigned int)bwlzhhuff[i]);
+	    fprintf(stderr,"\n");
+	    exit(0);
+	  }
+#endif
+	  if (verbose)
+	    {
+	      int i;
+	      fprintf(stderr,"Huffman data length is %d B.\n",huffdatalen);
+	      for (i=0; i<N_HUFFMAN_ALGO; i++)
+		fprintf(stderr,"Huffman dictionary for algorithm %s is %d B.\n",Ptngc_comp_get_huff_algo_name(i),nhufflen[i]-huffdatalen);
+	      fprintf(stderr,"Resulting algorithm: %s. Size=%d B\n",Ptngc_comp_get_huff_algo_name(huffalgo),bwlzhhufflen);
+	    }
+
+	  /* Store the number of huffman values in this block. */
+	  output[outdata++]=((unsigned int)nrle)&0xFFU;
+	  output[outdata++]=(((unsigned int)nrle)>>8)&0xFFU;
+	  output[outdata++]=(((unsigned int)nrle)>>16)&0xFFU;
+	  output[outdata++]=(((unsigned int)nrle)>>24)&0xFFU;
+
+	  /* Store the size of the huffman block. */
+	  output[outdata++]=((unsigned int)bwlzhhufflen)&0xFFU;
+	  output[outdata++]=(((unsigned int)bwlzhhufflen)>>8)&0xFFU;
+	  output[outdata++]=(((unsigned int)bwlzhhufflen)>>16)&0xFFU;
+	  output[outdata++]=(((unsigned int)bwlzhhufflen)>>24)&0xFFU;
+
+	  /* Store the huffman block. */
+	  memcpy(output+outdata,bwlzhhuff,bwlzhhufflen);
+	  outdata+=bwlzhhufflen;
+
+	  if (reducealgo==1)
+	    {
+	      /* Store the number of values in this block. */
+	      output[outdata++]=((unsigned int)noffsets)&0xFFU;
+	      output[outdata++]=(((unsigned int)noffsets)>>8)&0xFFU;
+	      output[outdata++]=(((unsigned int)noffsets)>>16)&0xFFU;
+	      output[outdata++]=(((unsigned int)noffsets)>>24)&0xFFU;
+
+	      if (noffsets>0)
+		{
+		  if (verbose)
+		    fprintf(stderr,"Huffman for offsets\n");
+      
+		  huffalgo=-1;
+		  Ptngc_comp_huff_compress_verbose(offsets,noffsets,bwlzhhuff,&bwlzhhufflen,&huffdatalen,nhufflen,&huffalgo,1);
+		  if (verbose)
+		    {
+		      int i;
+		      fprintf(stderr,"Huffman data length is %d B.\n",huffdatalen);
+		      for (i=0; i<N_HUFFMAN_ALGO; i++)
+			fprintf(stderr,"Huffman dictionary for algorithm %s is %d B.\n",Ptngc_comp_get_huff_algo_name(i),nhufflen[i]-huffdatalen);
+		      fprintf(stderr,"Resulting algorithm: %s. Size=%d B\n",Ptngc_comp_get_huff_algo_name(huffalgo),bwlzhhufflen);
+		    }
+
+		  /* If huffman was bad for these offsets, just store the offsets as pairs of bytes. */
+		  if (bwlzhhufflen<noffsets*2)
+		    {
+		      output[outdata++]=0;
+
+		      /* Store the size of the huffman block. */
+		      output[outdata++]=((unsigned int)bwlzhhufflen)&0xFFU;
+		      output[outdata++]=(((unsigned int)bwlzhhufflen)>>8)&0xFFU;
+		      output[outdata++]=(((unsigned int)bwlzhhufflen)>>16)&0xFFU;
+		      output[outdata++]=(((unsigned int)bwlzhhufflen)>>24)&0xFFU;
+
+		      /* Store the huffman block. */
+		      memcpy(output+outdata,bwlzhhuff,bwlzhhufflen);
+		      outdata+=bwlzhhufflen;
+		    }
+		  else
+		    {
+		      int i;
+		      output[outdata++]=1;
+		      for (i=0; i<noffsets; i++)
+			{
+			  output[outdata++]=((unsigned int)offsets[i])&0xFFU;
+			  output[outdata++]=(((unsigned int)offsets[i])>>8)&0xFFU;
+			}
+		      if (verbose)
+			fprintf(stderr,"Store raw offsets: %d B\n",noffsets*2);
+		    }
+		}
+
+#if 0
+	      {
+		int i,ndict;
+		FILE *f=fopen("len.dict","w");
+		Ptngc_comp_make_dict_hist(lens,nlens,dict,&ndict,hist);
+		for (i=0; i<ndict; i++)
+		  fprintf(f,"%d %d %d\n",i,dict[i],hist[i]);
+		fclose(f);
+		f=fopen("off.dict","w");
+		Ptngc_comp_make_dict_hist(offsets,noffsets,dict,&ndict,hist);
+		for (i=0; i<ndict; i++)
+		  fprintf(f,"%d %d %d\n",i,dict[i],hist[i]);
+		fclose(f);
+		f=fopen("len.time","w");
+		for (i=0; i<ndict; i++)
+		  fprintf(f,"%d\n",lens[i]);
+		fclose(f);
+		f=fopen("off.time","w");
+		for (i=0; i<ndict; i++)
+		  fprintf(f,"%d\n",offsets[i]);
+		fclose(f);
+	      }
+#endif
+
+	      if (verbose)
+		fprintf(stderr,"Huffman for lengths\n");
+      
+	      huffalgo=-1;
+	      Ptngc_comp_huff_compress_verbose(lens,nlens,bwlzhhuff,&bwlzhhufflen,&huffdatalen,nhufflen,&huffalgo,1);
+	      if (verbose)
+		{
+		  int i;
+		  fprintf(stderr,"Huffman data length is %d B.\n",huffdatalen);
+		  for (i=0; i<N_HUFFMAN_ALGO; i++)
+		    fprintf(stderr,"Huffman dictionary for algorithm %s is %d B.\n",Ptngc_comp_get_huff_algo_name(i),nhufflen[i]-huffdatalen);
+		  fprintf(stderr,"Resulting algorithm: %s. Size=%d B\n",Ptngc_comp_get_huff_algo_name(huffalgo),bwlzhhufflen);
+		}
+
+	      /* Store the number of values in this block. */
+	      output[outdata++]=((unsigned int)nlens)&0xFFU;
+	      output[outdata++]=(((unsigned int)nlens)>>8)&0xFFU;
+	      output[outdata++]=(((unsigned int)nlens)>>16)&0xFFU;
+	      output[outdata++]=(((unsigned int)nlens)>>24)&0xFFU;
+	  
+	      /* Store the size of the huffman block. */
+	      output[outdata++]=((unsigned int)bwlzhhufflen)&0xFFU;
+	      output[outdata++]=(((unsigned int)bwlzhhufflen)>>8)&0xFFU;
+	      output[outdata++]=(((unsigned int)bwlzhhufflen)>>16)&0xFFU;
+	      output[outdata++]=(((unsigned int)bwlzhhufflen)>>24)&0xFFU;
+	  
+	      /* Store the huffman block. */
+	      memcpy(output+outdata,bwlzhhuff,bwlzhhufflen);
+	      outdata+=bwlzhhufflen;
+	    }
+#ifdef PARTIAL_MTF3
+	}
+#endif
+    }
+
+  *output_len=outdata;
+  free(hist);
+  free(dict);
+  free(bwlzhhuff);
+#ifdef PARTIAL_MTF3
+  free(mtf3);
+#endif
+  free(tmpmem);
+}
+
+
+void DECLSPECDLLEXPORT bwlzh_compress(unsigned int *vals, int nvals,
+		  unsigned char *output, int *output_len)
+{
+  bwlzh_compress_gen(vals,nvals,output,output_len,1,0);
+}
+
+void DECLSPECDLLEXPORT bwlzh_compress_verbose(unsigned int *vals, int nvals,
+			  unsigned char *output, int *output_len)
+{
+  bwlzh_compress_gen(vals,nvals,output,output_len,1,1);
+}
+
+
+void DECLSPECDLLEXPORT bwlzh_compress_no_lz77(unsigned int *vals, int nvals,
+		  unsigned char *output, int *output_len)
+{
+  bwlzh_compress_gen(vals,nvals,output,output_len,0,0);
+}
+
+void DECLSPECDLLEXPORT bwlzh_compress_no_lz77_verbose(unsigned int *vals, int nvals,
+			  unsigned char *output, int *output_len)
+{
+  bwlzh_compress_gen(vals,nvals,output,output_len,0,1);
+}
+
+
+static void bwlzh_decompress_gen(unsigned char *input, int nvals,
+			       unsigned int *vals,
+			       int verbose)
+{
+  unsigned int *vals16;
+  int nvals16;
+  int ndict;
+  int bwt_index;
+  unsigned int *bwt=NULL; 
+  unsigned int *mtf=NULL;
+#ifdef PARTIAL_MTF3
+  unsigned char *mtf3=NULL;
+  int imtfinner;
+#endif
+  unsigned int *rle=NULL;
+  unsigned int *offsets=NULL;
+  unsigned int *lens=NULL;
+  unsigned int *dict=warnmalloc(0x20004*sizeof *dict);
+  unsigned int *hist=warnmalloc(0x20004*sizeof *hist);
+  int nrle, noffsets, nlens;
+  unsigned char *bwlzhhuff=NULL;
+  int bwlzhhufflen;
+  int max_vals_per_block=MAX_VALS_PER_BLOCK;
+  int valsleft;
+  int thisvals;
+  int valstart;
+  int inpdata=0;
+  int nvalsfile;
+
+  unsigned int *tmpmem=warnmalloc(max_vals_per_block*18*sizeof *tmpmem);
+
+#if 0
+  verbose=1;
+#endif
+
+
+  bwlzhhuff=warnmalloc(Ptngc_comp_huff_buflen(3*nvals));
+  vals16=tmpmem;
+  bwt=tmpmem+max_vals_per_block*3;
+  mtf=tmpmem+max_vals_per_block*6;
+  rle=tmpmem+max_vals_per_block*9;
+  offsets=tmpmem+max_vals_per_block*12;
+  lens=tmpmem+max_vals_per_block*15;
+#ifdef PARTIAL_MTF3
+  mtf3=warnmalloc(max_vals_per_block*3*3*sizeof *mtf3); /* 3 due to expansion of 32 bit to 16 bit, 3 due to up to 3 bytes
+							   per 16 value. */
+#endif
+
+  if (verbose)
+    {
+      fprintf(stderr,"Number of input values: %d\n",nvals);
+    }
+
+  /* Read the number of real values in the whole block. */
+  nvalsfile=(int)(((unsigned int)input[inpdata]) |
+		  (((unsigned int)input[inpdata+1])<<8) |
+		  (((unsigned int)input[inpdata+2])<<16) |
+		  (((unsigned int)input[inpdata+3])<<24));
+  inpdata+=4;
+
+  if (nvalsfile!=nvals)
+    {
+      fprintf(stderr,"BWLZH: The number of values found in the file is different from the number of values expected.\n");
+      exit(EXIT_FAILURE);
+    }
+  
+  valsleft=nvals;
+  valstart=0;
+  while (valsleft)
+    {
+      int valsnew;
+      int reducealgo;
+      /* Read the number of real values in this block. */
+      thisvals=(int)(((unsigned int)input[inpdata]) |
+		     (((unsigned int)input[inpdata+1])<<8) |
+		     (((unsigned int)input[inpdata+2])<<16) |
+		     (((unsigned int)input[inpdata+3])<<24));
+      inpdata+=4;
+      
+      valsleft-=thisvals;
+
+      /* Read the number of nvals16 values in this block. */
+      nvals16=(int)(((unsigned int)input[inpdata]) |
+		    (((unsigned int)input[inpdata+1])<<8) |
+		    (((unsigned int)input[inpdata+2])<<16) |
+		    (((unsigned int)input[inpdata+3])<<24));
+      inpdata+=4;
+
+      /* Read the BWT index. */
+      bwt_index=(int)(((unsigned int)input[inpdata]) |
+		      (((unsigned int)input[inpdata+1])<<8) |
+		      (((unsigned int)input[inpdata+2])<<16) |
+		      (((unsigned int)input[inpdata+3])<<24));
+      inpdata+=4;
+
+      if (thisvals>max_vals_per_block)
+	{
+	  /* More memory must be allocated for decompression. */
+	  max_vals_per_block=thisvals;
+	  if (verbose)
+	    fprintf(stderr,"Allocating more memory: %d B\n",(int)(max_vals_per_block*15*sizeof *tmpmem));
+	  tmpmem=warnrealloc(tmpmem,max_vals_per_block*18*sizeof *tmpmem);
+	  vals16=tmpmem;
+	  bwt=tmpmem+max_vals_per_block*3;
+	  mtf=tmpmem+max_vals_per_block*6;
+	  rle=tmpmem+max_vals_per_block*9;
+	  offsets=tmpmem+max_vals_per_block*12;
+	  lens=tmpmem+max_vals_per_block*15;
+#ifdef PARTIAL_MTF3
+	  mtf3=warnrealloc(mtf3,max_vals_per_block*3*3*sizeof *mtf3); /* 3 due to expansion of 32 bit to 16 bit, 3 due to up to 3 bytes
+									 per 16 value. */
+#endif	  
+	}
+
+#ifdef PARTIAL_MTF3
+      for (imtfinner=0; imtfinner<3; imtfinner++)
+	{
+	  int i;
+	  if (verbose)
+	    fprintf(stderr,"Doing partial MTF: %d\n",imtfinner);
+#endif
+
+	  reducealgo=(int)input[inpdata];
+	  inpdata++;
+
+	  /* Read the number of huffman values in this block. */
+	  nrle=(int)(((unsigned int)input[inpdata]) |
+		     (((unsigned int)input[inpdata+1])<<8) |
+		     (((unsigned int)input[inpdata+2])<<16) |
+		     (((unsigned int)input[inpdata+3])<<24));
+	  inpdata+=4;
+
+	  /* Read the size of the huffman block. */
+	  bwlzhhufflen=(int)(((unsigned int)input[inpdata]) |
+			     (((unsigned int)input[inpdata+1])<<8) |
+			     (((unsigned int)input[inpdata+2])<<16) |
+			     (((unsigned int)input[inpdata+3])<<24));
+	  inpdata+=4;
+
+	  if (verbose)
+	    fprintf(stderr,"Decompressing huffman block of length %d.\n",bwlzhhufflen);
+	  /* Decompress the huffman block. */
+	  Ptngc_comp_huff_decompress(input+inpdata,bwlzhhufflen,rle);
+	  inpdata+=bwlzhhufflen;
+
+	  if (reducealgo==1) /* LZ77 */
+	    {
+	      int offstore;
+	      /* Read the number of huffman values in this block. */
+	      noffsets=(int)(((unsigned int)input[inpdata]) |
+			     (((unsigned int)input[inpdata+1])<<8) |
+			     (((unsigned int)input[inpdata+2])<<16) |
+			     (((unsigned int)input[inpdata+3])<<24));
+	      inpdata+=4;
+
+	      if (noffsets>0)
+		{
+		  /* How are the offsets stored? */
+		  offstore=(int)input[inpdata++];
+		  if (offstore==0)
+		    {
+		      /* Read the size of the huffman block. */
+		      bwlzhhufflen=(int)(((unsigned int)input[inpdata]) |
+					 (((unsigned int)input[inpdata+1])<<8) |
+					 (((unsigned int)input[inpdata+2])<<16) |
+					 (((unsigned int)input[inpdata+3])<<24));
+		      inpdata+=4;
+		  
+		      if (verbose)
+			fprintf(stderr,"Decompressing offset huffman block.\n");
+		  
+		      /* Decompress the huffman block. */
+		      Ptngc_comp_huff_decompress(input+inpdata,bwlzhhufflen,offsets);
+		      inpdata+=bwlzhhufflen;
+		    }
+		  else
+		    {
+		      int i;
+		      if (verbose)
+			fprintf(stderr,"Reading offset block.\n");
+		      for (i=0; i<noffsets; i++)
+			{
+			  offsets[i]=(int)(((unsigned int)input[inpdata]) |
+					   (((unsigned int)input[inpdata+1])<<8));
+			  inpdata+=2;
+			}
+		    }
+		}
+#if 0
+	      {
+		int i;
+		for (i=0; i<nrle; i++)
+		  fprintf(stderr,"RLE %d: %d\n",i,rle[i]);
+		for (i=0; i<noffsets; i++)
+		  fprintf(stderr,"OFFSET %d: %d\n",i,offsets[i]);
+	      }
+#endif
+
+
+	      /* Read the number of huffman values in this block. */
+	      nlens=(int)(((unsigned int)input[inpdata]) |
+			  (((unsigned int)input[inpdata+1])<<8) |
+			  (((unsigned int)input[inpdata+2])<<16) |
+			  (((unsigned int)input[inpdata+3])<<24));
+	      inpdata+=4;
+
+	      /* Read the size of the huffman block. */
+	      bwlzhhufflen=(int)(((unsigned int)input[inpdata]) |
+				 (((unsigned int)input[inpdata+1])<<8) |
+				 (((unsigned int)input[inpdata+2])<<16) |
+				 (((unsigned int)input[inpdata+3])<<24));
+	      inpdata+=4;
+	  
+	      if (verbose)
+		fprintf(stderr,"Decompressing length huffman block.\n");
+	  
+	      /* Decompress the huffman block. */
+	      Ptngc_comp_huff_decompress(input+inpdata,bwlzhhufflen,lens);
+	      inpdata+=bwlzhhufflen;
+
+	      if (verbose)
+		fprintf(stderr,"Decompressing LZ77.\n");
+	  
+	      Ptngc_comp_from_lz77(rle,nrle,lens,nlens,offsets,noffsets,mtf,nvals16);
+	    }
+	  else if (reducealgo==0) /* RLE */
+	    {
+#ifdef SHOWIT
+	      printvals("rle",rle,nrle);
+#endif
+	  
+	      if (verbose)
+		fprintf(stderr,"Decompressing rle block.\n");
+	      Ptngc_comp_conv_from_rle(rle,mtf,nvals16);
+	    }
+
+#ifdef PARTIAL_MTF3
+	  for (i=0; i<nvals16; i++)
+	    mtf3[imtfinner*nvals16+i]=(unsigned char)mtf[i];
+	}
+#else
+#ifdef SHOWIT
+      printvals("mtf",mtf,nvals16);
+#endif
+      
+#endif
+
+
+      if (verbose)
+	fprintf(stderr,"Inverse MTF.\n");
+#ifdef PARTIAL_MTF3
+      Ptngc_comp_conv_from_mtf_partial3(mtf3,nvals16,bwt);
+#else
+#ifdef PARTIAL_MTF
+      Ptngc_comp_conv_from_mtf_partial(mtf,nvals16,bwt);
+#else
+      Ptngc_comp_canonical_dict(dict,&ndict);      
+      Ptngc_comp_conv_from_mtf(mtf,nvals16,dict,ndict,bwt);
+#endif
+#endif
+
+#ifdef SHOWIT
+      printvals("bwt",bwt,nvals16);
+      fprintf(stderr,"BWT INDEX is %d\n",bwt_index);
+#endif
+      
+      if (verbose)
+	fprintf(stderr,"Inverse BWT.\n");
+      Ptngc_comp_from_bwt(bwt,nvals16,bwt_index,vals16);
+
+#ifdef SHOWIT
+      printvals("vals16",vals16,nvals16);
+#endif
+
+      if (verbose)
+	fprintf(stderr,"Decompressing vals16 block.\n");
+      Ptngc_comp_conv_from_vals16(vals16,nvals16,vals+valstart,&valsnew);
+
+#ifdef SHOWIT
+      printvals("vals",vals+valstart,thisvals);
+#endif
+
+      if (valsnew!=thisvals)
+	{
+	  fprintf(stderr,"BWLZH: Block contained different number of values than expected.\n");
+	  exit(EXIT_FAILURE);
+	}
+      valstart+=thisvals;
+    }
+  free(hist);
+  free(dict);
+  free(bwlzhhuff);
+#ifdef PARTIAL_MTF3
+  free(mtf3);
+#endif
+  free(tmpmem);
+}
+
+
+void DECLSPECDLLEXPORT bwlzh_decompress(unsigned char *input, int nvals,
+		    unsigned int *vals)
+{
+  bwlzh_decompress_gen(input,nvals,vals,0);
+}
+
+void DECLSPECDLLEXPORT bwlzh_decompress_verbose(unsigned char *input, int nvals,
+			    unsigned int *vals)
+{
+  bwlzh_decompress_gen(input,nvals,vals,1);
+}
+