acms.c

ACMS压缩算法

/*  Arithmetic coding routines by A.C. Popat, 1990.  For theory, see
*      Popat, Ashok C., "Scalar Quantization with Arithmetic Coding,"
*      SM thesis, MIT Dept. of Elec. Eng. and Comp. Sci., Cambridge,
*      Massaachusetts, June 1990.
*/


/*
*  Copyright 1990 Massachusetts Institute of Technology
*  All rights reserved.
*  
*  Permission to use, copy, or modify this software and its documentation
*  for educational and research purposes only and without fee is hereby
*  granted, provided that this copyright notice appear on all copies and
*  supporting documentation.  For any other uses of this software, in
*  original or modified form, including but not limited to distribution
*  in whole or in part, specific prior permission must be obtained from
*  M.I.T. and the author.  These programs shall not be used, rewritten,
*  or adapted as the basis of a commercial software or hardware product
*  without first obtaining appropriate licenses from M.I.T.  M.I.T. makes
*  no representations about the suitability of this software for any
*  purpose.  It is provided "as is" without express or implied warranty.
*/



#include <stdio.h>
#include <math.h>

#define MAXK 512        /* max. allowable size of source alphabet */
#define MAXLEN 10000    /* max. allowable length of source sequence */
#define MAXCBITS 100000 /* max. allowable length of code bit sequence */
#define MAXPR 16        /* max. allowable precision */
#define MAXCCW 12       /* max. allowable width of carry-control field */

typedef unsigned long codereg;

static int k;       /* size of source alphabet */
static int pr;      /* precision of arithmetic operations (e.g., 16) */
static int ccw;     /* width of carry-control field  (e.g., 8) */
static codereg twotothe[32];  /* array of masks that isolate individual bits */
static double *p;           /* source probabilities */
static codereg q[MAXK];     /* rounded source probabilities */
static codereg cumq[MAXK];  /* cumulative rounded source probabilities */
static int verbose = 0;      /* set to zero to make silent */

int optprnd(); /* routine to round probs to admissible vals */
static void printulong();    /* debugging utility */
static int codereglen();     /* returns length, in bits, of type codereg */
static int bad_pr_ccw();     /* checks for validity of pr and ccw */
static void inittwotothe();  /* sets twotothe[i] = 2^i */
static double entropy();     /* returns entropy of source (not necessary
				for encoding and decoding) */

/* Macro for testing individual bits:
** synopsis: bitvalue = bittst(register,bitnumber); */
#define bittst(r,n) (((r)&twotothe[n]) > 0 ? 1:0)

/* Macro for setting value of individual bits:
** synopsis: bitset(register,bitnumber,bitvalue); */
#define bitset(r,n,b) if (b>0) r |= twotothe[n]; else r &= ~twotothe[n]

/***************************************************************************/
/* Macro for counting the number of consecutive ones in the high end of
** the carry-control field.  Needed by "carryalert." Result is in "nco": */
#define countnco nco=0; while( bittst(ex,ccfmsb-nco) ) nco++

/***************************************************************************/
/* Macro for putting a bit in the code string:
** synopsis: sendbit(bitvalue); */

#define sendbit(b) codebits[(*nbits)++] = b;\
                                                          \
                   /* if 1 and if there is no carry trap currently in ex,\
                   ** then update ncojt.  Otherwise, zero ncojt: */\
                                                          \
                   if ((b==1)&&(ctbptr>ccfmsb)) ncojt++;\
                   else ncojt = 0

/***************************************************************************/
/* Macro for conditionally inserting a carry-trap bit. */

