/* crunch.c, begun 12/6/91 using fragments from earlier modules */
#include<stdio.h>
#if NeXT
#else
#include  <malloc.h>
#endif
 typedef unsigned char byte;
 byte bits[8]={1,2,4,8,16,32,64,128};
 /*--------------------------------------------------------------------------*/
int anacrunch(x,array,slice,nx,ny, limit)
 /* compress 16 bit array into x (a byte array) using ny blocks each of size
 nx, bit slice size slice, returns # of bytes in x */
	 unsigned char *x;
	 int nx,ny,slice, limit;
	 short array[];
 {
 void	swapl();
 struct compresshead {
 int     tsize,nblocks,bsize;
 byte    slice_size,type; } *ch;
 short iq;
 unsigned nb,ixa,ixb;
 unsigned register i,j,r1,in;
 int r0,r2,r3,r4,r5,r6,r7,r8,mask,fac;
 int i2,k,ix,iy;
 unsigned char xq;
 union { int i; short w; unsigned char b[4]; } y;
 /* begin execution */
 limit = limit - 10;	/* some margin since we don't check all times */
 mask=1; for (i=0;i<slice;i++) mask=2*mask;
 fac=mask;       mask=mask-1; /* no inline expon. in C */
 /* determine the # of bytes to transfer to 32 bit int for fixed portion */
 if (slice == 0) nb=0; else { if (slice < 2 ) nb=1;
	 else { if (slice < 10) nb=2; else nb=3;    }};
 y.i=0;
 /* do the compression header */
 ch = (struct compresshead *) x;
 /* important note - can't use the sizeof(struct compresshead) because it
	 is 14 on some machines and rounded up to 16 on others */
 /*x = x + sizeof(struct compresshead);*/
 x = x + 14;
 ch->bsize = nx;  ch->nblocks = ny;  ch->slice_size = slice;  ch->type = 0;
 i=0;    r1=0;   in=0;
 for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
 /* load the first value, reverse bytes (VAX style)*/
#if LITTLEENDIAN
 y.w=array[in]   ;x[i]=y.b[0]    ;x[i+1]=y.b[1];
#else
 y.w=array[in]   ;x[i]=y.b[1]    ;x[i+1]=y.b[0];
#endif 
 /* in=in+1;*/
 r1=r1+16;
 r2=0;
 ixa=1+iy*nx;    ixb=(iy+1)*nx;
 /*for (ix=1; ix<nx; ix++)       {       */      /* start of ix (inner) loop */
 for (in=ixa; in<ixb; in++)      {               /* start of ix (inner) loop */
 /* first the fixed slice portion */
 y.i=array[in]-array[in-1];
 r3=(y.i>>slice);
 i=r1>>3;
 j=r1%8;
 if ( i > limit ) return -1;		/* bad news, went too far */
 /* now load nb bytes into x */
 /*low order byte of y.i is first in stream */
#if LITTLEENDIAN
 if (j == 0) {y.i=(y.i & mask); x[i]=y.b[0];}
	 else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[0];}
 if (nb>1) { x[i+1]=y.b[1]; if (nb>2) x[i+2]=y.b[2]; }
#else
 if (j == 0) {y.i=(y.i & mask); x[i]=y.b[3];}
	 else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[3];}
 if (nb>1) { x[i+1]=y.b[2]; if (nb>2) x[i+2]=y.b[1]; }
 
#endif 
 r1=r1+slice;       /* bump r1 pass the fixed part */
 i=r1>>3;                j=r1%8;
 /* note that r3 is the # of bits required minus 1 */
 if (r3==0) { if (j ==0 ) {x[i]=bits[j];} else {x[i]=x[i]|bits[j];}
 r1+=1;} else    {
 r3=2*r3;        if (r3<0) r3 = -r3-1;
 if (r3<31)  {
 r0=j+r3;        /* this is the bit that needs setting offset from x[i] */
 if (r0 < 8) { if (j == 0) x[i]=bits[r0]; else x[i]=x[i]|bits[r0];}
 else {if (j == 0) x[i]=0;  j=r0%8; if (r0 < 16) x[i+1]=bits[j];
 else { i2=i+r0/8; for (k=i+1;k<i2;k++) x[k]=0;  x[i2]=bits[j]; }
 }
 r1+=1+r3;
 } else {        /* big one exception, should be rare */
 /* does not need to be efficient, used rarely */
 /*printf("big one \n");*/
 if (j == 0) x[i]=0;     /* gotta zero the virgins */
 r0=j+31;        j=r0%8;         i2=i+r0/8;
 for (k=i+1;k<i2;k++) x[k]=0;    x[i2]=bits[j];
 /* recompute the difference and load 17 bits (always 3 bytes) */
 r1=r1+32;
 i=r1/8;
 j=r1%8;
 if (j == 0) x[i]=0;     /* gotta zero the virgins */
 y.i=((array[in]-array[in-1])& 0x1ffff) << j;
#if LITTLEENDIAN
 x[i]=x[i] | y.b[0]; x[i+1]=y.b[1];      x[i+2]=y.b[2];
#else
 x[i]=x[i] | y.b[3]; x[i+1]=y.b[2];      x[i+2]=y.b[1];
#endif
 r1=r1+17;       } /* end of big one exception */
		 }       /* end of (r3==0) conditional */
 /*in=in+1; */                   }       /* end of ix loop */
 /* some checks here */
 /* bump to next byte boundary */
 i=(r1+7)/8;     r1=8*i;                 }       /* end of iy loop */
 ch->tsize = i = i + 14;
 /* we have to put these in a form readable by the Vax (these may be used
	 by fcwrite) */
#if LITTLEENDIAN
#else
 swapl(&(ch->tsize),1); swapl(&(ch->bsize),1); swapl(&(ch->nblocks),1); 
#endif
 return  i;      /*return # of bytes used */
 }       /* end of routine */
 /*--------------------------------------------------------------------------*/
int anacrunch8(x,array,slice,nx,ny, limit)
 /* compress 8 bit array into x (a byte array) using ny blocks each of size
 nx, bit slice size slice, returns # of bytes in x */
 unsigned char *x;
 int nx,ny,slice, limit;
 byte array[];
 {
 void	swapl();
 struct compresshead {
 int     tsize,nblocks,bsize;
 byte    slice_size,type; } *ch;
 short iq;
 unsigned nb,ixa,ixb;
 unsigned register i,j,r1,in;
 int r0,r2,r3,r4,r5,r6,r7,r8,mask,fac;
 int i2,k,ix,iy;
 unsigned char xq;
 union { int i; short w; unsigned char b[4]; } y;
 /* begin execution */
 limit = limit - 10;	/* some margin since we don't check all times */
 mask=1; for (i=0;i<slice;i++) mask=2*mask;
 fac=mask;       mask=mask-1; /* no inline expon. in C */
 /* determine the # of bytes to transfer to 32 bit int for fixed portion */
 if (slice > 8) slice = 8;
 if (slice == 0) nb=0; else { if (slice < 2 ) nb=1;
	 else { if (slice < 10) nb=2; else nb=3;    }};
 y.i=0;
 /* do the compression header */
 ch = (struct compresshead *) x;
 /* important note - can't use the sizeof(struct compresshead) because it
	 is 14 on some machines and rounded up to 16 on others */
 /*x = x + sizeof(struct compresshead);*/
 x = x + 14;
 ch->bsize = nx;  ch->nblocks = ny;  ch->slice_size = slice;  ch->type = 1;
 i=0;    r1=0;   in=0;
 for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
 /* load the first value */
 x[i] = array[in];
 r1=r1+8;
 r2=0;
 ixa=1+iy*nx;    ixb=(iy+1)*nx;
 for (in=ixa; in<ixb; in++)      {               /* start of ix (inner) loop */
 /* first the fixed slice portion */
 y.i= (int) array[in]- (int) array[in-1];
 r3=(y.i>>slice);
 i=r1>>3;
 j=r1%8;
 if ( i > limit ) return -1;		/* bad news, went too far */
 /* now load nb bytes into x */
 /*low order byte of y.i is first in stream */
#if LITTLEENDIAN
 if (j == 0) {y.i=(y.i & mask); x[i]=y.b[0];}
	 else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[0];}
 if (nb>1) { x[i+1]=y.b[1]; }
#else
 if (j == 0) {y.i=(y.i & mask); x[i]=y.b[3];}
	 else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[3];}
 if (nb>1) { x[i+1]=y.b[2]; }
 
