encode.c

Arithmetic coding source code. A very good software for arithmetic coding.

    SHOW_MEM(L, ulong) /* offset */

    allocate(freq, uint, MAX_SYMBOL+1);
    allocate(syms, uint, MAX_SYMBOL+1);  // interp coding uses A[n]
    SHOW_MEM(MAX_SYMBOL+1, uint)
    SHOW_MEM(MAX_SYMBOL+1, uint)

    START_OUTPUT(out_file);
	OUTPUT_ULONG(out_file, MAGIC,sizeof(ulong)*8);
	num_header_bits += sizeof(ulong)*8;  

    if (very_verbose > 0) BLOCK_SUMMARY_HEADINGS;

    if (block_size > 0) {
        while ((b = fread(block, sizeof(uint), block_size, in_file)) > 0) {
            process_block(block, b, freq, syms, out_file);
	        num_blocks++;
        }
    } else {
        while (fread(block,sizeof(uint),1, in_file) == 1)  { /* block ahead */
            up = block;
            b = 0;
            while (*up > 0) {
                up++;
                b++;
	            if (b == current_array_size) {
                    current_array_size <<= 1;
                    if ((block = (uint *)realloc(block, current_array_size*sizeof(uint))) == NULL) {
                        fprintf(stderr,"Out of memory for block\n");
                        exit(-1);
                    }
                    up = block + b;
                }
                if (fread(up,sizeof(uint),1, in_file) != 1) {
                    fprintf(stderr,"WARNING: last block not terminated by 0\n");
                    *up = 0;
                    b--; /* exclude the 0 I add in 3 lines down */
                }
	        }
            b++; /* include the 0 */
            process_block(block, b, freq, syms, out_file);
	        num_blocks++; 
        } 
    }

    OUTPUT_ULONG(out_file, 0, LOG2_MAX_SYMBOL);  // last block indicator n = 0
    num_header_bits += LOG2_MAX_SYMBOL;

    FINISH_OUTPUT(out_file);

    if (verbose) print_summary_stats();
} /* one_pass_encoding() */

/*
** Fill freq[] with freqs
*/
uint
one_pass_freq_count(uint block[], uint b, uint freq[],uint syms[], uint ms) {
    uint n, *up;

        /* clear all elements up to max_symbol = ms */
    for(up = freq ; up <= freq + ms ; up++)         
        *up = 0;

    n = 0;

    for(up = block ; up < block + b ; up++) {
        if (freq[*up] == 0)
            syms[n++] = *up;
        freq[*up]++;
    }

    freq[EOF_SYMBOL] = 1;
    syms[n++]        = EOF_SYMBOL;

    return n;
}

/*
** INPUT:        freq[i] is the frequency of symbol i.
**               syms[0..n-1] is a list of the symbols in freq[] 
**               with a non-zero freq
** OUTPUT:       None.
** SIDE EFFECTS: Outputs some bits.
**               syms[i] contains the ordinal symbol # for symbol i.
**               freq[i] contains the codeword length for symbol i. 
**
** (1) Sort syms[0..n-1] using freq[syms[i]] as the key for syms[i]
** (2) Run in-place Huffman on freqs to overwrite freqs with codeword lens
** (3) Build cw_lens[] and then canonical coding arrays
** (4) Output n.
** (5) Sort syms in increasing order.
** (5) Interp code syms and then output unary coded of reversed max_cw-freq[i]
** (6) Overwrite syms with mapping.
*/

int pcmp(char *a, char *b) { return *((uint *)a) - *((uint *)b); }