#define carryalert countnco;\
                   if ((nco+ncojt) >= ccw) {\
                     /* carry trap needs to be set */\
                                                               \
		     nco = ccw - ncojt; \
                                                               \
                     /* establish position of carry-trap bit */\
                                                               \
		     ctbptr = ccfmsb + ncojt - ccw + 1;\
                                                               \
		     /* if carry trap bit is outside ex, then \
                     ** it must be in position ccfmsb+1, so   \
		     ** simply send zero: */\
                                                               \
                     if (ctbptr>ccfmsb) {sendbit(0);}\
                                                               \
		     /* othersize, zero the trap bit and send \
                     ** the most significant bit in ccf, which \
                     ** must be a 1: */\
                                                               \
                     else { \
                       bitset(ex,ctbptr,0);\
      		       {sendbit(1);}\
		     } \
		     }


/***************************************************************************/
/* Macro for shifting ex left, emitting the shifted bit, and conditionally
** setting carry traps. */

#define shiftexout \
      sendbit(bittst(ex,ccfmsb));\
      ex += ex;         /* shift ex left by one */ \
      ex &= efullmask;  /* keep unused part of register clean  */ \
      if (ctbptr <= ccfmsb) ctbptr++; /* if carry trap in ex, shift it */\
                                                                         \
      /* if carry trap now outside ex, then conditionally set new trap: */\
      if (ctbptr > ccfmsb) {carryalert;}\


/***************************************************************************/
void encode (alphabetsize, probs, letters, nletters, precision,
	     ccwidth, codebits, nbits)
     int alphabetsize, precision, ccwidth, letters[], codebits[];
     double probs[];
     long nletters, *nbits;
{
  /* Encoder registers: */
  codereg ex,ew,ed;

  /* Encoder masks:  
  **         efullmask isolates the low-order (pr+ccw) bits of a register */
  codereg efullmask;

  int status;      /* status holds ret'd val from optprnd */
  long n;          /* n is the source letter index */
  int let;         /* let holds current letter */

  /* the following variables are used in handling carries */
  int ncojt;       /* number of consecutive ones just transmitted */
  int nco;         /* number of consecutive ones in high end of ccf. */
  int ctbptr;      /* position of carry-trap bit */
  int ccfmsb;      /* position of most signif. bit in carry-control field */

  int i, roundaccum; /* debug */

  /* set file-local variables: */
  p = &probs[0];
  k = alphabetsize;
  pr = precision;
  ccw = ccwidth;

  if(bad_pr_ccw()) {
    fprintf(stderr,"pr and/or ccw have bad values in encode\n");
    exit(-1);
  }

  /* Initialize masks */
  inittwotothe();
  efullmask = twotothe[pr+ccw] - 1;
  ccfmsb = pr + ccw - 1;
  
  /* Round probabilities to admissible values: */
  status = optprnd();
  if (status!=0) {
    fprintf(stderr,"optprnd failed in encode\n");
    exit(-1);
  }


  /* set ex=0 and ew = 0.111111...: */
  ex = 0;
  ew = twotothe[pr] - 1;

  *nbits = 0;          /* no code bits generated yet */
  ncojt = 0;           /* no consecutive ones just transmitted */
  nco = 0;             /* no ones yet in the high end of ccf */
  ctbptr = ccfmsb + 1; /* initially, the carry-trap bit is not in ccf */

  /** main encoding loop: */
  for (n=0;n<nletters;n++) {

    let = letters[n];     /* current source letter */

    /* add trunc(ew*cumq[let]) to ex: */
    ex += (ew*cumq[let])>>pr;
    
    /* set ed = truncated product of ew and q[let]: */
    ed = ew*q[let];
    ed = (ed>>pr) << pr;     /* truncate */
    if (ed==0) {
      fprintf(stderr,"subinterval width is zero in encode\n");
      exit(-1);
    }

    /* normalize ed, and in the process, emit code bits: */
    while (bittst(ed,2*pr-1) == 0) {
      shiftexout;
      ed += ed;         /* shift ed left by one for next iteration */
    }
    /* truncate ed to pr bits to get new ew.  since the radix point of ed
    ** is to the left of bit (2*pr-1) and the radix point of ew is
    ** to the left of bit pr-1, obtain ew by shifting ed right by pr bits: */
    ew = ed >> pr;
    /* end of the encoding operation for the nth letter */
  }

  /* conclude by flushing out the code register ex. */
  for (n=0;n<(pr+ccw);n++) {
    shiftexout;
  }    
}  