#endif 
 r1=r1+slice;       /* bump r1 pass the fixed part */
 i=r1>>3;                j=r1%8;
 /* note that r3 is the # of bits required minus 1 */
 if (r3==0) { if (j ==0 ) {x[i]=bits[j];} else {x[i]=x[i]|bits[j];}
 r1+=1;} else    {
 r3=2*r3;        if (r3<0) r3 = -r3-1;
 if (r3<31)  {
 r0=j+r3;        /* this is the bit that needs setting offset from x[i] */
 if (r0 < 8) { if (j == 0) x[i]=bits[r0]; else x[i]=x[i]|bits[r0];}
 else {if (j == 0) x[i]=0;  j=r0%8; if (r0 < 16) x[i+1]=bits[j];
 else { i2=i+r0/8; for (k=i+1;k<i2;k++) x[k]=0;  x[i2]=bits[j]; }
 }
 r1+=1+r3;
 } else {        /* big one exception, should be rare */
 /* does not need to be efficient, used rarely */
 /* printf("big one \n"); */
 if (j == 0) x[i]=0;     /* gotta zero the virgins */
 r0=j+31;        j=r0%8;         i2=i+r0/8;
 for (k=i+1;k<i2;k++) x[k]=0;    x[i2]=bits[j];
 /* recompute the difference and load 9 bits (always 2 bytes) */
 r1=r1+32;
 i=r1/8;
 j=r1%8;
 if (j == 0) x[i]=0;     /* gotta zero the virgins */
 y.i=((array[in]-array[in-1])& 0x1ff) << j;
#if LITTLEENDIAN
 x[i]=x[i] | y.b[0]; x[i+1]=y.b[1];
#else
 x[i]=x[i] | y.b[3]; x[i+1]=y.b[2];
#endif
 r1=r1+9;       } /* end of big one exception */
		 }       /* end of (r3==0) conditional */
		 }       /* end of ix loop */
 /* some checks here */
 /* bump to next byte boundary */
 i=(r1+7)/8;     r1=8*i;                 }       /* end of iy loop */
 ch->tsize = i = i + 14;
 /* we have to put these in a form readable by the Vax (these may be used
	 by fcwrite) */
#if LITTLEENDIAN
#else
 swapl(&(ch->tsize),1); swapl(&(ch->bsize),1); swapl(&(ch->nblocks),1); 
#endif
 return  i;      /*return # of bytes used */
 }       /* end of routine */
 /*--------------------------------------------------------------------------*/
int anacrunch32(x,array,slice,nx,ny, limit)
 /* compress 32 bit array into x (a byte array) using ny blocks each of size
 nx, bit slice size slice, returns # of bytes in x */
	 unsigned char *x;
	 int nx,ny,slice, limit;
	 int array[];
 {
 void	swapl();
 struct compresshead {
 int     tsize,nblocks,bsize;
 byte    slice_size,type; } *ch;
 short iq;
 unsigned int nb,ixa,ixb,big=0;
 unsigned register i,j,r1,in;
 int r0,r2,r4,fac;
 long long	r3, mask, y64;
 int i2,k,ix,iy;
 unsigned char xq;
 union { int i; short w; unsigned char b[4]; } y;
 union { long long l64; unsigned char b[8];  } yy;
 /* begin execution */
 limit = limit - 10;	/* some margin since we don't check all times */
 mask=1; for (i=0;i<slice;i++) mask=2*mask;
 fac=mask;       mask=mask-1; /* no inline expon. in C */
 /* determine the # of bytes to transfer to 32 bit int for fixed portion */
 nb = (slice + 14)/8;	/* range 1 to 5 */
 if (slice == 0) nb=0;	/* but slice = 0 a special case */

 y.i=0;
 /* do the compression header */
 ch = (struct compresshead *) x;
 /* important note - can't use the sizeof(struct compresshead) because it
	 is 14 on some machines and rounded up to 16 on others */
 /*x = x + sizeof(struct compresshead);*/
 x = x + 14;
 ch->bsize = nx;  ch->nblocks = ny;  ch->slice_size = slice;  ch->type = 4;
 i=0;    r1=0;   in=0;
 for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
 /* load the first value, reverse bytes (VAX style)*/
#if LITTLEENDIAN
 y.i=array[in]; x[i]=y.b[0]; x[i+1]=y.b[1]; x[i+2]=y.b[2]; x[i+3]=y.b[3];
#else
 y.i=array[in]; x[i]=y.b[3]; x[i+1]=y.b[2]; x[i+2]=y.b[1]; x[i+3]=y.b[0];
#endif 
 r1=r1+32;
 r2=0;
 ixa=1+iy*nx;    ixb=(iy+1)*nx;
 for (in=ixa; in<ixb; in++)      {               /* start of ix (inner) loop */
 /* first the fixed slice portion */
 y64 = (long long) array[in] - (long long) array[in-1];
 r3 = (y64>>slice);
 i=r1>>3;
 j=r1%8;
 if ( i > limit ) return -1;		/* bad news, went too far */
 /* now load nb bytes into x */
 /*low order byte of y.i is first in stream */
 if (j == 0) {
	y64=(y64 & mask);	x[i]= (byte) y64;
	/* since we started at bit 0, spillover to the next byte is
	determined as follows (and is unlikely since slice gt 8 is unusual */
	if (slice > 8) { x[i+1]= (byte) (y64 >> 8);
	 if (slice > 16) { x[i+2]= (byte) (y64 >> 16);
	  if (slice > 24) { x[i+3]= (byte) (y64 >> 24);}}}
} else {
	y64=(y64 & mask)<<j;	x[i]=x[i] | (byte) y64;
	/* spillover more likely here */
	if (nb>1) { x[i+1] = (byte) (y64 >> 8);
	 if (nb>2) { x[i+2] = (byte) (y64 >> 16);
	  if (nb>3) { x[i+3] = (byte) (y64 >> 24);
	   if (nb>4) { x[i+4] = (byte) (y64 >> 32);}}}}
}

 r1=r1+slice;       /* bump r1 pass the fixed part */
 i=r1>>3;                j=r1%8;
 /* note that r3 is the # of bits required minus 1 */
 if (r3==0) { if (j ==0 ) {x[i]= *bits;} else {x[i]=x[i]| bits[j];}
 r1++;} else    {
 r3=2*r3;        if (r3<0) r3 = -r3-1;
 if (r3<31)  {
 r0=j+r3;        /* this is the bit that needs setting offset from x[i] */
 if (r0 < 8) { if (j == 0) x[i]=bits[r0]; else x[i]=x[i]|bits[r0];}
   else  {
    if (j == 0) x[i]=0;
     j=r0%8;
     if (r0 < 16) x[i+1]=bits[j]; else {
 	i2=i+r0/8;
	for (k=i+1;k<i2;k++) x[k]=0;
	x[i2]=bits[j]; }
   }
   r1+=1+r3;
 } else {        /* big one exception, should be rare */
 /* does not need to be efficient, used rarely */
 /* printf("big one for I*4\n"); */
 big++;
 if (j == 0) x[i]=0;     /* gotta zero the virgins */
 r0=j+31;        j=r0%8;         i2=i+r0/8;
 for (k=i+1;k<i2;k++) x[k]=0;    x[i2]=bits[j];
 /* recompute the difference and load 33 bits (always 5 bytes) */
 r1=r1+32;
 i=r1/8;
 j=r1%8;
 if (j == 0) x[i]=0;     /* gotta zero the virgins */
 yy.l64=(((long long) array[in]- (long long) array[in-1])& 0x1ffffffff) << j;
#if LITTLEENDIAN
 x[i]=x[i] | yy.b[0]; x[i+1]=yy.b[1]; x[i+2]=yy.b[2];
 x[i+3]=yy.b[3]; x[i+4]=yy.b[4];
#else
 x[i]=x[i] | yy.b[7]; x[i+1]=yy.b[6]; x[i+2]=yy.b[5];
 x[i+3]=yy.b[4]; x[i+4]=yy.b[3];
#endif
 r1=r1+33;       } /* end of big one exception */
		 }       /* end of (r3==0) conditional */
 /*in=in+1; */                   }       /* end of ix loop */
 /* some checks here */
 /* bump to next byte boundary */
 i=(r1+7)/8;     r1=8*i;                 }       /* end of iy loop */
 ch->tsize = i = i + 14;
 /* we have to put these in a form readable by the Vax (these may be used
	 by fcwrite) */