void
build_codes(FILE *out_file, uint syms[], uint freq[], uint n) {
    uint i, *p;
    uint max_codeword_length, min_codeword_length;
    uint cw_lens[L+1];

//{uint i;
//fprintf(stderr,"*********************************************************\n");
//fprintf(stderr,"n         : %u\n",n);
//fprintf(stderr,"syms: ");
//for(i=0;i<n;i++) fprintf(stderr,"%4u ",syms[i]);
//fprintf(stderr,"\n");
//fprintf(stderr,"freq: ");
//for(i=0;i<11;i++) fprintf(stderr,"%4u ",freq[i]);
//fprintf(stderr,"\n\n");
//}
    indirect_sort(freq, syms, syms, n);

    calculate_minimum_redundancy(freq, syms, n);

//{uint i;
//fprintf(stderr,"freq: ");
//for(i=0;i<6;i++) fprintf(stderr,"%4u ",freq[i]);
//fprintf(stderr,"\n\n");
//}

        // calculcate max_codeword_length and set cw_lens[]
    for(i = 0 ; i <= L ; i++) 
        cw_lens[i] = 0;
    min_codeword_length = max_codeword_length = freq[syms[0]];
    for(p = syms ; p < syms + n ; p++) {
        if (freq[*p] > max_codeword_length)
            max_codeword_length = freq[*p];
        cw_lens[freq[*p]]++;
    }
    
    #ifdef OUTPUT_PRELUDE_CODELENGTHS
    {
        int i;
        fprintf(stderr,"%d ",max_codeword_length);
        for(i = 1 ; i <= max_codeword_length ; i++)
            fprintf(stderr,"%d ",cw_lens[i]);
        fprintf(stderr,"\n");
    }
    #endif

    build_canonical_arrays(cw_lens, max_codeword_length);

//fprintf(stderr,"*********************************************************\n");
//fprintf(stderr,"n: %10u\n",n);
//fprintf(stderr,"max_cw_len: %5u\n",max_codeword_length);
//{uint i;
//fprintf(stderr,"cw_lens : \n");
//for(i=0;i<=max_codeword_length;i++)
//fprintf(stderr,"%u\n",cw_lens[i]);
//fprintf(stderr,"\n");
//}

    OUTPUT_ULONG(out_file, n, LOG2_MAX_SYMBOL);
    OUTPUT_ULONG(out_file, max_codeword_length, LOG2_L);
    num_header_bits += LOG2_L + LOG2_MAX_SYMBOL;

    nqsort((char *)syms, n, sizeof(uint), pcmp);

    #ifdef OUTPUT_PRELUDE_SUBALPHABET_GAPS
    {
        int i;
        fprintf(stderr,"%d ",n);
        fprintf(stderr,"%d ",syms[0]);
        for(i = 1 ; i < n ; i++)
            fprintf(stderr,"%d ",syms[i]-syms[i-1]);
        fprintf(stderr,"\n");
    }
    #endif

    for(p = syms ; p < syms + n; p++) {
        OUTPUT_UNARY_CODE(out_file, max_codeword_length - freq[*p]);
        num_unary_bits += max_codeword_length - freq[*p] + 1;
    }
    interp_encode(out_file, syms, n);

    generate_mapping(cw_lens, syms, freq,  max_codeword_length, n);
} /* build_codes() */

/*
** Build lj_base[] and offset from the codelens in A[0..n-1]
** A[] need not be sorted.
**
** Return cw_lens[] a freq count of codeword lengths.
*/
void
build_canonical_arrays(uint cw_lens[], uint max_cw_length){
    ulong *q;
    uint  *p;

        // build offset
    q = offset;
    *q = 0;
    for(p = cw_lens + 1; p < cw_lens + max_cw_length ; p++, q++)
        *(q+1) = *q + *p;        

        // generate the min_code array
        // min_code[i] = (min_code[i+1] + cw_lens[i+2]) >>1
    q  = min_code + max_cw_length-1;
    *q = 0;
    for(q--, p = cw_lens + max_cw_length ; q >= min_code ; q--, p--)
        *q = (*(q+1) + *p) >> 1; 

        // generate the lj_base array
    q  = lj_base;
    unsigned long *pp = min_code;
    int left_shift = (sizeof(ulong) << 3) - 1;
    for(p = cw_lens + 1; q < lj_base + max_cw_length; p++, q++, pp++, left_shift--)
        if (*p == 0)
            *q = *(q-1);
        else
            *q = (*pp) << left_shift;
    for( p = cw_lens + 1, q = lj_base ; *p == 0 ; p++, q++)
        *q = MAX_ULONG;

//{uint i;
//for(i = 0 ; i < max_cw_length ; i++)
//    fprintf(stderr,"%3d %8lu %8lx %8lx\n",i,offset[i],min_code[i],lj_base[i]);
//fprintf(stderr,"\n");
//for(i = 1 ; i <= L ; i++)
//    fprintf(stderr,"%d ",cw_lens[i]);
//fprintf(stderr,"\n");
//}
} // build_canonical_arrays()

/*
** INPUT: syms[0..n-1] lists symbol numbers
**        freq[i] contains the codeword length of symbol i
**        cw_lens[1..max_cw_length] is the number of codewords of length i
**
** OUTPUT: None
**
** SIDE EFFECTS: syms[0..max_symbol] is overwritten with canonical code mapping.
**               cw_lens[] is destroyed.
*/
void
generate_mapping(uint cw_lens[],     uint syms[], uint freq[], 
                 uint max_cw_length, uint n) {
    int i;

    for( i = 1 ; i <= (int)max_cw_length ; i++)
        cw_lens[i] += cw_lens[i-1];     

    for(i = n - 1 ; i >= 0 ; i--)
        syms[syms[i]] = cw_lens[freq[syms[i]] - 1]++;

//{uint i;
//fprintf(stderr,"mapping\n");
//for(i = 0 ; i <= n; i++)
    //fprintf(stderr,"%8u %8u\n",syms[i],i);
//fprintf(stderr,"\n");
//}
} /* generate_mapping() */