#define shiftdxin \
      dx += dx;         /* shift dx left by one */ \
      if (m<nbits) dx += codebits[m++]; \
      if (bwacptr < ccw) bwacptr++; \
      if (bwacptr >= ccw) { \
         if ((dx&dccmask) == dccmask) { \
            dx += codebits[m++];\
            bwacptr = 0; \
	 }\
      }

void decode (alphabetsize, probs, letters, nletters, precision,
	     ccwidth, codebits, nbits)
     int alphabetsize, precision, ccwidth, letters[], codebits[];
     double probs[];
     long nletters, nbits;
{
  codereg dx,dw,dd,dy,dz;

  /* ADD PROCESSING STEP COMMENTS HERE FOR DECODE */

  /* decoder masks:  
  **         dccmask isolates the carry-control field (bits o thru ccw-1)
  **         (these masks are initialized in "init()" ):  */
  codereg dccmask;

  int status;   /* status holds ret'd val from init */
  long n;          /* n is the source letter index */
  int let;         /* let holds current letter */
  int m;        /* number of code bits used so far */
  int l;        /* loop index for candidate letters */
  int a,b;      /* indices for binary search */

  /* the following variables are used in checking for processing carries: */
  int bwacptr;  /* position of bit with carry added to it */

  /* set file-local variables: */
  p = &probs[0];
  k = alphabetsize;
  pr = precision;
  ccw = ccwidth;

  if(bad_pr_ccw()) {
    fprintf(stderr,"pr and/or ccw have bad values in decode\n");
    exit(-1);
  }

  /* Initialize masks */
  inittwotothe();
  dccmask = twotothe[ccw] - 1;

  /* Round probabilities: */
  status = optprnd();
  if (status!=0) {
    fprintf(stderr,"optprnd failed in decode\n");
    exit(-1);
  }

  /* set dx=0 and dw = 0.111111...: */
  dx = 0;
  dw = twotothe[pr] - 1;

  m = 0;          /* no code bits have been processed yet */
  bwacptr = ccw;  /* bit with added carry is not initially in ccf */

  /* shift in pr+2*ccw code bits: */
  for (n=0;n<(pr+ccw+ccw);n++) {
    shiftdxin;
  }

  /** main decoding loop: */
  for (n=0;n<nletters;n++) {

    /* decode current letter by determining what was added to ex. */
    dy = dx >> ccw;       /* get arith field of dx in dy */

    /* simple binary search... */
    a = 0;
    b = k-1;
    if (dz = (cumq[b]*dw) >> pr , dy < dz) {
      while ( l = (a + b) >> 1, l > a) {
	if (dz = (cumq[l]*dw) >> pr , dy < dz) b = l;
	else a = l;
      }
      let = a;
    }
    else let = b;

    letters[n] = let;
    dx -= ( (cumq[let]*dw) >> pr) << ccw;

    /* set dd = double precision product of dw and q[let]: */
    dd = dw*q[let];

    dd = (dd>>pr) << pr;     /* quantize */
    if (dd==0) {
      fprintf(stderr,"subinterval width is zero in decode\n");
      exit(-1);
    }
    /* normalize dd, and in the process, shift in code bits: */

    while (bittst(dd,2*pr-1) == 0) {
      shiftdxin;
      dd += dd;         /* shift dd left by one */
    }

    dw = dd >> pr;

    /* end of the decoding operation for the nth letter */
  }
  
}  

#define BIGNUM 1.0e20     /* huge value to assign "infinite" cost */