#if LITTLEENDIAN
#else
 swapl(&(ch->tsize),1); swapl(&(ch->bsize),1); swapl(&(ch->nblocks),1); 
#endif
 /* printf("number of big ones for this I*4 = %d\n", big); */
 return  i;      /*return # of bytes used */
 }       /* end of routine */
 /*--------------------------------------------------------------------------*/
void swapl(x,n)
 char x[];
 int     n;	/* the number of longs (4 bytes each) to reverse */
 { int   i1,i2,i3,i4,i;
   char  xq;
 i1=0; i2=1; i3=2; i4=3;
 for (i=0;i<n;i++) { xq=x[i1];
	 x[i1]=x[i4];    x[i4]=xq;       xq=x[i2];
	 x[i2]=x[i3];    x[i3]=xq;
	 i1+=4;  i2+=4;  i3+=4;  i4+=4;  }
 }
 /*--------------------------------------------------------------------------*/
void swapd(x,n)
 char x[];
 int     n;	/* the number of F*8 (8 bytes each) to reverse */
 { int   i1,i2,i3,i4,i5,i6,i7,i8,i;
   char  xq;
 i1=0; i2=1; i3=2; i4=3; i5=5; i6=5; i7=6; i8=7;
 for (i=0;i<n;i++) { xq=x[i1];
	 x[i1]=x[i8];    x[i8]=xq;       xq=x[i2];
	 x[i2]=x[i7];    x[i7]=xq;	 xq=x[i3];
	 x[i3]=x[i6];	 x[i6]=xq;	 xq=x[i4];
	 x[i4]=x[i5];	 x[i5]=xq;
	 i1+=8;  i2+=8;  i3+=8;  i4+=8; i5+=8; i6+=8; i7+=8; i8+=8;  }
 }
 /*--------------------------------------------------------------------------*/
int swapb(x,n)
 char x[];
 int     n;
 {
 int   i;
 char  xq;
 if (n < 2 ) { printf("error in swapb count, n = %d\n",n);  return -1; }
 if (n%2 != 0 ) n--;
 for (i=0;i<n;i += 2) { xq=x[i];  x[i] = x[i+1];  x[i+1] = xq; }
 return 1;
 }
 /*--------------------------------------------------------------------------*/
int anadecrunch(x,array,r9,nx,ny)
 /* decompress a bit stream in x; result is n I*2 elements, put in array;
	 bit slice size r9 */
 unsigned char x[];
 int nx,ny,r9;
 short array[];
 {
 short iq;
 int r0,r1,r2,r3,r4,r5,r6,r7,r8,nb,mask;
 int j,in,i,k,ix,iy;
 unsigned char xq;
 union { int i; short w; unsigned char b[4]; } y;
 /* begin execution */
 mask=1; for (i=0;i<r9;i++) mask=2*mask; mask=mask-1;
 /*printf("slice width = %d\n",r9);*/
 /*printf ("mask = %x, %d\n",mask,mask);*/
 /* determine the # of bytes to transfer to 32 bit int for fixed portion */
 if (r9 == 0) nb=0; else { if (r9 < 2 ) nb=1;
	 else { if (r9 < 10) nb=2; else nb=3;    }};
 y.i=0;
 i=0;    r1=0;   in=0;
 for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
 /* get the first value */
#if LITTLEENDIAN
 y.b[0]=x[i];  y.b[1]=x[i+1];    iq=y.w; array[in++]=iq;
#else
 y.b[0]=x[i+1];  y.b[1]=x[i];    iq=y.w; array[in++]=iq;
#endif
 /*printf("first value = %d 0x%x\n",iq, iq);*/
 r1=r1+16;
 r2=0;
 for (ix=1; ix<nx; ix++) {
 /* first the fixed slice portion */
 i=r1/8;         j=r1%8;
#if LITTLEENDIAN
 y.b[0]=x[i];
#else
 y.b[3]=x[i];
#endif
 /* test effect on timing */
#if LITTLEENDIAN
 if (nb>1) { y.b[1]=x[i+1]; if (nb>2) y.b[2]=x[i+2]; }
#else
 if (nb>1) { y.b[2]=x[i+1]; if (nb>2) y.b[1]=x[i+2]; }
#endif
 /* shift and mask out the bit slice */
 r2=(y.i>>j) & mask;
 /*printf("r2 = %x, %d\n",r2,r2);*/
 /* the variable bit portion, find the first set bit */
 r1=r1+r9;       /* bump r1 pass the fixed part */
 i=r1/8;         j=r1%8;
 if ((xq=x[i]>>j) != 0) {
 /* caught it on first byte, find the bit */
 if ((xq&1) != 0) r0=1; else {
 if ((xq&2) != 0) r0=2; else {
 if ((xq&4) != 0) r0=3; else {
 if ((xq&8) != 0) r0=4; else {
 if ((xq&16) != 0) r0=5; else {
 if ((xq&32) != 0) r0=6; else {
 if ((xq&64) != 0) r0=7; else {
 if ((xq&128) != 0) r0=8; }}}}}}}}       else {
 /* not in first byte (or part of one) checked, carry on, first count bits in
	 that first byte */
 r0=8-j;
 /* check up to 4 more bytes, if not found than an error */
 for (k=i+1;k<i+5;k++) { if ( (xq=x[k]) != 0 ) {
 /* caught it here, find the bit and then jump from loop */
 if ((xq&1) != 0) r0+=1; else {
 if ((xq&2) != 0) r0+=2; else {
 if ((xq&4) != 0) r0+=3; else {
 if ((xq&8) != 0) r0+=4; else {
 if ((xq&16) != 0) r0+=5; else {
 if ((xq&32) != 0) r0+=6; else {
 if ((xq&64) != 0) r0+=7; else {
 if ((xq&128) != 0) r0+=8; }}}}}}} break; } else { r0=r0+8; 
				 /* add 8 bits for each all zero byte */
 if (r0 > 32) { printf("DECRUNCH -- bad bit sequence, cannot continue\n");
	 printf("i = %d, r1 = %d, ix= %d, iy = %d\n",i,r1,ix,iy);
	 return -1; }       }       }       }
 r1=r1+r0;       /* update pointer */
			 /* r0 even or odd determines sign of difference */
 if ((r0&1) != 0) { 
							 /* positive case */
 r0=(r0/2)<<r9;  iq=iq+r2;       iq=iq+r0;       array[in]=iq;
  } else
 { if (r0 == 32) { 
	 /* a long one, yank out the next 17 bits and use as difference */
 i=r1/8;         j=r1%8;
#if LITTLEENDIAN
 y.b[0]=x[i];
 y.b[1]=x[i+1]; y.b[2]=x[i+2];
#else
 y.b[3]=x[i];
 y.b[2]=x[i+1]; y.b[1]=x[i+2];
#endif
 /* shift and mask out the 17 bit slice */
 r2=(y.i>>j) & 0x1ffff;
 r1=r1+17;
 /* if the top bit was set, do a sign extend, note that 32 bit arithmetic used*/
 if ( (r2& 0x10000) != 0 ) r2=r2 | 0xffff0000;
 r4=array[in-1]; r4=r4+r2; array[in]=r4; iq=r4;
 } else {
						 /* minus case (normal) */
 r0=(-r0/2)<<r9; iq=iq+r2;       iq=iq+r0;       array[in]=iq;
 }}
 in=in+1;                                }   	    /* end of ix loop */
 i=(r1+7)/8;     r1=8*i;                 }   	    /* end of iy loop */
 return 1;
 }  						     /* end of routine */
 /*--------------------------------------------------------------------------*/
