stats.c

source code for arithmatic coding

    arithmetic_encode(low_w, high_w, pContext->total);

    if (symbol != 1)        /* If not the special ESC / NOT_KNOWN symbol */
    {
            /*update the singleton count if symbol was previously a singleton */
        if (pContext->type == DYNAMIC && high-low == pContext->incr)
            pContext->nSingletons -= pContext->incr;

            /* increment the symbol's frequency count */
        INCR_SYMBOL_PROB(pContext, symbol, low, high, pContext->incr);
    }

    adjust_zero_freq(pContext);

    while (pContext->total > Max_frequency)
        halve_context(pContext);

    if (symbol == 1) return NOT_KNOWN;
    return 0;
}/* encode() */




/*
 *
 * decode function is passed a context, and returns a symbol
 *
 */
int 
decode(context *pContext)
{
    int    symbol;
    int p, m, e;
    freq_value low, high, target;

    int            n = pContext->max_length;
    freq_value total = pContext->total;
    freq_value    *M = pContext->tree;
 
    target = arithmetic_decode_target(total);

#ifdef MOST_PROB_AT_END
	/* Check if most probable symbol (shortcut decode)
	 */
    if (target >= total - pContext->most_freq_count) {
	  arithmetic_decode( total - pContext->most_freq_count, total, total);
	  symbol = pContext->most_freq_symbol;
	  low    = pContext->most_freq_pos;
	  high   = low + pContext->most_freq_count;

      INCR_SYMBOL_PROB(pContext, symbol, low, high, pContext->incr);

      if (symbol != 1) 
        if (pContext->type == DYNAMIC && high-low == pContext->incr)
            pContext->nSingletons -= pContext->incr;
	}
    else
	/* Not MPS, have to decode slowly */
  {
	if (target >= pContext->most_freq_pos)
	    target += pContext->most_freq_count;
#endif
    p = 1; low = 0;
    while ( ((p << 1) <= pContext->max_length ) && (M[p] <= target)) {
        target -= M[p];
        low    += M[p];
        p      <<= 1;
    }

    symbol = p;
    m = p >> 1;
    e = 0;

    while (m >= 1) {
        if (symbol + m <= n) {
            e += M[symbol + m];
            if (M[symbol] - e <= target) {
                target    -= M[symbol] - e;
                low       += M[symbol] - e;
                if (symbol != 1) M[symbol] += pContext->incr;
                symbol    += m;
                e         =  0;
            }
        }
        m >>= 1;
    }
    if (symbol!= 1) M[symbol] += pContext->incr;

    if (symbol & 1)
        high = low + pContext->tree[symbol];
    else {
        GET_COUNT(pContext, symbol, high);
        high += low;
    }
    if (symbol != 1) high -= pContext->incr;

#ifdef MOST_PROB_AT_END
    if (low >= pContext->most_freq_pos)   /* Ie: Was moved */
	    arithmetic_decode(low  - pContext->most_freq_count,
			  high - pContext->most_freq_count,
			  total);
    else
#endif
        arithmetic_decode(low, high, total);

    /* update the singleton count if symbol was previously a singleton */
    if (symbol != 1)
    {
      if (pContext->type == DYNAMIC && high-low == pContext->incr)
        pContext->nSingletons -= pContext->incr;

        pContext->total += pContext->incr;

        if (symbol == pContext->most_freq_symbol)
	        pContext->most_freq_count += pContext->incr;
        else if (high-low+pContext->incr > pContext->most_freq_count) { 
            pContext->most_freq_symbol = symbol;
	        pContext->most_freq_count  = high - low + pContext->incr;
	        pContext->most_freq_pos    = low;
	    } else if (symbol < pContext->most_freq_symbol)
	        pContext->most_freq_pos += pContext->incr;
    }

#ifdef MOST_PROB_AT_END
  }  /* If not MPS */
#endif

    adjust_zero_freq(pContext);

    /* halve all frequencies if necessary */
    while (pContext->total > Max_frequency)
        halve_context(pContext);

    if (symbol == 1) return NOT_KNOWN;

    return symbol-2;
}/* decode() */



/*
 *
 * Get the low and high limits of the frequency interval
 * occupied by a symbol.
 *
 */
static void 
get_interval(context *pContext, freq_value *pLow, freq_value *pHigh, int symbol)
{
    freq_value low, count;
    int p, q;
    freq_value *tree = pContext->tree;

        /* go too far */
    for(p = 1, low = 0 ; p < symbol ; ) {
        low  += tree[p], 
        p   <<= 1;
    }

        /* subtract off the extra freqs from low */
    q = symbol;
    while ((q != p) && (q <= pContext->max_length)) {
        low -= tree[q], 
        q   = FORW(q);
    }

    GET_COUNT(pContext, symbol, count);

    *pLow = low;
    *pHigh = low + count;
}
 

/*
 *
 * Halve_context is responsible for halving all the frequency counts in a 
 * context.
 * Halves context in linear time by converting tree to list of freqs
 * then back again.
 *
 * It ensures the most probable symbol size and range stay updated.
 * If halving the context gives rise to a sudden drop in the
 * ZERO_FREQ_PROB(), and if it was the MPS, it will stay recorded as the most
 * probable symbol even if it isn't.  This may cause slight compression
 * inefficiency.  (The ZERO_FREQ_PROB() as implemented does not have this
 * characteristic, and in this case the inefficiency cannot occur)
 */