extern int optprnd()
{
  double benefit[MAXK], maxbenefit, penalty[MAXK], minpenalty,
                  cost[MAXK], mincost;
  codereg sumq;
  int i,j,l;

#ifdef OPT_RND
  for (i=0;i<k;i++) q[i] = 2;
#else
  for (i=0;i<k;i++) {
    q[i] = twotothe[pr]*p[i];
    if (q[i] < 2) q[i] = 2;
  }
#endif  

  sumq=0; for(i=0;i<k;i++) sumq += q[i];

  /* calculate relative penalty of reducing q[i] by 1: */
  for (i=0;i<k;i++)
    penalty[i] = -p[i] * log(((double) q[i]-1)/((double) q[i]));
     
  /* keep subtracting 1's until sumq = twotothe[pr]: */
  while (sumq > twotothe[pr]) {
    
    /* find j = arg min (penalty[]): */
    minpenalty = BIGNUM;
    for (i=0;i<k;i++) {
      if ((penalty[i] < minpenalty) && (q[i] > 2)) {
	minpenalty = penalty[i];
	j = i;
      }
    }
    /* decrease q[j] by 1 and compute penalty of decreasing it further: */
    q[j] -= 1;
    sumq -= 1;
    penalty[j] = -p[j] * log(((double) q[j]-1)/((double) q[j]));
  }
  

  /* calculate relative benefit of increasing q[i] by 1: */
  for (i=0;i<k;i++) benefit[i] = p[i] * log(((double) q[i]+1)/((double) q[i]));

  /* keep adding 1's until sumq = twotothe[pr]: */
  while (sumq < twotothe[pr]) {
    
    /* find j = arg max(benefit[]): */
    maxbenefit = 0.0;
    for (i=0;i<k;i++) {
      if (benefit[i] > maxbenefit) {
	maxbenefit = benefit[i];
	j = i;
      }
    }
    /* increase q[j] by 1 and compute benefit of increasing it further: */
    q[j] += 1;
    sumq += 1;
    benefit[j] = p[j] * log(((double) q[j]+1)/((double) q[j]));
  }
  cumq[0] = 0;
  for (i=1;i<k;i++) cumq[i] = cumq[i-1] + q[i-1];
  return(0);
}
#undef BIGNUM    

/************************* some useful routines ****************************/

static void printulong(str,r)   /* prints out a register as 1's and 0's */
     codereg r;
     char *str;
{
  int i;
  if (verbose) {
    printf("%s ",str);
    for (i=31;i>=0;i--) {
      printf("%d",bittst(r,i));
      if ((i==pr)||(i==ccw)) printf(" ");
    }
    printf("    %9.6lf",((double) r)/((double) twotothe[pr]));
    printf("\n");
  }
}

static int bad_pr_ccw()       /* checks validity of pr and ccw */
{
  int rw;
  rw = codereglen();
  if (pr > rw/2) {
    fprintf(stderr,"pr must not have more bits than half a register\n");
    exit(-1);
  }
  if (pr+ccw > rw) {
    fprintf(stderr,"sum of pr and ccw too high in bad_pr_ccw:");
    fprintf(stderr,"%d + %d >  %d (register width)\n",pr,ccw,rw);
    return(-1);
  }
  else return(0);
}

static void inittwotothe()    /* initialize array twotothe[] */
{
  int i;
  twotothe[0] = 1;
  for (i=1;i<32;i++) twotothe[i] = 2*twotothe[i-1];    /* = 2^i */
}
  
static int codereglen()   /* determines the length, in bits, of type codereg */
{
  codereg r = ~0;                /* initialize r to all ones */
  int len = 0;
  while (r <<= 1) len++;         /* see how many shifts are require to clear */
  return(len+1);
}

#define SMALLNUM 1.0e-20
static double entropy()     /* returns entropy of source */
{
  int i;
  double sum=0.0;
  for (i=0;i<k;i++) if (p[i] > SMALLNUM) sum -= p[i]*log(p[i]);
  return(sum/log(2.0));
}
#undef SMALLNUM