int anadecrunch8(x,array,r9,nx,ny)
 /* decompress a bit stream in x; result is n I*1 elements, put in array;
	  bit slice size r9 */
 byte x[];
 int nx,ny,r9;
 byte array[];
				 /* byte version, modified from I*2 version
				 r. shine 6/5/91 */
 {
 byte iq;
 int r0,r1,r2,r3,r4,r5,r6,r7,r8,nb,mask;
 int j,in,i,k,ix,iy;
 byte xq;
 union { int i; short w; unsigned char b[4]; } y;
 /* begin execution */
 mask=1; for (i=0;i<r9;i++) mask=2*mask; mask=mask-1;
 /* determine the # of bytes to transfer to 32 bit int for fixed portion */
 if (r9 == 0) nb=0; else { if (r9 < 2 ) nb=1;
	 else { if (r9 < 10) nb=2; else nb=3;    }};
 y.i=0;
 i=0;    r1=0;   in=0;
 for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
					 /* get the first value */
 iq=x[i];        array[in++]=iq;
 r1=r1+8;
 r2=0;
 for (ix=1; ix<nx; ix++) {
					 /* first the fixed slice portion */
 i=r1/8;         j=r1%8;
#if LITTLEENDIAN
 y.b[0]=x[i];
 if (nb>1) { y.b[1]=x[i+1]; if (nb>2) y.b[2]=x[i+2]; }
#else
 y.b[3]=x[i];
 if (nb>1) { y.b[2]=x[i+1]; if (nb>2) y.b[1]=x[i+2]; }
#endif
 /* shift and mask out the bit slice */
 r2=(y.i>>j) & mask;
			 /* the variable bit portion, find the first set bit */
 r1=r1+r9;     				  /* bump r1 pass the fixed part */
 i=r1/8;         j=r1%8;
 if ((xq=x[i]>>j) != 0) {
 /* caught it on first byte, find the bit */
 if ((xq&1) != 0) r0=1; else {
 if ((xq&2) != 0) r0=2; else {
 if ((xq&4) != 0) r0=3; else {
 if ((xq&8) != 0) r0=4; else {
 if ((xq&16) != 0) r0=5; else {
 if ((xq&32) != 0) r0=6; else {
 if ((xq&64) != 0) r0=7; else {
 if ((xq&128) != 0) r0=8; }}}}}}}}       else {
 /* not in first byte (or part of one) checked, carry on, first count bits in
	 that first byte */
 r0=8-j;
 /* check up to 4 more bytes, if not found than an error */
 for (k=i+1;k<i+5;k++) { if ( (xq=x[k]) != 0 ) {
 /* caught it here, find the bit and then jump from loop */
 if ((xq&1) != 0) r0+=1; else {
 if ((xq&2) != 0) r0+=2; else {
 if ((xq&4) != 0) r0+=3; else {
 if ((xq&8) != 0) r0+=4; else {
 if ((xq&16) != 0) r0+=5; else {
 if ((xq&32) != 0) r0+=6; else {
 if ((xq&64) != 0) r0+=7; else {
 if ((xq&128) != 0) r0+=8; }}}}}}} break; } else { r0=r0+8; 
 /* add 8 bits for each all zero byte */
 if (r0 > 32) { printf("DECRUNCH -- bad bit sequence, cannot continue");
	  return -1; }       }       }       }
 r1=r1+r0;       /* update pointer */
 /* r0 even or odd determines sign of difference */
 if ((r0&1) != 0) { 
						 /* positive case */
 r0=(r0/2)<<r9;  iq=iq+r2;       iq=iq+r0;       array[in]=iq;
  } else
 { if (r0 == 32) { 
 /* a long one, yank out the next 9 bits and use as difference */
 i=r1/8;         j=r1%8;
#if LITTLEENDIAN
 y.b[0]=x[i];
 y.b[1]=x[i+1];
#else
 y.b[3]=x[i];
 y.b[2]=x[i+1];
#endif
 /* shift and mask out the 9 bit slice */
 r2=(y.i>>j) & 0x1ff;
 r1=r1+9;
 /* if the top bit was set, do a sign extend, note that 32 bit arithmetic used*/
 if ( (r2& 0x100) != 0 ) r2=r2 | 0xffffff00;
 r4=array[in-1]; r4=r4+r2; array[in]=r4; iq=r4;
 } else {
 /* minus case (normal) */
 r0=(-r0/2)<<r9; iq=iq+r2;       iq=iq+r0;       array[in]=iq;
 } }
 in=in+1;                                }       /* end of ix loop */
 i=(r1+7)/8;     r1=8*i;                 }       /* end of iy loop */
 return 1;
 }       /* end of routine */
 /*--------------------------------------------------------------------------*/
int anadecrunch32(x,array,r9,nx,ny)
 /* decompress a bit stream in x; result is n I*4 elements, put in array;
	 bit slice size r9 */
 unsigned char x[];
 int nx,ny,r9;
 int array[];
 {
 int  iq;
 int r0,r1,r2,r4,nb;
 int j,in,i,k,ix,iy, mask;
 long long	y64;
 unsigned char xq;
 union { int i; short w; unsigned char b[4]; } y;
 union { long long l64; unsigned char b[8];  } yy;
 /* begin execution */
 mask=1; for (i=0;i<r9;i++) mask=2*mask; mask=mask-1;
 /* printf("slice width = %d\n",r9); */
 /* printf ("mask = %x, %d\n",mask,mask); */
 /* determine the # of bytes to transfer to 32 bit int for fixed portion */
 nb = (r9 + 14)/8;	/* range 1 to 5 */
 if (r9 == 0) nb=0;	/* but slice = 0 a special case */
 /* printf("nb = %d\n", nb); */
 y.i=0;
 i=0;    r1=0;   in=0;
 for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
 /* get the first value, 4 bytes */
#if LITTLEENDIAN
 y.b[0]=x[i];  y.b[1]=x[i+1];  y.b[2]=x[i+2];  y.b[3]=x[i+3];
 iq = array[in++]=y.i;
#else
 y.b[0]=x[i+3];  y.b[1]=x[i+2];  y.b[2]=x[i+1];  y.b[3]=x[i];
 iq = array[in++]=y.i;
#endif
 /*printf("first value = %d 0x%x\n",iq, iq);*/
 r1=r1+32;
 r2=0;
 for (ix=1; ix<nx; ix++) {
 /* first the fixed slice portion */
 i=r1/8;         j=r1%8;
#if LITTLEENDIAN
 yy.b[0]=x[i];
 if (nb>1) { yy.b[1]=x[i+1]; if (nb>2) { yy.b[2]=x[i+2];
  if (nb>3) { yy.b[3]=x[i+3];
  if (nb>4) { yy.b[4]=x[i+4];
  }}}}
#else
 yy.b[7]=x[i];
 if (nb>1) { yy.b[6]=x[i+1]; if (nb>2) { yy.b[5]=x[i+2];
  if (nb>3) { yy.b[4]=x[i+3];
  if (nb>4) { yy.b[3]=x[i+4];
  }}}}
#endif
 /* shift and mask out the bit slice */
 r2= (int) ((yy.l64>>j) & mask);
 /*printf("r2 = %x, %d\n",r2,r2);*/
 /* the variable bit portion, find the first set bit */
 r1=r1+r9;       /* bump r1 pass the fixed part */
 i=r1/8;         j=r1%8;
 if ((xq= (x[i]>>j) ) != 0) {
 /* caught it on first byte, find the bit */
 if ((xq&1) != 0) r0=1; else {
 if ((xq&2) != 0) r0=2; else {
 if ((xq&4) != 0) r0=3; else {
 if ((xq&8) != 0) r0=4; else {
 if ((xq&16) != 0) r0=5; else {
 if ((xq&32) != 0) r0=6; else {
 if ((xq&64) != 0) r0=7; else {
 if ((xq&128) != 0) r0=8; }}}}}}}}       else {
 /* not in first byte (or part of one) checked, carry on, first count bits in
	 that first byte */
 r0=8-j;
 /* check up to 4 more bytes, if not found than an error */
 for (k=i+1;k<i+5;k++) { if ( (xq=x[k]) != 0 ) {
 /* caught it here, find the bit and then jump from loop */
 if ((xq&1) != 0) r0+=1; else {
 if ((xq&2) != 0) r0+=2; else {
 if ((xq&4) != 0) r0+=3; else {
 if ((xq&8) != 0) r0+=4; else {
 if ((xq&16) != 0) r0+=5; else {
 if ((xq&32) != 0) r0+=6; else {
 if ((xq&64) != 0) r0+=7; else {
 if ((xq&128) != 0) r0+=8; }}}}}}} break; } else { r0=r0+8; 
				 /* add 8 bits for each all zero byte */
 if (r0 > 32) { printf("DECRUNCH -- bad bit sequence, cannot continue\n");
	 printf("i = %d, r1 = %d, ix= %d, iy = %d\n",i,r1,ix,iy);
	 return -1; }       }       }       }
 r1=r1+r0;       /* update pointer */
			 /* r0 even or odd determines sign of difference */
 /*printf("r0 = %d\n", r0);*/
 if ((r0&1) != 0) { 
							 /* positive case */
 /*printf("plus case, r0, r2, iq = %d %d %d\n", r0, r2, iq);*/
 r0=(r0/2)<<r9;  iq=iq+r2;       iq=iq+r0;       array[in]=iq;
 /*printf("r0 now = %d\n", r0);*/
  } else
 { if (r0 == 32) { 
	 /* a long one, yank out the next 33 bits and use as difference */
 i=r1/8;         j=r1%8;
#if LITTLEENDIAN
 yy.b[0]=x[i];
 yy.b[1]=x[i+1]; yy.b[2]=x[i+2]; yy.b[3]=x[i+3]; yy.b[4]=x[i+4];
#else
 yy.b[7]=x[i];
 yy.b[6]=x[i+1]; yy.b[5]=x[i+2]; yy.b[4]=x[i+3]; yy.b[3]=x[i+4];
#endif
 /* shift and mask out the 33 bit slice */
 y64=(yy.l64>>j) & 0x1ffffffff;
 r1=r1+33;
 /* if the top bit was set, do a sign extend, note that 64 bit arithmetic used*/
 if ( (y64 & 0x100000000) != 0 ) y64 = y64 | 0xffffffff00000000;
 y64 = y64 + (long long) array[in-1];
 iq = array[in]= (long) y64;
 } else {
						 /* minus case (normal) */
 /*printf("minus case, r0, r2, iq = %d %d %d\n", r0, r2, iq);*/
 r0=(-r0/2)<<r9; iq=iq+r2;       iq=iq+r0;       array[in]=iq;
 }}
 in=in+1;                                }   	    /* end of ix loop */
 i=(r1+7)/8;     r1=8*i;                 }   	    /* end of iy loop */
 return 1;
 }  						     /* end of routine */
 /*--------------------------------------------------------------------------*/