static void
halve_context(context *pContext)
{
    int  shifts, p, symbol;
    freq_value incr;

    pContext->incr = (pContext->incr + MIN_INCR) >> 1;    /* halve increment */
    if (pContext->incr < MIN_INCR) 
            pContext->incr = MIN_INCR;
    pContext->nSingletons = incr = pContext->incr;
    
        /*
        ** Convert Moffat tree to array of freqs
        */
    for (shifts=0 , p = pContext->max_length ; p > 1 ; shifts++ ) p >>= 1;
    p  = 1 << shifts;      /* p is now to 2^floor(log_2 pContext->max_length) */
    while( p > 1 ) {
        symbol = p;
        while (symbol + (p >> 1) <= pContext->max_length ) {
            pContext->tree[symbol] -= pContext->tree[symbol + (p >> 1)];
            symbol                 += p;
        }
        p >>= 1;
    }

        /*
        ** Halve the counts (ignore tree[1] as it will be changed soon)
        */
    pContext->total = 0;
    for (p = 2; p <= pContext->max_length; p++) {
        pContext->tree[p] = (pContext->tree[p] + 1) >> 1;
        pContext->total  += pContext->tree[p];
        if (pContext->tree[p] == incr)
            pContext->nSingletons += incr;
    }

        /*
        ** Convert array of freqs to Moffat tree
        */
    for (p = 2; p <= pContext->max_length; ) {
        symbol = p;
        while (symbol + (p >> 1) <= pContext->max_length ) {
            pContext->tree[symbol] += pContext->tree[symbol + (p >> 1)];
            symbol                 += p;
        }
        p <<= 1;
    }

    if (pContext->type == STATIC)
        pContext->nSingletons = 0;

    pContext->tree[1] = ZERO_FREQ_PROB(pContext);
    pContext->total  += ZERO_FREQ_PROB(pContext);

    /* Recalc new most_freq_symbol info if it exists
     * (since roundoff may mean not exactly half of previous value)
     */

    if (pContext->most_freq_symbol != -1) 
      { freq_value low, high;
        get_interval(pContext, &low, &high, pContext->most_freq_symbol);
        pContext->most_freq_count = high-low;
        pContext->most_freq_pos = low;
      }

    adjust_zero_freq(pContext);
}/* halve_context() */


/*
 *
 * free memory allocated for a context and initialize empty context
 * of original size
 *
 */
void 
purge_context(context *pContext)
{
    int i;

    free(pContext->tree);
    
    /* malloc new tree of original size */
    if ((pContext->tree = (freq_value *)malloc((pContext->initial_size + 1)
                        * sizeof(freq_value))) == NULL)
    {
    fprintf(stderr, "stats: not enough memory to create context\n");
    exit(1);
    }
    pContext->length = 1;
    pContext->total = 0;
    pContext->nSymbols = 1;        /* Start with escape symbol */

    pContext->most_freq_symbol = -1;    /* Indicates no such symbol */
    pContext->most_freq_count = 0;
    pContext->most_freq_pos = 0;

    pContext->max_length = pContext->initial_size;
    for (i = 0; i <= pContext->initial_size; i++)
	    pContext->tree[i] = 0;
                      /* increment is initially 2 ^ f */
    pContext->incr = (freq_value) 1 << F_bits;
    pContext->nSingletons = 0;

    init_zero_freq(pContext);
    adjust_zero_freq(pContext);
}

/******************************************************************************
*
* functions for binary contexts
*
******************************************************************************/

/*
 *
 * create a binary_context
 * binary contexts consist of two counts and an increment which
 * is normalized
 *
 */
binary_context *create_binary_context(void)
{
    binary_context *pContext;

#ifdef VARY_NBITS
    Max_frequency = ((freq_value) 1 << F_bits);
#endif

    pContext = (binary_context *) malloc(sizeof(binary_context));
    if (pContext == NULL)
    {
    fprintf(stderr, "stats: not enough memory to create context\n");
    exit(1);
    }
                        /* start with incr=2^(f-1) */
    pContext->incr = (freq_value) 1 << (F_bits - 1);
    pContext->c0 = pContext->incr;
    pContext->c1 = pContext->incr;
    return pContext;
}



/*
 *
 * encode a binary symbol using special binary arithmetic
 * coding functions
 * returns 0 if successful
 *
 */
int
binary_encode(binary_context *pContext, int bit)
{
    binary_arithmetic_encode(pContext->c0, pContext->c1, bit);

    /* increment symbol count */
    if (bit == 0)
    pContext->c0 += pContext->incr;
    else
    pContext->c1 += pContext->incr;

    /* halve frequencies if necessary */
    if (pContext->c0 + pContext->c1 > Max_frequency)
    {
    pContext->c0 = (pContext->c0 + 1) >> 1;
    pContext->c1 = (pContext->c1 + 1) >> 1;
    pContext->incr = (pContext->incr + MIN_INCR) >> 1;
    }
    return 0;
}    



/*
 *
 * decode a binary symbol using specialised binary arithmetic
 * coding functions
 *
 */
int
binary_decode(binary_context *pContext)
{
    int bit;

    bit = binary_arithmetic_decode(pContext->c0, pContext->c1);

    /* increment symbol count */
    if (bit == 0)
    pContext->c0 += pContext->incr;
    else
    pContext->c1 += pContext->incr;

    /* halve frequencies if necessary */
    if (pContext->c0 + pContext->c1 > Max_frequency)
    {
    pContext->c0 = (pContext->c0 + 1) >> 1;
    pContext->c1 = (pContext->c1 + 1) >> 1;
    pContext->incr = (pContext->incr + MIN_INCR) >> 1;
    }    
    return bit;
}