int anacrunchrun(x,array,slice,nx,ny, limit)
 /* compress 16 bit array into x (a byte array) using ny blocks each of size
 nx, bit slice size slice, returns # of bytes in x */
	 unsigned char *x;
	 int nx,ny,slice, limit;
	 short array[];
 {
 void	swapl();
 struct compresshead {
 int     tsize,nblocks,bsize;
 byte    slice_size,type; } *ch;
 short iq, *p;
 unsigned nb,ixa,ixb;
 unsigned register i,j,r1;
 int	r0,r2,r3,r4,r5,r6,r7,r8,mask,fac, ii, nrun, lrun, ic;
 int	*dif, *d, nc, zq, yq, test, *dd;
 enum {RUN, LITERAL } state;
 int	i2,k,ix,iy;
 unsigned	char xq;
 union { int i; short w; unsigned char b[4]; } y;
 /* begin execution */
 limit = limit - 10;	/* some margin since we don't check all times */
 mask=1; for (i=0;i<slice;i++) mask=2*mask;
 fac=mask;       mask=mask-1; /* no inline expon. in C */
 /* determine the # of bytes to transfer to 32 bit int for fixed portion */
 if (slice == 0) nb=0; else { if (slice < 2 ) nb=1;
	 else { if (slice < 10) nb=2; else nb=3;    }};
 y.i=0;
 /* do the compression header */
 ch = (struct compresshead *) x;
 /* important note - can't use the sizeof(struct compresshead) because it
	 is 14 on some machines and rounded up to 16 on others */
 x = x + 14;
 ch->bsize = nx;  ch->nblocks = ny;  ch->slice_size = slice;  ch->type = 2;
 i=0;    r1=0;
 dif = (int *) malloc(nx*4);			/* line buffer */
 for (iy=0;iy<ny;iy++) {                 	/* start of iy (outer) loop */
 /* load the first value, reverse bytes (VAX style)*/
#if LITTLEENDIAN
 y.w=array[iy*nx]   ;x[i++]=y.b[0]    ;x[i++]=y.b[1];
#else
 y.w=array[iy*nx]   ;x[i++]=y.b[1]    ;x[i++]=y.b[0];
#endif 
  /* compute and store the first differences for this line */
 p = (array+nx*iy);	nc=nx-1;
 d=dif; yq=(int) *p++;	zq=(int) *p++;
 while (nc--) { *d++ = zq - yq; yq = zq; zq = (int) *p++; }
 r1=r1+16;
 r2=0;
 p = (array+nx*iy);	nc=nx-1;
 d=dif;
 ic = i++;			/* count position */
 r1=r1+8;		/* to cover first count */
 lrun = 0;		/* literal count) */
 while (1) {
 /* look for a run */
 /* printf("*d, *(d+1) = %d %d, nc = %d\n",*d, *(d+1), nc );*/
 y.i = *d++;
 if (nc > 1) {
 while ( y.i == *d ) {	/* at least a run of 2 */
 dd = d+1;	nrun = 2;
 while ( nc-- > 2 && y.i == *dd) {
 /* printf("run!, y.i, *dd = %d %d, nc = %d\n", y.i, *dd, nc ); */
 nrun++;  dd++; }
 /* short runs are not worth it, make the legal limit 4 */
 /* printf("nrun = %d, nc = %d\n", nrun,nc);*/
 if ( nrun >= 4 ) {	/* code the run */
 /* a previous literal ? */
 if (lrun != 0) {
 /* printf("previous was literal, ic, i = %d %d\n", ic,i);*/
 x[ic] = lrun;	i = (r1+7)/8;	lrun = 0;
 /* printf("now, i = %d\n",i );*/
 } else i=ic;
 while (nrun  > 128 )	{ /* a big one, multiple runs */
 /* need only 2 bytes to represent run, runs can't be 17 bits */
 if (nrun == 129)	/* beware the dreaded 129 */
   { x[i++] = 0x82; nrun -= 127;} else { x[i++] = 0x81;	nrun -= 128; }
#if LITTLEENDIAN
 x[i++]=y.b[0];	x[i++]=y.b[1];
#else
 x[i++]=y.b[3];	x[i++]=y.b[2];
#endif
 }
 /* printf("encoding run, nrun = %d, i=%d, iy = %d\n",nrun,i,iy); */
#if LITTLEENDIAN
 x[i++] = -(nrun-1);	x[i++]=y.b[0];	x[i++]=y.b[1];
#else
 x[i++] = -(nrun-1);	x[i++]=y.b[3];	x[i++]=y.b[2];
#endif
 /* prepare for a literal and for next run check */
 nc--;
 if (nc <= 0) goto ended_on_run;
 lrun = 0;	ic = i++;	r1 = 8*i;	d = dd;
 y.i = *d++;
 /* printf("after a run, new y.i = %d, new *d = %d\n", y.i, *d);*/
 } else { nc = nc + nrun -1;	break; }
 }	/* not a run, do next literal, assume setup for literals */
 } else if (nc <= 0) break;
 nc--;
 /* increment the literal count */
 /* printf("literal, lrun = %d, nc = %d, ic,i = %d %d\n", lrun, nc, ic,i);*/
 if (++lrun > 127)	{		/* need a new literal run count */
 x[ic] = 127;
 /* printf("ic = %d, i,r1 = %d %d\n", ic,i,r1); */
 /* bump to next byte boundary */
 i = (r1+7)/8;	ic = i++;     r1 = 8*i;		lrun = 1;
 /* printf("next ic = %d\n", ic); */
 }
 /* first the fixed slice portion */
 /* printf("y.i = %d\n", y.i);*/
 r3=(y.i>>slice);
 i=r1>>3;	/* byte number */
 j=r1 & 7;		/* bit number */
 if ( i > limit ) return -1;			/* bad news, went too far */
 /* now load nb bytes into x */
 /*low order byte of y.i is first in stream */
#if LITTLEENDIAN
 if (j == 0) {y.i=(y.i & mask); x[i]=y.b[0];}
	 else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[0];}
 if (nb>1) { x[i+1]=y.b[1]; if (nb>2) x[i+2]=y.b[2]; }
#else
 if (j == 0) {y.i=(y.i & mask); x[i]=y.b[3];}
	 else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[3];}
 if (nb>1) { x[i+1]=y.b[2]; if (nb>2) x[i+2]=y.b[1]; }
#endif 
 r1=r1+slice;       			/* bump r1 pass the fixed part */
 i=r1>>3;                j=r1 & 7;
 /* note that r3 is the # of bits required minus 1 */
 /* printf("r3 = %d\n", r3);*/
 if (r3==0) { if (j ==0 ) {x[i]=bits[j];} else {x[i]=x[i]|bits[j];}
 r1+=1;} else    {
 r3=2*r3;        if (r3<0) r3 = -r3-1;
 if (r3<31)  {
 r0=j+r3;        /* this is the bit that needs setting offset from x[i] */
 if (r0 < 8) { if (j == 0) x[i]=bits[r0]; else x[i]=x[i]|bits[r0];}
 else {if (j == 0) x[i]=0;  j=r0%8; if (r0 < 16) x[i+1]=bits[j];
 else { i2=i+r0/8; for (k=i+1;k<i2;k++) x[k]=0;  x[i2]=bits[j]; }
 }
 r1+=1+r3;
 } else {        /* big one exception, should be rare */
 /* does not need to be efficient, used rarely */
 /* printf("big one \n");*/
 if (j == 0) x[i]=0;     /* gotta zero the virgins */
 r0=j+31;        j=r0%8;         i2=i+r0/8;
 for (k=i+1;k<i2;k++) x[k]=0;    x[i2]=bits[j];
 /* recompute the difference and load 17 bits (always 3 bytes) */
 r1=r1+32;
 i=r1/8;
 j=r1%8;
 if (j == 0) x[i]=0;     /* gotta zero the virgins */
 y.i=((*(d-1))& 0x1ffff) << j;
#if LITTLEENDIAN
 x[i]=x[i] | y.b[0]; x[i+1]=y.b[1];      x[i+2]=y.b[2];
#else
 x[i]=x[i] | y.b[3]; x[i+1]=y.b[2];      x[i+2]=y.b[1];
#endif
 r1=r1+17;       } /* end of big one exception */
		 }       /* end of (r3==0) conditional */
 /* printf("end of literal, i,r1 = %d %d\n", i,r1);*/
 /* printf(" x[r1/8] = %d\n",  x[r1/8]);*/
		 }       /* end of ix loop */
 /* some checks here */
 /* bump to next byte boundary */
 /* a final literal ? */
 /* printf("final lrun = %d, ic = %d\n", lrun, ic);*/
 if (lrun != 0) { x[ic] = lrun;	lrun = 0; }
 i=(r1+7)/8;
 ended_on_run:
 r1=8*i;
		  }       /* end of iy loop */
 ch->tsize = i = i + 14;
 /* we have to put these in a form readable by the Vax (these may be used
	 by fcwrite) */
#if LITTLEENDIAN
#else
 swapl(&(ch->tsize),1); swapl(&(ch->bsize),1); swapl(&(ch->nblocks),1); 
#endif
 free(dif);
 return  i;      /*return # of bytes used */
 }       /* end of routine */
 /*--------------------------------------------------------------------------*/
int anacrunchrun8(x,array,slice,nx,ny, limit)
 /* compress 8 bit array into x (a byte array) using ny blocks each of size
 nx, bit slice size slice, returns # of bytes in x */
 unsigned char *x;
 int nx,ny,slice, limit;
 byte array[];
 {
 void	swapl();
 struct compresshead {
 int     tsize,nblocks,bsize;
 byte    slice_size,type; } *ch;
 byte	*p;
 unsigned nb,ixa,ixb;
 unsigned register i,j,r1;
 int	r0,r2,r3,r4,r5,r6,r7,r8,mask,fac, ii, nrun, lrun, ic;
 int	*dif, *d, nc, zq, yq, test, *dd;
 int	i2,k,ix,iy;
 unsigned	char xq;
 union { int i; short w; unsigned char b[4]; } y;
 /* begin execution */
 limit = limit - 10;	/* some margin since we don't check all times */
 mask=1; for (i=0;i<slice;i++) mask=2*mask;
 fac=mask;       mask=mask-1; /* no inline expon. in C */
 /* determine the # of bytes to transfer to 32 bit int for fixed portion */
 if (slice == 0) nb=0; else { if (slice < 2 ) nb=1;
	 else { if (slice < 10) nb=2; else nb=3;    }};
 y.i=0;
 /* do the compression header */
 ch = (struct compresshead *) x;
 /* important note - can't use the sizeof(struct compresshead) because it
	 is 14 on some machines and rounded up to 16 on others */
 x = x + 14;
 ch->bsize = nx;  ch->nblocks = ny;  ch->slice_size = slice;  ch->type = 3;
 i=0;    r1=0;
 dif = (int *) malloc(nx*4);			/* line buffer */
 for (iy=0;iy<ny;iy++) {                 	/* start of iy (outer) loop */
 /* load the first value */
 x[i++] = array[iy*nx];
 
 /* compute and store the first differences for this line */
 p = (array+nx*iy);	nc=nx-1;
 d=dif; yq=(int) *p++;	zq=(int) *p++;
 while (nc--) { *d++ = zq - yq; yq = zq; zq = (int) *p++; }
 r1=r1+8;
 r2=0;
 p = (array+nx*iy);	nc=nx-1;
 d=dif;
 ic = i++;			/* count position */
 r1=r1+8;		/* to cover first count */
 lrun = 0;		/* literal count) */
 while (1) {
 /* look for a run */
 /* printf("*d, *(d+1) = %d %d, nc = %d\n",*d, *(d+1), nc );*/
 y.i = *d++;
 if (nc > 1) {
 while ( y.i == *d ) {	/* at least a run of 2 */
 dd = d+1;	nrun = 2;
 while ( nc-- > 2 && y.i == *dd) {
 /* printf("run!, y.i, *dd = %d %d, nc = %d\n", y.i, *dd, nc ); */
 nrun++;  dd++; }
 /* short runs are not worth it, make the legal limit 4 */
 /* printf("nrun = %d, nc = %d\n", nrun,nc);*/
 if ( nrun >= 4 ) {	/* code the run */
 /* a previous literal ? */
 if (lrun != 0) {
 /* printf("previous was literal, ic, i = %d %d\n", ic,i);*/
 x[ic] = lrun;	i = (r1+7)/8;	lrun = 0;
 /* printf("now, i = %d\n",i );*/
 } else i=ic;
 while (nrun  > 128 )	{ /* a big one, multiple runs */
 /* printf("big run, nrun = %d\n", nrun); */
 /* need only 2 bytes to represent run, runs can't be 17 bits */
 if (nrun == 129)	/* beware the dreaded 129 */
   { x[i++] = 0x82; nrun -= 127;} else { x[i++] = 0x81;	nrun -= 128; }
#if LITTLEENDIAN
 x[i++]=y.b[0];	x[i++]=y.b[1];
#else
 x[i++]=y.b[3];	x[i++]=y.b[2];
#endif
 }
 /* printf("encoding run, nrun = %d, i=%d, iy = %d\n",nrun,i,iy); */
#if LITTLEENDIAN
 x[i++] = -(nrun-1);	x[i++]=y.b[0];	x[i++]=y.b[1];
#else
 x[i++] = -(nrun-1);	x[i++]=y.b[3];	x[i++]=y.b[2];
#endif
 /* prepare for a literal and for next run check */
 nc--;
 if (nc <= 0) goto ended_on_run;
 lrun = 0;	ic = i++;	r1 = 8*i;	d = dd;
 y.i = *d++;
 /* printf("after a run, new y.i = %d, new *d = %d\n", y.i, *d);*/
 } else { nc = nc + nrun -1;	break; }
 }	/* not a run, do next literal, assume setup for literals */
 } else if (nc <= 0) break;
 nc--;
 /* increment the literal count */
 /* printf("literal, lrun = %d, nc = %d, ic,i = %d %d\n", lrun, nc, ic,i);*/
 if (++lrun > 127)	{		/* need a new literal run count */
 x[ic] = 127;
 /* printf("ic = %d, i,r1 = %d %d\n", ic,i,r1); */
 /* bump to next byte boundary */
 i = (r1+7)/8;	ic = i++;     r1 = 8*i;		lrun = 1;
 /* printf("next ic = %d\n", ic); */
 }
 /* first the fixed slice portion */
 /* printf("y.i = %d\n", y.i);*/
 r3=(y.i>>slice);
 i=r1>>3;	/* byte number */
 j=r1 & 7;		/* bit number */
 if ( i > limit ) return -1;			/* bad news, went too far */
 /* now load nb bytes into x */
 /*low order byte of y.i is first in stream */
#if LITTLEENDIAN
 if (j == 0) {y.i=(y.i & mask); x[i]=y.b[0];}
	 else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[0];}
 if (nb>1) { x[i+1]=y.b[1]; if (nb>2) x[i+2]=y.b[2]; }
#else
 if (j == 0) {y.i=(y.i & mask); x[i]=y.b[3];}
	 else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[3];}
 if (nb>1) { x[i+1]=y.b[2]; if (nb>2) x[i+2]=y.b[1]; }
#endif
 
 r1=r1+slice;       			/* bump r1 pass the fixed part */
 i=r1>>3;                j=r1 & 7;
 /* note that r3 is the # of bits required minus 1 */
 if (r3==0) { if (j ==0 ) {x[i]=bits[j];} else {x[i]=x[i]|bits[j];}
 r1+=1;} else    {
 r3=2*r3;        if (r3<0) r3 = -r3-1;
 if (r3<31)  {
 r0=j+r3;        /* this is the bit that needs setting offset from x[i] */
 if (r0 < 8) { if (j == 0) x[i]=bits[r0]; else x[i]=x[i]|bits[r0];}
 else {if (j == 0) x[i]=0;  j=r0%8; if (r0 < 16) x[i+1]=bits[j];
 else { i2=i+r0/8; for (k=i+1;k<i2;k++) x[k]=0;  x[i2]=bits[j]; }
 }
 r1+=1+r3;
 } else {        /* big one exception, should be rare */
 /* does not need to be efficient, used rarely */
 /* printf("big one \n");*/
 if (j == 0) x[i]=0;     /* gotta zero the virgins */
 r0=j+31;        j=r0%8;         i2=i+r0/8;
 for (k=i+1;k<i2;k++) x[k]=0;    x[i2]=bits[j];
 /* recompute the difference and load 9 bits (always 2 bytes) */
 r1=r1+32;
 i=r1/8;
 j=r1%8;
 if (j == 0) x[i]=0;     /* gotta zero the virgins */
 y.i=((*(d-1))& 0x1ff) << j;
#if LITTLEENDIAN
 x[i]=x[i] | y.b[0]; x[i+1]=y.b[1];
#else
 x[i]=x[i] | y.b[3]; x[i+1]=y.b[2];
#endif
 r1=r1+9;       } /* end of big one exception */
		 }       /* end of (r3==0) conditional */
 /* printf("end of literal, i,r1 = %d %d\n", i,r1);*/
 /* printf(" x[r1/8] = %d\n",  x[r1/8]);*/
		 }       /* end of ix loop */
 /* some checks here */
 /* bump to next byte boundary */
 /* a final literal ? */
 /* printf("final lrun = %d, ic = %d\n", lrun, ic);*/
 if (lrun != 0) { x[ic] = lrun;	lrun = 0; }
 i=(r1+7)/8;
 ended_on_run:
 r1=8*i;
		  }       /* end of iy loop */
 ch->tsize = i = i + 14;
 /* we have to put these in a form readable by the Vax (these may be used
	 by fcwrite) */
#if LITTLEENDIAN
#else
 swapl(&(ch->tsize),1); swapl(&(ch->bsize),1); swapl(&(ch->nblocks),1); 
#endif
 free(dif);
 return  i;      /*return # of bytes used */
 }
 /*--------------------------------------------------------------------------*/
int anadecrunchrun(x,array,r9,nx,ny)
 /* decompress a bit stream in x; result is n I*2 elements, put in array;
	 bit slice size r9 */
 /* this version handles the run length encoding used in anacrunchrun */
 unsigned char x[];
 int nx,ny,r9;
 short array[];
 {
 short iq;
 int r0,r1,r2,r3,r4,r5,r6,r7,r8,nb,mask,nrun,n,nc;
 int j,in,i,k,ix,iy;
 unsigned char xq;
 union { int i; short w; unsigned char b[4]; } y;
 /* begin execution */
 mask=1; for (i=0;i<r9;i++) mask=2*mask; mask=mask-1;
 /*printf("slice width = %d\n",r9);*/
 /*printf ("mask = %x, %d\n",mask,mask);*/
 /* determine the # of bytes to transfer to 32 bit int for fixed portion */
 if (r9 == 0) nb=0; else { if (r9 < 2 ) nb=1;
	 else { if (r9 < 10) nb=2; else nb=3;    }};
 y.i=0;
 i=0;    r1=0;   in=0;
 for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
 /* get the first value, assume bytes reversed */
#if LITTLEENDIAN
 y.b[0]=x[i++];	y.b[1]=x[i++];    iq=y.w; array[in++]=iq;
#else
 y.b[1]=x[i++];	y.b[0]=x[i++];    iq=y.w; array[in++]=iq;
#endif
 /* printf("first value = %d 0x%x\n",iq, iq); */
 r1=r1+16;
 r2=0;
 nc=nx-1;
 while (nc>0) {
 /* look at the next run length code */
 /* printf("i = %d\n", i); */
 nrun = (int) x[i++];
 /* printf("nrun = %d\n", nrun); */
 if (nrun > 127) {	/* a run of a constant difference */
 n = 255 - nrun + 2;
 nc = nc - n;
#if LITTLEENDIAN
 y.b[0]=x[i++];	y.b[1]=x[i++];
#else
 y.b[1]=x[i++];	y.b[0]=x[i++];
#endif
 /* printf("increment (run) = %d\n", y.w); */
 while (n--) {
 array[in] = array[in-1] + y.w;	in++;
 }
 iq = array[in-1];	r1=8*i;
 } else {	/* a literal */
 r1 = 8 * i;
 nc = nc - nrun;
 while(nrun--) {
 /* first the fixed slice portion */
 i=r1/8;         j=r1%8;
 /* printf("start literal, i,r1 = %d %d\n", i,r1); */
#if LITTLEENDIAN
 y.b[0]=x[i];
 /* test effect on timing */
 if (nb>1) { y.b[1]=x[i+1]; if (nb>2) y.b[2]=x[i+2]; }
#else
 y.b[3]=x[i];
 /* test effect on timing */
 if (nb>1) { y.b[2]=x[i+1]; if (nb>2) y.b[1]=x[i+2]; }
#endif
 /* shift and mask out the bit slice */
 r2=(y.i>>j) & mask;
 /* printf("r2 = %x, %d\n",r2,r2);*/
 /* the variable bit portion, find the first set bit */
 r1=r1+r9;       /* bump r1 pass the fixed part */
 i=r1/8;         j=r1%8;
 if ((xq=x[i]>>j) != 0) {
 /* caught it on first byte, find the bit */
 if ((xq&1) != 0) r0=1; else {
 if ((xq&2) != 0) r0=2; else {
 if ((xq&4) != 0) r0=3; else {
 if ((xq&8) != 0) r0=4; else {
 if ((xq&16) != 0) r0=5; else {
 if ((xq&32) != 0) r0=6; else {
 if ((xq&64) != 0) r0=7; else {
 if ((xq&128) != 0) r0=8; }}}}}}}}       else {
 /* not in first byte (or part of one) checked, carry on, first count bits in
	 that first byte */
 r0=8-j;
 /* check up to 4 more bytes, if not found than an error */
 for (k=i+1;k<i+5;k++) { if ( (xq=x[k]) != 0 ) {
 /* caught it here, find the bit and then jump from loop */
 if ((xq&1) != 0) r0+=1; else {
 if ((xq&2) != 0) r0+=2; else {
 if ((xq&4) != 0) r0+=3; else {
 if ((xq&8) != 0) r0+=4; else {
 if ((xq&16) != 0) r0+=5; else {
 if ((xq&32) != 0) r0+=6; else {
 if ((xq&64) != 0) r0+=7; else {
 if ((xq&128) != 0) r0+=8; }}}}}}} break; } else { r0=r0+8; 
				 /* add 8 bits for each all zero byte */
 if (r0 > 32) { printf("DECRUNCH -- bad bit sequence, cannot continue\n");
	 printf("i = %d, r1 = %d, ix= %d, iy = %d\n",i,r1,ix,iy);
	 return -1; }       }       }       }
 r1=r1+r0;       /* update pointer */
			 /* r0 even or odd determines sign of difference */
 if ((r0&1) != 0) { 
							 /* positive case */
 r0=(r0/2)<<r9;  iq=iq+r2;       iq=iq+r0;       array[in]=iq;
 /* printf("r0,r2,iq = %d %d %d\n", r0,r2,iq);*/
  } else
 { if (r0 == 32) { 
	 /* a long one, yank out the next 17 bits and use as difference */
 i=r1/8;         j=r1%8;
#if LITTLEENDIAN
 y.b[0]=x[i];
 y.b[1]=x[i+1]; y.b[2]=x[i+2];
#else
 y.b[3]=x[i];
 y.b[2]=x[i+1]; y.b[1]=x[i+2];
#endif
 /* shift and mask out the 17 bit slice */
 r2=(y.i>>j) & 0x1ffff;
 r1=r1+17;
 /* if the top bit was set, do a sign extend, note that 32 bit arithmetic used*/
 if ( (r2& 0x10000) != 0 ) r2=r2 | 0xffff0000;
 /* printf("big one, r2 = %d, array[in-1] = %d\n", r2, array[in-1]);*/
 r4=array[in-1]; r4=r4+r2; array[in]=r4; iq=r4;
 } else {
						 /* minus case (normal) */
 r0=(-r0/2)<<r9; iq=iq+r2;       iq=iq+r0;       array[in]=iq;
 /* printf("r0,r2,iq = %d %d %d\n", r0,r2,iq); */
 }}
 /* printf("end literal, i,r1 = %d %d\n", i,r1); */
 in=in+1;
 }
 /* printf("r1, i after nrun exhausted = %d %d\n",r1,i); */
 i=(r1+7)/8;     r1=8*i;
 /* printf("new i = %d\n",i); */
 }
 }   	    /* end of ix loop */
 if (nc < 0) {
  printf("bad loop in decrunchrun, nc=%d, iy=%d, in= %d\n",nc,iy,in);  return -1; }
 
 i=(r1+7)/8;     r1=8*i;                 }   	    /* end of iy loop */
 return 1;
 }  						     /* end of routine */
 /*--------------------------------------------------------------------------*/
int anadecrunchrun8(x,array,r9,nx,ny)
 /* decompress a bit stream in x; result is n I*2 elements, put in array;
	 bit slice size r9 */
 /* this version handles the run length encoding used in anacrunchrun */
 byte x[];
 int nx,ny,r9;
 byte array[];
 {
 byte iq;
 int r0,r1,r2,r3,r4,r5,r6,r7,r8,nb,mask,nrun,n,nc;
 int j,in,i,k,ix,iy;
 unsigned char xq;
 union { int i; short w; unsigned char b[4]; } y;
 /* begin execution */
 mask=1; for (i=0;i<r9;i++) mask=2*mask; mask=mask-1;
 /* determine the # of bytes to transfer to 32 bit int for fixed portion */
 if (r9 == 0) nb=0; else { if (r9 < 2 ) nb=1;
	 else { if (r9 < 10) nb=2; else nb=3;    }};
 y.i=0;
 i=0;    r1=0;   in=0;
 for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
 /* get the first value */
 iq=x[i++];	array[in++]=iq;
 /* printf("first value = %d 0x%x\n",iq, iq); */
 r1=r1+16;
 r2=0;
 nc=nx-1;
 /* printf("nc = %d\n", nc); */
 while (nc>0) {
 /* look at the next run length code */
 /* printf("i = %d\n", i); */
 nrun = (int) x[i++];
 /* printf("nrun = %d\n", nrun); */
 if (nrun > 127) {	/* a run of a constant difference */
 n = 255 - nrun + 2;
 nc = nc - n;
#if LITTLEENDIAN
 y.b[0]=x[i++];	y.b[1]=x[i++];
#else
 y.b[1]=x[i++];	y.b[0]=x[i++];
#endif
 /* printf("increment (run) = %d of length %d\n", y.w,n); */
 while (n--) {
 array[in] = array[in-1] + y.w;	in++;
 }
 iq = array[in-1];	r1=8*i;
 } else {	/* a literal */
 r1 = 8 * i;
 nc = nc - nrun;
 while(nrun--) {
 /* first the fixed slice portion */
 i=r1/8;         j=r1%8;
 /* printf("start literal, i,r1 = %d %d\n", i,r1); */
#if LITTLEENDIAN
 y.b[0]=x[i];
 /* test effect on timing */
 if (nb>1) { y.b[1]=x[i+1]; if (nb>2) y.b[2]=x[i+2]; }
#else
 y.b[3]=x[i];
 /* test effect on timing */
 if (nb>1) { y.b[2]=x[i+1]; if (nb>2) y.b[1]=x[i+2]; }
#endif
 /* shift and mask out the bit slice */
 r2=(y.i>>j) & mask;
 /* the variable bit portion, find the first set bit */
 r1=r1+r9;       /* bump r1 pass the fixed part */
 i=r1/8;         j=r1%8;
 if ((xq=x[i]>>j) != 0) {
 /* caught it on first byte, find the bit */
 if ((xq&1) != 0) r0=1; else {
 if ((xq&2) != 0) r0=2; else {
 if ((xq&4) != 0) r0=3; else {
 if ((xq&8) != 0) r0=4; else {
 if ((xq&16) != 0) r0=5; else {
 if ((xq&32) != 0) r0=6; else {
 if ((xq&64) != 0) r0=7; else {
 if ((xq&128) != 0) r0=8; }}}}}}}}       else {
 /* not in first byte (or part of one) checked, carry on, first count bits in
	 that first byte */
 r0=8-j;
 /* check up to 4 more bytes, if not found than an error */
 for (k=i+1;k<i+5;k++) { if ( (xq=x[k]) != 0 ) {
 /* caught it here, find the bit and then jump from loop */
 if ((xq&1) != 0) r0+=1; else {
 if ((xq&2) != 0) r0+=2; else {
 if ((xq&4) != 0) r0+=3; else {
 if ((xq&8) != 0) r0+=4; else {
 if ((xq&16) != 0) r0+=5; else {
 if ((xq&32) != 0) r0+=6; else {
 if ((xq&64) != 0) r0+=7; else {
 if ((xq&128) != 0) r0+=8; }}}}}}} break; } else { r0=r0+8; 
				 /* add 8 bits for each all zero byte */
 if (r0 > 32) { printf("DECRUNCH -- bad bit sequence, cannot continue\n");
	 printf("i = %d, r1 = %d, ix= %d, iy = %d\n",i,r1,ix,iy);
	 return -1; }       }       }       }
 r1=r1+r0;       /* update pointer */
			 /* r0 even or odd determines sign of difference */
 if ((r0&1) != 0) { 
							 /* positive case */
 r0=(r0/2)<<r9;  iq=iq+r2;       iq=iq+r0;       array[in]=iq;
 /* printf("r0,r2,iq = %d %d %d\n", r0,r2,iq);*/
  } else
 { if (r0 == 32) { 
	 /* a long one, yank out the next 9 bits and use as difference */
 i=r1/8;         j=r1%8;
#if LITTLEENDIAN
 y.b[0]=x[i];
 y.b[1]=x[i+1];
#else
 y.b[3]=x[i];
 y.b[2]=x[i+1];
#endif
 /* shift and mask out the 9 bit slice */
 r2=(y.i>>j) & 0x1ff;
 r1=r1+9;
 /* if the top bit was set, do a sign extend, note that 32 bit arithmetic used*/
 if ( (r2& 0x100) != 0 ) r2=r2 | 0xffffff00;
 /* printf("long one decoded, r2 = %d, array[in-1]=%d\n", r2, array[in-1]); */
 r4=array[in-1]; r4=r4+r2; array[in]=r4; iq=r4;
 } else {
						 /* minus case (normal) */
 r0=(-r0/2)<<r9; iq=iq+r2;       iq=iq+r0;       array[in]=iq;
 /* printf("r0,r2,iq = %d %d %d\n", r0,r2,iq); */
 }}
 /* printf("end literal, i,r1 = %d %d\n", i,r1); */
 in=in+1;
 }
 /* printf("r1, i after literal nrun exhausted = %d %d\n",r1,i); */
 i=(r1+7)/8;     r1=8*i;
 /* printf("new i = %d\n",i); */
 }
 }   	    /* end of ix loop */
 if (nc < 0) {
  printf("bad loop in decrunchrun8, nc=%d, iy=%d, in= %d\n",nc,iy,in);
  return -1; }
 
 i=(r1+7)/8;     r1=8*i;                 }   	    /* end of iy loop */
 return 1;
  }  						     /* end of routine */
 /*--------------------------------------------------------------------------*/