annotate srf2fastq/io_lib-1.12.2/io_lib/deflate_interlaced.c @ 0:d901c9f41a6a default tip

Migrated tool version 1.0.1 from old tool shed archive to new tool shed repository

/*
 * This code implements "interlaced deflate", and is *based* on a
 * simplistic implementation of the Deflate algorithm.
 * See http://www.ietf.org/rfc/rfc1951.txt for details on the original
 * Deflate algorithm.
 *
 * It differs from RFC1951 in two important ways:
 *
 * 1) It only supports the huffman encoding step and does not attempt
 *    to do any LZ-style string matching to generate distance codes.
 *    (These generally do not improve data compression for our desired
 *    use.)
 *
 * 2) It optionally allows interleaving of multiple huffman trees for
 *    a single data stream. NB: when multiple codes are used this is
 *    incompatible with RFC1951.
 *
 * It has been written here, instead of using zlib, so that we can separate
 * out the encoding of the huffman tree from the compression of the data
 * stream into separate memory sections with the intent to optimise
 * compression of very small blocks of data by sharing one set of frequency
 * tables (ie huffman tree) with multiple sets of compressed data blocks.
 *
 * James Bonfield, 2007
 */

#define NDEBUG /* disable asserts for production use */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <unistd.h>
#include <assert.h>
#include <unistd.h>

#include "io_lib/deflate_interlaced.h"

#ifndef MIN
#    define MIN(a,b) ((a) < (b) ? (a) : (b))
#endif
#ifndef MAX
#    define MAX(a,b) ((a) > (b) ? (a) : (b))
#endif

/* #define TEST_MAIN */

/*
 * ---------------------------------------------------------------------------
 * Local structs & defines
 */

/* Used in tree construction only */
typedef struct node {
    int count;
    int sym; /* char or SYM_EOF */
    struct node *parent;
    struct node *next;
} node_t;

#define SYM_EOF 256

static void output_code_set(FILE *fp, huffman_codes_t *codes);
static void output_code_set2(FILE *fp, huffman_codes_t *codes);
int next_symbol(block_t *in, int *htab);

/*
 * ---------------------------------------------------------------------------
 * Our standard precomputed tables, for DNA, English text, etc.
 */

/* DNA */
/*
 * A   00
 * C   01
 * G   110
 * T   10
 * N   1110
 * EOF 11110
 * ?   11111*
 */
static huffman_code_t codes_dna[] = {
    {'A',  2}, {'C',  2}, {'T',  2}, {'G',  3}, {'N',  4}, {  0,  5},
    {SYM_EOF,  6},
    {  1, 13}, {  2, 13}, {  3, 13}, {  4, 13}, {  5, 13}, {  6, 13},
    {  7, 14}, {  8, 14}, {  9, 14}, { 10, 14}, { 11, 14}, { 12, 14},
    { 13, 14}, { 14, 14}, { 15, 14}, { 16, 14}, { 17, 14}, { 18, 14},
    { 19, 14}, { 20, 14}, { 21, 14}, { 22, 14}, { 23, 14}, { 24, 14},
    { 25, 14}, { 26, 14}, { 27, 14}, { 28, 14}, { 29, 14}, { 30, 14},
    { 31, 14}, { 32, 14}, { 33, 14}, { 34, 14}, { 35, 14}, { 36, 14},
    { 37, 14}, { 38, 14}, { 39, 14}, { 40, 14}, { 41, 14}, { 42, 14},
    { 43, 14}, { 44, 14}, { 45, 14}, { 46, 14}, { 47, 14}, {'0', 14},
    {'1', 14}, {'2', 14}, {'3', 14}, {'4', 14}, {'5', 14}, {'6', 14},
    {'7', 14}, {'8', 14}, {'9', 14}, { 58, 14}, { 59, 14}, { 60, 14},
    { 61, 14}, { 62, 14}, { 63, 14}, { 64, 14}, {'B', 14}, {'D', 14},
    {'E', 14}, {'F', 14}, {'H', 14}, {'I', 14}, {'J', 14}, {'K', 14},
    {'L', 14}, {'M', 14}, {'O', 14}, {'P', 14}, {'Q', 14}, {'R', 14},
    {'S', 14}, {'U', 14}, {'V', 14}, {'W', 14}, {'X', 14}, {'Y', 14},
    {'Z', 14}, { 91, 14}, { 92, 14}, { 93, 14}, { 94, 14}, { 95, 14},
    { 96, 14}, {'a', 14}, {'b', 14}, {'c', 14}, {'d', 14}, {'e', 14},
    {'f', 14}, {'g', 14}, {'h', 14}, {'i', 14}, {'j', 14}, {'k', 14},
    {'l', 14}, {'m', 14}, {'n', 14}, {'o', 14}, {'p', 14}, {'q', 14},
    {'r', 14}, {'s', 14}, {'t', 14}, {'u', 14}, {'v', 14}, {'w', 14},
    {'x', 14}, {'y', 14}, {'z', 14}, {123, 14}, {124, 14}, {125, 14},
    {126, 14}, {127, 14}, {128, 14}, {129, 14}, {130, 14}, {131, 14},
    {132, 14}, {133, 14}, {134, 14}, {135, 14}, {136, 14}, {137, 14},
    {138, 14}, {139, 14}, {140, 14}, {141, 14}, {142, 14}, {143, 14},
    {144, 14}, {145, 14}, {146, 14}, {147, 14}, {148, 14}, {149, 14},
    {150, 14}, {151, 14}, {152, 14}, {153, 14}, {154, 14}, {155, 14},
    {156, 14}, {157, 14}, {158, 14}, {159, 14}, {160, 14}, {161, 14},
    {162, 14}, {163, 14}, {164, 14}, {165, 14}, {166, 14}, {167, 14},
    {168, 14}, {169, 14}, {170, 14}, {171, 14}, {172, 14}, {173, 14},
    {174, 14}, {175, 14}, {176, 14}, {177, 14}, {178, 14}, {179, 14},
    {180, 14}, {181, 14}, {182, 14}, {183, 14}, {184, 14}, {185, 14},
    {186, 14}, {187, 14}, {188, 14}, {189, 14}, {190, 14}, {191, 14},
    {192, 14}, {193, 14}, {194, 14}, {195, 14}, {196, 14}, {197, 14},
    {198, 14}, {199, 14}, {200, 14}, {201, 14}, {202, 14}, {203, 14},
    {204, 14}, {205, 14}, {206, 14}, {207, 14}, {208, 14}, {209, 14},
    {210, 14}, {211, 14}, {212, 14}, {213, 14}, {214, 14}, {215, 14},
    {216, 14}, {217, 14}, {218, 14}, {219, 14}, {220, 14}, {221, 14},
    {222, 14}, {223, 14}, {224, 14}, {225, 14}, {226, 14}, {227, 14},
    {228, 14}, {229, 14}, {230, 14}, {231, 14}, {232, 14}, {233, 14},
    {234, 14}, {235, 14}, {236, 14}, {237, 14}, {238, 14}, {239, 14},
    {240, 14}, {241, 14}, {242, 14}, {243, 14}, {244, 14}, {245, 14},
    {246, 14}, {247, 14}, {248, 14}, {249, 14}, {250, 14}, {251, 14},
    {252, 14}, {253, 14}, {254, 14}, {255, 14},
};

/* DNA with a few ambiguity codes */
static huffman_code_t codes_dna_ambig[] = {
    {'A',  2}, {'C',  2}, {'T',  2}, {'G',  3}, {'N',  4}, {  0,  7},
    { 45,  7}, {'B',  8}, {'D',  8}, {'H',  8}, {'K',  8}, {'M',  8},
    {'R',  8}, {'S',  8}, {'V',  8}, {'W',  8}, {'Y',  8}, {SYM_EOF, 11},
    {226, 14}, {  1, 15}, {  2, 15}, {  3, 15}, {  4, 15}, {  5, 15},
    {  6, 15}, {  7, 15}, {  8, 15}, {  9, 15}, { 10, 15}, { 11, 15},
    { 12, 15}, { 13, 15}, { 14, 15}, { 15, 15}, { 16, 15}, { 17, 15},
    { 18, 15}, { 19, 15}, { 20, 15}, { 21, 15}, { 22, 15}, { 23, 15},
    { 24, 15}, { 25, 15}, { 26, 15}, { 27, 15}, { 28, 15}, { 29, 15},
    { 30, 15}, { 31, 15}, { 32, 15}, { 33, 15}, { 34, 15}, { 35, 15},
    { 36, 15}, { 37, 15}, { 38, 15}, { 39, 15}, { 40, 15}, { 41, 15},
    { 42, 15}, { 43, 15}, { 44, 15}, { 46, 15}, { 47, 15}, {'0', 15},
    {'1', 15}, {'2', 15}, {'3', 15}, {'4', 15}, {'5', 15}, {'6', 15},
    {'7', 15}, {'8', 15}, {'9', 15}, { 58, 15}, { 59, 15}, { 60, 15},
    { 61, 15}, { 62, 15}, { 63, 15}, { 64, 15}, {'E', 15}, {'F', 15},
    {'I', 15}, {'J', 15}, {'L', 15}, {'O', 15}, {'P', 15}, {'Q', 15},
    {'U', 15}, {'X', 15}, {'Z', 15}, { 91, 15}, { 92, 15}, { 93, 15},
    { 94, 15}, { 95, 15}, { 96, 15}, {'a', 15}, {'b', 15}, {'c', 15},
    {'d', 15}, {'e', 15}, {'f', 15}, {'g', 15}, {'h', 15}, {'i', 15},
    {'j', 15}, {'k', 15}, {'l', 15}, {'m', 15}, {'n', 15}, {'o', 15},
    {'p', 15}, {'q', 15}, {'r', 15}, {'s', 15}, {'t', 15}, {'u', 15},
    {'v', 15}, {'w', 15}, {'x', 15}, {'y', 15}, {'z', 15}, {123, 15},
    {124, 15}, {125, 15}, {126, 15}, {127, 15}, {128, 15}, {129, 15},
    {130, 15}, {131, 15}, {132, 15}, {133, 15}, {134, 15}, {135, 15},
    {136, 15}, {137, 15}, {138, 15}, {139, 15}, {140, 15}, {141, 15},
    {142, 15}, {143, 15}, {144, 15}, {145, 15}, {146, 15}, {147, 15},
    {148, 15}, {149, 15}, {150, 15}, {151, 15}, {152, 15}, {153, 15},
    {154, 15}, {155, 15}, {156, 15}, {157, 15}, {158, 15}, {159, 15},
    {160, 15}, {161, 15}, {162, 15}, {163, 15}, {164, 15}, {165, 15},
    {166, 15}, {167, 15}, {168, 15}, {169, 15}, {170, 15}, {171, 15},
    {172, 15}, {173, 15}, {174, 15}, {175, 15}, {176, 15}, {177, 15},
    {178, 15}, {179, 15}, {180, 15}, {181, 15}, {182, 15}, {183, 15},
    {184, 15}, {185, 15}, {186, 15}, {187, 15}, {188, 15}, {189, 15},
    {190, 15}, {191, 15}, {192, 15}, {193, 15}, {194, 15}, {195, 15},
    {196, 15}, {197, 15}, {198, 15}, {199, 15}, {200, 15}, {201, 15},
    {202, 15}, {203, 15}, {204, 15}, {205, 15}, {206, 15}, {207, 15},
    {208, 15}, {209, 15}, {210, 15}, {211, 15}, {212, 15}, {213, 15},
    {214, 15}, {215, 15}, {216, 15}, {217, 15}, {218, 15}, {219, 15},
    {220, 15}, {221, 15}, {222, 15}, {223, 15}, {224, 15}, {225, 15},
    {227, 15}, {228, 15}, {229, 15}, {230, 15}, {231, 15}, {232, 15},
    {233, 15}, {234, 15}, {235, 15}, {236, 15}, {237, 15}, {238, 15},
    {239, 15}, {240, 15}, {241, 15}, {242, 15}, {243, 15}, {244, 15},
    {245, 15}, {246, 15}, {247, 15}, {248, 15}, {249, 15}, {250, 15},
    {251, 15}, {252, 15}, {253, 15}, {254, 15}, {255, 15},
};

/* English text */
static huffman_code_t codes_english[] = {
    { 32,  3}, {'e',  3}, {'a',  4}, {'i',  4}, {'n',  4}, {'o',  4},
    {'s',  4}, {'t',  4}, {'d',  5}, {'h',  5}, {'l',  5}, {'r',  5},
    {'u',  5}, { 10,  6}, { 13,  6}, { 44,  6}, {'c',  6}, {'f',  6},
    {'g',  6}, {'m',  6}, {'p',  6}, {'w',  6}, {'y',  6}, { 46,  7},
    {'b',  7}, {'v',  7}, { 34,  8}, {'I',  8}, {'k',  8}, { 45,  9},
    {'A',  9}, {'N',  9}, {'T',  9}, { 39, 10}, { 59, 10}, { 63, 10},
    {'B', 10}, {'C', 10}, {'E', 10}, {'H', 10}, {'M', 10}, {'S', 10},
    {'W', 10}, {'x', 10}, { 33, 11}, {'0', 11}, {'1', 11}, {'F', 11},
    {'G', 11}, {  0, 15}, {  1, 15}, {  2, 15}, {  3, 15}, {  4, 15},
    {  5, 15}, {  6, 15}, {  7, 15}, {  8, 15}, {  9, 15}, { 11, 15},
    { 12, 15}, { 14, 15}, { 15, 15}, { 16, 15}, { 17, 15}, { 18, 15},
    { 19, 15}, { 20, 15}, { 21, 15}, { 22, 15}, { 23, 15}, { 24, 15},
    { 25, 15}, { 26, 15}, { 27, 15}, { 28, 15}, { 29, 15}, { 30, 15},
    { 31, 15}, { 35, 15}, { 36, 15}, { 37, 15}, { 38, 15}, { 40, 15},
    { 41, 15}, { 42, 15}, { 43, 15}, { 47, 15}, {'2', 15}, {'3', 15},
    {'4', 15}, {'5', 15}, {'6', 15}, {'7', 15}, {'8', 15}, {'9', 15},
    { 58, 15}, { 60, 15}, { 61, 15}, { 62, 15}, { 64, 15}, {'D', 15},
    {'J', 15}, {'K', 15}, {'L', 15}, {'O', 15}, {'P', 15}, {'Q', 15},
    {'R', 15}, {'U', 15}, {'V', 15}, {'X', 15}, {'Y', 15}, {'Z', 15},
    { 91, 15}, { 92, 15}, { 93, 15}, { 94, 15}, { 95, 15}, { 96, 15},
    {'j', 15}, {'q', 15}, {'z', 15}, {123, 15}, {124, 15}, {125, 15},
    {126, 15}, {127, 15}, {128, 15}, {129, 15}, {130, 15}, {131, 15},
    {132, 15}, {133, 15}, {134, 15}, {135, 15}, {136, 15}, {137, 15},
    {138, 15}, {139, 15}, {140, 15}, {141, 15}, {142, 15}, {143, 15},
    {144, 15}, {145, 15}, {146, 15}, {147, 15}, {148, 15}, {149, 15},
    {150, 15}, {151, 15}, {152, 15}, {153, 15}, {154, 15}, {155, 15},
    {156, 15}, {157, 15}, {158, 15}, {159, 15}, {160, 15}, {161, 15},
    {162, 15}, {163, 15}, {164, 15}, {165, 15}, {166, 15}, {167, 15},
    {168, 15}, {169, 15}, {170, 15}, {171, 15}, {172, 15}, {173, 15},
    {174, 15}, {175, 15}, {176, 15}, {177, 15}, {178, 15}, {179, 15},
    {180, 15}, {181, 15}, {182, 15}, {183, 15}, {184, 15}, {185, 15},
    {186, 15}, {187, 15}, {188, 15}, {189, 15}, {190, 15}, {191, 15},
    {192, 15}, {193, 15}, {194, 15}, {195, 15}, {196, 15}, {197, 15},
    {198, 15}, {199, 15}, {200, 15}, {201, 15}, {202, 15}, {203, 15},
    {204, 15}, {205, 15}, {206, 15}, {207, 15}, {208, 15}, {209, 15},
    {210, 15}, {211, 15}, {212, 15}, {213, 15}, {214, 15}, {215, 15},
    {216, 15}, {217, 15}, {218, 15}, {219, 15}, {220, 15}, {221, 15},
    {222, 15}, {223, 15}, {224, 15}, {225, 15}, {226, 15}, {227, 15},
    {228, 15}, {229, 15}, {230, 15}, {231, 15}, {232, 15}, {233, 15},
    {234, 15}, {235, 15}, {236, 15}, {237, 15}, {238, 15}, {239, 15},
    {240, 15}, {241, 15}, {242, 15}, {243, 15}, {244, 15}, {245, 15},
    {246, 15}, {247, 15}, {248, 15}, {249, 15}, {250, 15}, {251, 15},
    {252, 15}, {253, 15}, {254, 15}, {255, 15}, {SYM_EOF, 15},
};

static huffman_codeset_t *static_codeset[NCODES_STATIC];

/*
 * ---------------------------------------------------------------------------
 * Block_t structure support
 */

/*
 * Allocates and returns a new block_t struct of a specified default size.
 * A default 'data' pointer may be passed in, in which it must have
 * been created using malloc(size). Otherwise if data is NULL then
 * size indicates the amount of memory to allocate. Size maybe zero to
 * defer allocation.
 *
 * Returns newly created block_t* on success
 *         NULL on failure
 */
block_t *block_create(unsigned char *data, size_t size) {
    block_t *b = (block_t *)malloc(sizeof(*b));
    if (!b)
	return NULL;

    b->data = data;
    b->alloc = size;
    b->byte = 0;
    b->bit = 0;

    if (size && !data && NULL == (b->data = calloc(size, 1))) {
	free(b);
	return NULL;
    }

    return b;
}

/*
 * Deallocates memory created by block_create().
 * keep_data is a boolean which if true requests that the data held within
 * the block should not be deallocated as it is in use elsewhere.
 */
void block_destroy(block_t *blk, int keep_data) {
    if (!blk)
	return;

    if (!keep_data && blk->data)
	free(blk->data);
    
    free(blk);
}

/*
 * Ensures a block_t holds at least 'size' bytes.
 * Newly allocated data is initialised to zero.
 *
 * Returns 0 on success
 *        -1 on failure, leaving block pointing to the existing data
 */
int block_resize(block_t *blk, size_t size) {
    unsigned char *newp = NULL;

    if (!blk)
	return -1;

    /* Grow size to next power of 2, if we're growing */
    if (size > blk->alloc) {
	size--;
	size |= size >> 1;
	size |= size >> 2;
	size |= size >> 4;
	size |= size >> 8;
	size |= size >> 16;
	size++;
    }

    if (NULL == (newp = realloc(blk->data, size)))
	return -1;
    else
	blk->data = newp;

    if (size > blk->alloc)
	memset(&blk->data[blk->alloc], 0, size - blk->alloc);
    blk->alloc = size;

    return 0;
}


/*
 * ---------------------------------------------------------------------------
 * Tree building and code generation functions
 */

/*
 * Reverses the order of bits in the bottom nbits of val.
 * Returns the bit-reverse value.
 */
unsigned int bit_reverse(unsigned int val, int nbits) {
    unsigned int new = 0, i;

    for (i = 0; i < nbits; i++) {
	new = (new << 1) | (val & 1);
	val >>= 1;
    }

    return new;
}


/*
 * Generates canonical huffman codes given a set of symbol bit lengths.
 * The results are stored within the supplied huffman_codes_t struct.
 *
 * Returns 0 on success
 *        -1 on failure
 */
static int canonical_codes(huffman_codes_t *c) {
    int i, j;
    unsigned int code, last_len;
    int clens[33];
    int offs[33];
    huffman_code_t ctmp[258];
    signed int symtab[258];

    /* Sort by bit-length, subfield symbol - much faster than qsort() */
    for (i = 0; i < 258; i++)
	symtab[i] = -1;
    for (i = 0; i < c->ncodes; i++)
	symtab[c->codes[i].symbol] = i;
    for (i = 0; i <= 32; i++)
	offs[i] = clens[i] = 0;
    for (i = 0; i < c->ncodes; i++)
	clens[c->codes[i].nbits]++;
    for (i = 1; i <= 32; i++)
	offs[i] = offs[i-1] + clens[i-1];
    for (i = 0; i < 258; i++) {
	if (symtab[i] != -1)
	    ctmp[offs[c->codes[symtab[i]].nbits]++] = c->codes[symtab[i]];
    }
    memcpy(c->codes, ctmp, c->ncodes * sizeof(huffman_code_t));

    /*
     * Force all codes to be <= max_code_len. This is needed due to the
     * 15-bit length limitation of Deflate literal codes and the 7-bit 
     * limit of the code bit-length table.
     */
    /* Find first point of failure */
    for (i = 0; i < c->ncodes; i++) {
	if (c->codes[i].nbits > c->max_code_len)
	    break;
    }

    /*
     * From here on we shrink the length of the current code by increasing
     * the length of an earlier symbol, at last_code.
     */
    if (i != c->ncodes) {
	int delta = 0;

	/*
	fprintf(stderr, "=== REORDERING %d ===\n", c->code_set);
	output_code_set(stderr, c);
	output_code_set2(stderr, c);
	*/

	for (; i < c->ncodes; i++) {
	    int k, cur_len;

	    c->codes[i].nbits -= delta;
	    if (c->codes[i].nbits <= c->max_code_len)
		continue;

	    for (j = i; j >= 0 && c->codes[j].nbits >= c->max_code_len; j--)
		;
	    if (j < 0) {
		fprintf(stderr,
			"Too many symbols to fit in bit-length requirements\n");
		fprintf(stderr, "=== FAILING ===\n");
		output_code_set(stderr, c);
		output_code_set2(stderr, c);
		abort();
	    }

	    /*
	    fprintf(stderr, "Changing code %d/%d to len %d\n",
		    c->codes[i].symbol, c->codes[j].symbol,
		    c->codes[j].nbits+1);
	    */
	    cur_len = c->codes[i].nbits;
	    c->codes[i].nbits = ++c->codes[j].nbits;

	    /*
	     * Shrink the next code by one, or if none at that bit-length
	     * the next 2, and so on
	     */
	    delta = 1;
	    for (k = i+1; delta && k < c->ncodes; k++) {
		while (c->codes[k].nbits > cur_len) {
		    delta *= 2;
		    cur_len++;
		}
		c->codes[k].nbits--;
		delta--;
	    }
	    assert(delta == 0);
	}

	/*
	fprintf(stderr, "=== REORDERED TO %d ===\n", c->code_set);
	output_code_set(stderr, c);
	output_code_set2(stderr, c);
	*/

	/* Ordering is shot - regenerate via brute force way */
	return canonical_codes(c);
    }


    /* Generate codes */
    code = last_len = 0; /* stop warning */
    for (i = 0; i < c->ncodes; i++) {
	int nbits = c->codes[i].nbits;

	if (i == 0) {
	    code = 0;
	    last_len = nbits;
	} else {
	    code++;
	}
	if (nbits > last_len) {
	    code <<= (nbits - last_len);
	    last_len = nbits;
	}
	c->codes[i].code = bit_reverse(code, nbits);
    }

    /* Reindex so the symbol is the primary index into codes */
    for (i = 0; i <= 257; i++) {
	c->lookup[i].nbits = 0;
    }
    for (i = 0; i < c->ncodes; i++) {
        c->lookup[c->codes[i].symbol] = c->codes[i];
    }

    return 0;
}

static int node_compar2(const void *vp1, const void *vp2) {
    const node_t *n1 = *(const node_t **)vp1;
    const node_t *n2 = *(const node_t **)vp2;

    /*
     * The sort order is vital here. This needs to return the same collating
     * order on all systems so that differing qsort() functions will not
     * swap around symbols with the same bit lengths, hence we sort by both
     * fields to force a unique stable ordering.
     */
    if (n1->count != n2->count)
	return n1->count - n2->count;
    else
	return n2->sym - n1->sym;
}

/*
 * Computes the huffman bit-lengths for a data set. We don't care
 * about the actual tree, just how deep the symbols end up.
 *
 * Huffman trees are constructed by constructing a set of nodes
 * initially containing the symbol and it's frequency. We then merge
 * the two least used nodes to produce a new node with a combined
 * frequency. Repeat until one root node is left.
 *
 * data/len is the input data to analyse.
 *
 * 'eof' is a boolean to indicate whether the EOF symbol should be included
 * in the symbols produced.
 *
 * all_codes is a boolean to indicate whether we should include symbols not
 * found in the input data set. (This was used to create the static lookup
 * tables.)
 *
 * Returns huffman_codes_t* on success
 *         NULL on failure
 */
huffman_codes_t *calc_bit_lengths(unsigned char *data, int len,
				  int eof, int max_code_len, int all_codes,
				  int start, int skip) {
    int i, ncodes;
    node_t nodes[258+257], *head, *new = &nodes[258];
    node_t *n2[258+257];
    int map[258];
    huffman_codes_t *c;
    int hist[256];

    if (NULL == (c = (huffman_codes_t *)malloc(sizeof(*c))))
	return NULL;
    c->codes_static = 0;
    c->max_code_len = max_code_len;

    /* Count frequencies of symbols */
    memset(hist, 0, 256*sizeof(*hist));
    /* Calc freqs */
    for (i = start; i < len; i+=skip) {
	hist[data[i]]++;
    }


    /*
     * Initialise nodes. We build a map of ASCII character code to node
     * number. (By default it's a simple 1:1 mapping unless legal_chars is
     * defined.)
     */
    ncodes = 0;
    if (eof) {
	nodes[ncodes].sym = SYM_EOF;
	nodes[ncodes].count = eof;
	nodes[ncodes].parent = NULL;
	n2[ncodes] = &nodes[ncodes];
	map[SYM_EOF] = ncodes++;
    }

    /* All 256 chars existing at a minimal level */
    if (all_codes) {
	for (i = 0; i < 256; i++) {
	    nodes[ncodes].sym = i;
	    nodes[ncodes].count = hist[i];
	    nodes[ncodes].parent = NULL;
	    n2[ncodes] = &nodes[ncodes];
	    map[i] = ncodes++;
	}
    } else {
	/* Only include non-zero symbols */
	for (i = 0; i < 256; i++) {
	    if (hist[i] == 0)
		continue;
	    nodes[ncodes].sym = i;
	    nodes[ncodes].count = hist[i];
	    nodes[ncodes].parent = NULL;
	    n2[ncodes] = &nodes[ncodes];
	    map[i] = ncodes++;
	}
    }

    /* Sort by counts, smallest first and form a sorted linked list */
    qsort(n2, ncodes, sizeof(*n2), node_compar2);

    /* Skip symbols that do not occur, unless all_codes is true */
    for (i = 0; i < ncodes; i++) {
	n2[i]->next = i+1 < ncodes ? n2[i+1] : NULL;
    }

    /* Repeatedly merge two smallest values */
    head = n2[0];
    while (head && head->next) {
	node_t *after = head->next, *n;
	int sum = head->count + head->next->count;
	
	for (n = head->next->next; n; after = n, n = n->next) {
	    if (sum < n->count)
		break;
	}

	/* Produce a new summation node and link it in place */
	after->next = new;
	new->next = n;
	new->sym = '?';
	new->count = sum;
	new->parent = NULL;
	head->parent = new;
	head->next->parent = new;
	head = head->next->next;

	new++;
    }

    /* Walk up tree computing the bit-lengths for our symbols */
    c->ncodes = ncodes;
    c->codes = (huffman_code_t *)malloc(c->ncodes * sizeof(*c->codes));
    if (NULL == c->codes) {
	free(c);
	return NULL;
    }

    for (i = 0; i < ncodes; i++) {
	int len = 0;
	node_t *n;
	for (n = n2[i]->parent; n; n = n->parent) {
	    len++;
	}

	c->codes[i].symbol = n2[i]->sym;
	c->codes[i].freq   = n2[i]->count;
	c->codes[i].nbits  = len ? len : 1; /* special case, nul input */
    }

    if (0 != canonical_codes(c)) {
	free(c);
	return NULL;
    }

    return c;
}

/*
 * A special case of the generate_code_set() function below, but for
 * creating predefined code sets from bit-length arrays. Useful for
 * code that wants to use a predetermined huffman tree.
 *
 *
 * Returns huffman_codes_t* on success; free using huffman_codes_destroy().
 *         NULL on failure.
 */
huffman_codeset_t *codes2codeset(huffman_code_t *codes, int ncodes,
				 int code_num) {
    huffman_codeset_t *cs;
    huffman_codes_t *c;

    if (NULL == (cs = (huffman_codeset_t *)malloc(sizeof(*cs))))
	return NULL;

    if (NULL == (c = (huffman_codes_t *)malloc(sizeof(*c))))
	return NULL;

    cs->codes = (huffman_codes_t **)malloc(sizeof(*cs->codes));
    cs->codes[0] = c;
    cs->ncodes = 1;
    cs->code_set = code_num;
    cs->blk = NULL;
    cs->bit_num = 0;
    cs->decode_t = NULL;
    cs->decode_J4 = NULL;
    
    c->codes_static = 1;
    c->max_code_len = MAX_CODE_LEN;

    c->codes = codes;
    c->ncodes = ncodes;

    cs->bit_num = 0; /* FIXME: need to know this */

    canonical_codes(c);

    return cs;
}

/*
 * Initialises and returns a huffman_codes_t struct from a specified code_set.
 * If code_set is not one of the standard predefined values then the
 * input data is analysed using calc_bit_lengths() above to produce the
 * optimal set of huffman codes, otherwise we return predefined values.
 *
 * 'eof' is a boolean to indicate whether the EOF symbol should be included
 * in the symbols produced.
 *
 * all_codes is a boolean to indicate whether we should include symbols not
 * found in the input data set. (This was used to create the static lookup
 * tables.)
 *
 * Returns huffman_codes_t* on success; free using huffman_codes_destroy().
 *         NULL on failure.
 */
huffman_codeset_t *generate_code_set(int code_set, int ncodes,
				     unsigned char *data, int len,
				     int eof, int max_code_len,
				     int all_codes) {
    huffman_codeset_t *cs;

    /*
     * Either INLINE or a CODE_USER+ set of codes.
     * => analyse the data and compute a new set of bit-lengths & codes.
     */
    if (code_set >= 128 || code_set == CODE_INLINE) { 
	int i;

	if (NULL == (cs = (huffman_codeset_t *)malloc(sizeof(*cs))))
	    return NULL;

	cs->code_set = code_set;
	cs->ncodes = ncodes;
	cs->codes = (huffman_codes_t **)malloc(cs->ncodes*sizeof(*cs->codes));
	cs->blk = NULL;
	cs->bit_num = 0;
	cs->decode_t = NULL;
	cs->decode_J4 = NULL;

	for (i = 0; i < ncodes; i++) {
	    /*
	     * If requested, include EOF all code sets, but at a
	     * frequency of only '1' occurrance where we predict it
	     * not to be needed.
	     */
	    if (eof && (len+i)%ncodes)
		eof = 1;
	    cs->codes[i] = calc_bit_lengths(data, len, eof, max_code_len,
					    all_codes, i, ncodes);
	    cs->codes[i]->codes_static = 0;
	    if (NULL == cs->codes[i]) {
		/* FIXME: tidy up */
		return NULL;
	    }

	    canonical_codes(cs->codes[i]);
	}

    /*
     * Otherwise we use the determined codes at the top of this file, such
     * as codes_dna and codes_english.
     */
    } else {
	if (code_set < 1 || code_set >= NCODES_STATIC) {
	    fprintf(stderr, "Unknown huffman code set '%d'\n", code_set);
	    return NULL;
	}

	/* If our global codeset hasn't been initialised yet, do so */
	if (!static_codeset[code_set]) {
	    huffman_codes_t *c = (huffman_codes_t *)malloc(sizeof(*c));

	    if (NULL == (cs = (huffman_codeset_t *)malloc(sizeof(*cs))))
		return NULL;

	    cs->codes = (huffman_codes_t **)malloc(sizeof(*cs->codes));
	    cs->codes[0] = c;
	    cs->ncodes = 1;
	    cs->code_set = code_set;
	    cs->blk = NULL;
	    cs->bit_num = 0;
	    cs->decode_t = NULL;
	    cs->decode_J4 = NULL;

	    c->codes_static = 1;
	    c->max_code_len = MAX_CODE_LEN;

	    switch(code_set) {
	    case CODE_DNA:
		c->codes = codes_dna;
		c->ncodes = sizeof(codes_dna)/sizeof(*c->codes);
		cs->bit_num = 5;
		break;

	    case CODE_DNA_AMBIG:
		c->codes = codes_dna_ambig;
		c->ncodes = sizeof(codes_dna_ambig)/sizeof(*c->codes);
		cs->bit_num = 1;
		break;

	    case CODE_ENGLISH:
		c->codes = codes_english;
		c->ncodes = sizeof(codes_english)/sizeof(*c->codes);
		cs->bit_num = 1;
		break;

	    default:
		fprintf(stderr, "Unknown huffman code set '%d'\n", code_set);
		return NULL;
	    }

	    canonical_codes(c);

	    static_codeset[code_set] = cs;
	}

	cs = static_codeset[code_set];
    }

    return cs;
}

void huffman_codes_destroy(huffman_codes_t *c) {
    if (!c)
	return;

    if (!c->codes_static && c->codes)
	free(c->codes);

    free(c);
}

void huffman_codeset_destroy(huffman_codeset_t *cs) {
    int i;

    if (!cs)
	return;

    /* If this codeset is one of the predefined global ones we do nothing */
    if (cs->ncodes == 1 && cs->codes[0]->codes_static)
	return;

    for (i = 0; i < cs->ncodes; i++)
	huffman_codes_destroy(cs->codes[i]);
    if (cs->codes)
	free(cs->codes);
    if (cs->blk)
	block_destroy(cs->blk, 0);
    if (cs->decode_t)
	free(cs->decode_t);
    if (cs->decode_J4)
	free(cs->decode_J4);

    free(cs);
}

/*
 * ---------------------------------------------------------------------------
 * Encoding and decoding related functions
 */

/*
 * Can store up to 24-bits worth of data encoded in an integer value
 * Possibly we'd want to have a less optimal store_bits function when dealing
 * with nbits > 24, but for now we assume the codes generated are never
 * that big. (Given this is only possible with 121392 or more
 * characters with exactly the correct frequency distribution we check
 * for it elsewhere.)
 */
static void store_bits(block_t *block, unsigned int val, int nbits) {
    /* fprintf(stderr, " store_bits: %02x %d\n", val, nbits); */

#if 1
    {
	unsigned int curr = block->data[block->byte];
	curr |= (val & ((1 << nbits)-1)) << block->bit;
	block->bit += nbits;
	while (block->bit >= 8) {
	    block->data[block->byte++] = curr & 0xff;
	    curr >>= 8;
	    block->bit -= 8;
	}
	block->data[block->byte] = curr & 0xff;
    }
    return;

#else

    {
      /* Slow, but simple */
      unsigned int mask = 1;
      int bit = 1 << block->bit;
      do {
	if (val & mask)
	    block->data[block->byte] |= bit;
	/*
	 * Data should be zeroed anyway, so this is not needed.
	 *
	 * else
	 *    block->data[block->byte] &= ~bit;
	*/

	if (++block->bit == 8) {
	    block->bit = 0;
	    block->byte++;
	    block->data[block->byte] = 0;
	    bit = 1;
	} else {
	    bit <<= 1;
	}
	mask <<= 1;
      } while(--nbits);
    }

#endif
}

/*
 * Reads up to 24-bits worth of data and returns. Updates the block
 * byte and bit values to indicate the current 'read' position.
 *
 * Returns unsigned value on success (>=0)
 *         -1 on failure
 */
static signed int get_bits(block_t *block, int nbits) {
    unsigned int val, bnum = 0;

    if (block->byte*8 + block->bit + nbits > block->alloc * 8)
	return -1;

    /* Fetch the partial byte of data */
    val = (block->data[block->byte]) >> block->bit;
    bnum = 8 - block->bit;

    /* And additional entire bytes worth as required */
    while (bnum <= nbits) {
	val |= block->data[++block->byte] << bnum;
	bnum += 8;
    }

    block->bit = (block->bit + nbits) % 8;
    return val & ((1 << nbits) - 1);
}

/* stores nbytes bytes, padding to align on the next byte boundary */
void store_bytes(block_t *block, unsigned char *val, int nbytes) {
    /* Align */
    if (block->bit != 0) {
	block->byte++;
	block->bit = 0;
    }

    /* Resize */
    block_resize(block, block->byte + nbytes + 1);

    /* Store */
    memcpy(&block->data[block->byte], val, nbytes);
    block->byte += nbytes;
}


/*
 * Encodes the huffman symbol bit-lengths as a serialised block of data
 * suitable for storing in a ZTR "ZLBH" chunk. This uses the Deflate
 * storage format defined in RFC1951.
 *
 * Returns: 0 on success
 *         -1 on failure
 */
int store_codes_single(block_t *out, huffman_codes_t *c) {
    int i;
    unsigned char bl_code[257]; /* bit-length codes and for codes 16-18 */
    unsigned char bl_opt[257];  /*     the operand to the blcode */
    unsigned char sorted_codes[258];
    int bl_freq[19]; /* frequency of bit-length codes produced */
    int bl_count;
    huffman_codes_t *bl_cds = NULL;
    int hclen_order[] = {
	16, 17, 18, 0, 8, 7, 9, 6, 10, 5, 11, 4, 12, 3, 13, 2, 14, 1, 15
    };
    int hlit, hdist, hclen, hcmap[19];

    /* output_code_set (stderr, c); */

    if (out->alloc < out->byte + 1000) {
	out->alloc = out->byte + 1000;
	if (NULL == (out->data = realloc(out->data, out->alloc)))
	    return -1;
    }

    /*
     *-----------------------------------------------------------------
     * Reformat the dynamic code bit-lengths into an alphabet of 19 
     * "code length" symbols as defined in RFC1951.
     */
    memset(sorted_codes, 0, 258);
    for (i = 0; i < c->ncodes; i++) {
	sorted_codes[c->codes[i].symbol] = c->codes[i].nbits;
    }
    for (i = 0; i < 19; i++)
	bl_freq[i] = 0;

    bl_count = 0;
    for (i = 0; i < 257; ) {
	int j = i+1, n;
	int v = sorted_codes[i];
	while (j < 257 && sorted_codes[j] == v)
	    j++;

	n = j-i; /* n = run-length */
	/* fprintf(stderr, "value=%d, run_len=%d\n", v, n); */
	if (v == 0) {
	    /* bit-len zero => no code and uses code 17/18 for run len */
	    while (n > 0) {
		while (n >= 11) {
		    bl_freq[18]++;
		    bl_code[bl_count] = 18;
		    bl_opt[bl_count] = MIN(n, 138)-11;
		    n -= bl_opt[bl_count]+11;
		    bl_count++;
		}

		while (n >= 3) {
		    bl_freq[17]++;
		    bl_code[bl_count] = 17;
		    bl_opt[bl_count] = MIN(n, 10)-3;
		    n -= bl_opt[bl_count]+3;
		    bl_count++;
		}

		switch (n) {
		case 2: bl_code[bl_count++] = 0; bl_freq[0]++; n--;
		case 1: bl_code[bl_count++] = 0; bl_freq[0]++; n--;
		}
	    }
	} else if (v <= 15) {
	    /* non-zero code, uses code 16 for run-len */
	    if (n >= 4) {
		bl_freq[v]++;
		bl_code[bl_count++] = v;
		n--;
		while (n >= 3) {
		    bl_freq[16]++;
		    bl_code[bl_count] = 16;
		    bl_opt[bl_count] = MIN(n, 6)-3;
		    n -= bl_opt[bl_count]+3;
		    bl_count++;
		}
	    }

	    switch(n) {
	    case 3: bl_code[bl_count++] = v; bl_freq[v]++; n--;
	    case 2: bl_code[bl_count++] = v; bl_freq[v]++; n--;
	    case 1: bl_code[bl_count++] = v; bl_freq[v]++; n--;
	    }
	} else {
	    fprintf(stderr, "Unsupported code length: %d\n", v);
	}

	i = j;
    }
    hlit = 257;

    /* Add a single distance code of zero bits. This means that there
     * are no distance codes used at all.
     */
    bl_code[bl_count++] = 0;
    bl_freq[0]++;
    hdist = 1;

    /* Produce new huffman codes for our code-length symbols. */
    bl_cds = calc_bit_lengths(bl_code, bl_count, 0, 7, 0, 0, 1);
    /* output_code_set (stderr, bl_cds); */

    /*
     *-----------------------------------------------------------------
     * Output the "code length" bit-lengths, 3 bits at a time.
     *
     * Compute how many HCLEN code length values we need, using the
     * predefined order in the RFC.
     */
    for (hclen = 19; hclen > 0; hclen--) {
	if (bl_freq[hclen_order[hclen-1]])
	    break;
    }

    store_bits(out, hlit-257,  5);
    store_bits(out, hdist-1, 5);
    store_bits(out, hclen-4, 4);

    for (i = 0; i < 19; i++)
	hcmap[i] = -1;
    for (i = 0; i < bl_cds->ncodes; i++)
	hcmap[bl_cds->codes[i].symbol] = i;
    for (i = 0; i < hclen; i++) {
	if (hcmap[hclen_order[i]] >= 0) {
	    store_bits(out, bl_cds->codes[hcmap[hclen_order[i]]].nbits, 3);
	} else {
	    store_bits(out, 0, 3);
	}
    }


    /*
     *----------------------------------------------------------------
     * Finally output the original bit-lengths using the code-len codes.
     */
    for (i = 0; i < bl_count; i++) {
	huffman_code_t *c = &bl_cds->codes[hcmap[bl_code[i]]];
	store_bits(out, c->code, c->nbits); 
	/*
	fprintf(stderr, "bl_code %d (opt %d), code %d/%d\n",
		bl_code[i], bl_opt[i], c->code, c->nbits);
	*/
	switch(bl_code[i]) {
	case 18:
	    store_bits(out, bl_opt[i], 7);
	    break;
	case 17:
	    store_bits(out, bl_opt[i], 3);
	    break;
	case 16:
	    store_bits(out, bl_opt[i], 2);
	    break;
	}
    }

    huffman_codes_destroy(bl_cds);

    return 0;
}

/*
 * A wrapper around store_codes_single to output either a single or multiple
 * huffman codes to a block.
 *
 * This also creates a new block and fills out the block header appropriately.
 *
 * Returns 0 on success
 *        -1 on failure
 */
int store_codes(block_t *out, huffman_codeset_t *cs, int last_block) {
    int i;

    if (out->alloc < out->byte + 1000) {
	out->alloc = out->byte + 1000;
	if (NULL == (out->data = realloc(out->data, out->alloc)))
	    return -1;
    }

    /* Header details */
    store_bits(out, last_block != 0,  1); /* last block */
    if (cs->ncodes == 1) {
	store_bits(out, 2,  2); /* dynamic huffman */
    } else {
	int nbits = 0;
	store_bits(out, 3,  2); /* multiple tree dynamic huffman */
	while (1<<nbits <= cs->ncodes-1)
	    nbits++;
	store_bits(out, nbits-1, 4);
	store_bits(out, cs->ncodes-1, nbits);
    }

    for (i = 0; i < cs->ncodes; i++) {
	if (-1 == store_codes_single(out, cs->codes[i]))
	    return -1;
    }

    return 0;
}

/*
 * This is the opposite of the store_codes() function. It loads generates
 * huffman_codes_t structs from the a serialised data stream as presented
 * in the above format.
 *
 * The input data is the data-string. On return the number of bytes
 * consumed will be returned in *len_used (if non NULL).
 * This is to allow stripping off of the huffman codes from a longer
 * array of data (ie probably followed by the STHUFF encoded chunk
 * itself).
 *
 * Returns: malloced huffman_codes_t structure on success.
 *          NULL on failure.
 */
huffman_codes_t *restore_codes_single(block_t *block) {
    int hlit, hdist, hclen;
    int hclen_order[19] = {
	16, 17, 18, 0, 8, 7, 9, 6, 10, 5, 11, 4, 12, 3, 13, 2, 14, 1, 15
    };
    int hc_bitlen[19], i;
    huffman_codes_t *bl_cds, *cds;
    int sym, sym_val, last_len;
    int htab[256];

    hlit   = get_bits(block, 5)+257;
    hdist  = get_bits(block, 5)+1;
    hclen  = get_bits(block, 4)+4;

    /*
    fprintf(stderr, "bfinal = %d, btype=%d\n", *bfinal, btype);
    fprintf(stderr, "hlit=0x%x, hdist=0x%x, hclen=0x%x\n",
	    hlit, hdist, hclen);
    */

    /* Read HCLEN code-lengths and construct huffman codes from them */
    for (i = 0; i < hclen; i++)
	hc_bitlen[hclen_order[i]] = get_bits(block, 3);
    for (; i < 19; i++)
	hc_bitlen[hclen_order[i]] = 0;

    bl_cds = (huffman_codes_t *)malloc(sizeof(*bl_cds));
    bl_cds->codes_static = 0;
    bl_cds->ncodes = 0;
    bl_cds->codes = (huffman_code_t *)malloc(19 * sizeof(*bl_cds->codes));
    bl_cds->max_code_len = 7;
    for (i = 0; i < 19; i++) {
	if (hc_bitlen[i]) {
	    bl_cds->codes[bl_cds->ncodes].symbol = i;
	    bl_cds->codes[bl_cds->ncodes].nbits = hc_bitlen[i];
	    bl_cds->ncodes++;
	}
    }

    canonical_codes(bl_cds);
    /* output_code_set (stderr, bl_cds); */

    /* Build a lookup table of possible codes to symbols */
    for (i = 0; i < 256; i++)
	htab[i] = -1;
    for (i = 0; i < bl_cds->ncodes; i++) {
	htab[bit_reverse(bl_cds->codes[i].code, bl_cds->codes[i].nbits)
	     | (1<<bl_cds->codes[i].nbits)]
	    = bl_cds->codes[i].symbol;
    }

    /* Now decode the next HLIT literal codes using bl_cds */
    cds = (huffman_codes_t *)malloc(sizeof(*cds));
    cds->codes_static = 0;
    cds->ncodes = 0;
    cds->codes = (huffman_code_t *)malloc(257 * sizeof(*cds->codes));
    cds->max_code_len = 15;
    sym_val = last_len = 0;
    while ((sym = next_symbol(block, htab)) != -1) {
	int count;
	/* fprintf(stderr, "LIT Sym=%d\n", sym); */

	switch(sym) {
	case 16:
	    count = get_bits(block, 2) + 3;
	    /* fprintf(stderr, "   +%d\n", count); */
	    for (i = 0; i < count; i++) {
		cds->codes[cds->ncodes].symbol = sym_val++;
		cds->codes[cds->ncodes++].nbits = last_len;
	    }
	    break;
	    
	case 17:
	    count = get_bits(block, 3) + 3;
	    /* fprintf(stderr, "   +%d\n", count); */
	    sym_val += count;
	    last_len = 0;
	    break;
	    
	case 18:
	    count = get_bits(block, 7) + 11;
	    /* fprintf(stderr, "   +%d\n", count); */
	    sym_val += count;
	    last_len = 0;
	    break;

	case 0:
	    sym_val++;
	    last_len = 0;
	    break;

	default:
	    cds->codes[cds->ncodes].symbol = sym_val++;
	    last_len = cds->codes[cds->ncodes++].nbits = sym;
	}

	if (sym_val >= hlit)
	    break;
    }
    assert(sym != -1);
    assert(cds->ncodes <= 257);

    /* Skip HDIST codes. Hopefully only 1 of zero length */
    sym_val = 0;
    while ((sym = next_symbol(block, htab)) != -1) {
	/* fprintf(stderr, "DIST Sym=%d\n", sym); */

	switch(sym) {
	case 16:
	    sym_val += get_bits(block, 2) + 3;
	    break;
	    
	case 17:
	    sym_val += get_bits(block, 3) + 3;
	    break;
	    
	case 18:
	    sym_val += get_bits(block, 7) + 11;
	    break;

	default:
	    sym_val++;
	}

	if (sym_val >= hdist)
	    break;
    }
    assert(sym != -1);

    huffman_codes_destroy(bl_cds);
    canonical_codes(cds);
    /* output_code_set(stderr, cds); */

    return cds;
}

huffman_codeset_t *restore_codes(block_t *block, int *bfinal) {
    int btype;
    huffman_codeset_t *cs;
    
    /* Header details */
    if (bfinal)
	*bfinal = get_bits(block, 1);
    else
	get_bits(block, 1);
    btype  = get_bits(block, 2);

    cs = (huffman_codeset_t *)malloc(sizeof(*cs));
    cs->code_set = 0;
    cs->blk = NULL;
    cs->bit_num = 0;
    cs->decode_t = NULL;
    cs->decode_J4 = NULL;

    if (btype == 2) {
	/* Standard Deflate algorithm */
	cs->ncodes = 1;
	cs->codes = (huffman_codes_t **)malloc(cs->ncodes*sizeof(*cs->codes));
	cs->codes[0] = restore_codes_single(block);
    } else if (btype == 3) {
	/* Deflate extension - multiple codes */
	int nbits, i;
	nbits  = get_bits(block, 4) + 1;
	cs->ncodes = get_bits(block, nbits) + 1;
	cs->codes = (huffman_codes_t **)malloc(cs->ncodes*sizeof(*cs->codes));
	for (i = 0; i < cs->ncodes; i++) {
	    cs->codes[i] = restore_codes_single(block);
	}
    } else {
	fprintf(stderr, "restore_codes() only implemented for "
		"BTYPE == DYNAMIC HUFFMAN and INTERLACED HUFFMAN\n");
	return NULL;
    }

    cs->bit_num = block->bit;

    return cs;
}

/*
 * Given a multiple sets of huffman codes and a block of data this
 * compresses and returns the data block. It iterates around each set
 * of huffman codes in a cyclic fashion encoding each byte with the
 * appropriate huffman codes.
 *
 * Returns: 0 on success
 *          -1 on failure
 */
int huffman_multi_encode(block_t *blk, huffman_codeset_t *cs,
			 int code_set, unsigned char *data, int len) {
    int i, nc;
    huffman_code_t *lookup;
    huffman_codes_t **c;

    if (!cs) {
	/* No codes known, so derive our own */
	fprintf(stderr, "FIXME: use generate_code_set() to build our own codes here\n");
	return -1;
    } else {
	c  = cs->codes;
	nc = cs->ncodes;
    }

    /*
     * The maximum size to encode len bytes is < 9 bits per symbol
     * (not quite 8 due to an EOF symbol) plus the overhead of the bit-length
     * tree. That in turn, with alternating 8/9 bit-lengths would max out
     * as 258*8 + 5+5+4 + 19*3 + 258*5 bits (429 bytes), but in practice
     * I'm not even sure if it's possible to construct such a set of code
     * lengths that would compress that poor.
     *
     * This of course assumes we're using appropriate compression codes.
     * Given a user may give a completely inappropriate code we have to
     * assume every symbol is actually 15 bits instead of < 9 on average.
     *
     * We ensure blk here is large enough for the worst case scenario so we
     * don't incur overheads in store_bits().
     */
    if (blk->alloc <= 429 + 2*len + 2 + blk->byte) {
	blk->alloc  = 429 + 2*len + 2 + blk->byte;
	blk->data = realloc(blk->data, blk->alloc);
	if (!blk->data)
	    return -1;
    }

    /*
     * Splitting this special case out is worth it as it's approx 30% faster.
     * Also note that the nc > 1 case is faster with a separate counter and
     * test than using modulus (by a factor of 2). It could be sped up further
     * for powers of 2 using bitwise AND, but the difference is not huge.
     */
    if (nc == 1) {
	lookup = c[0]->lookup;
	for (i = 0; i < len; i++) {
	    assert(lookup[data[i]].nbits > 0);
	    store_bits(blk, lookup[data[i]].code,
		       lookup[data[i]].nbits);
	}
    } else {
	int count = 0;
	for (i = 0; i < len; i++) {
	    lookup = c[count]->lookup;
	    assert(lookup[data[i]].nbits > 0);
	    store_bits(blk, lookup[data[i]].code,
		       lookup[data[i]].nbits);
	    if (++count == nc)
		count = 0;
	}
    }

    lookup = c[i%nc]->lookup;
    store_bits(blk, lookup[SYM_EOF].code, lookup[SYM_EOF].nbits);

    assert(blk->alloc > blk->byte);

    /* We probably massively overallocated, so return some of it back */
    blk->data = realloc(blk->data, blk->byte+1);
    blk->alloc = blk->byte+1;

    return 0;
}

/*
 * The opposite of huffman_encode().
 * Decode a huffman stream from 'block' using huffman codes 'c'.
 *
 * Returns: allocated block_t pointer on success
 *          NULL on failure.
 *
 * Method 1
 * --------
 *
 * At any node in our tree we can precompute a lookup table so that upon
 * reading the next 'k' bits we know the new node we'd end up in and what
 * symbols to export.
 * Then decoding simply works in fixed sets of k bits at a time.
 *
 * We use k=4 for efficient table space (they fit neatly in cache) and ease
 * of decoding 4-bits at a time. k=8 is about 20% faster as reading the input
 * byte by byte is easy, but the setup time is substantially longer
 * (16x at a guess) and the lookup tables no longer fit in the L1 cache.
 */
block_t *huffman_decode(block_t *in, huffman_codes_t *c) {
    block_t *out;
    htree_t t[513]; /* number of internal nodes */
    int i, j, n;
    int new_node, node_num;
    h_jump4_t J4[513][16];
    unsigned char *cp;

    if (NULL == (out = block_create(NULL, 8*in->alloc+8))) {
	block_destroy(in, 0);
	return NULL;
    }

    /* Construct the tree from the codes */
    new_node = 1;
    t[0].l[0] = t[0].l[1] = -1;
    t[0].c[0] = t[0].c[1] = 0;
    for (i = 0; i < c->ncodes; i++) {
	int n = 0;
	unsigned int v = c->codes[i].code;

	for (j = 0; j < c->codes[i].nbits-1; j++) {
	    int b = v & 1;
	    if (t[n].c[b]) {
		n = t[n].c[b];
	    } else {
		n = (t[n].c[b] = new_node++);
		t[n].c[0] = t[n].c[1] = 0;
		t[n].l[0] = t[n].l[1] = -1;
	    }
	    v >>= 1;
	}
	/* last bit */
	t[n].l[v & 1] = c->codes[i].symbol;
    }

    /* Build the 16 wide lookup table per node */
    for (n = 0; n < new_node; n++) {
	for (j = 0; j < 16; j++) {
	    unsigned int v = j;
	    int n2 = n;
	    h_jump4_t *hj = &J4[n][j];
	    hj->nsymbols = 0;
	    hj->top_bit = 0;
	    for (i = 0; i < 4; i++) {
		int b = v & 1;
		if (t[n2].l[b] >= 0) {
		    hj->symbol[hj->nsymbols++] = t[n2].l[b];
		    if (t[n2].l[b] == SYM_EOF)
			if (!hj->top_bit)
			    hj->top_bit |= 1 << (hj->nsymbols-1);
		}
		n2 = t[n2].c[b];
		v >>= 1;
	    }
	    hj->jump = n2;
	}
    }

#if 0
    /* Debug output */
    for (n = 0; n < new_node; n++) {
	printf("Node %d, c[]={%d,%d}, l[]={%d,%d}\n",
	       n, t[n].c[0], t[n].c[1], t[n].l[0], t[n].l[1]);
	for (i = 0; i < 256; i++) {
	    printf("\t%02x %s =>%02d, ", i, print_8rev(i), J4[n][i].jump);
	    for (k = 0; k < J4[n][i].nsymbols; k++) {
		if (isprint(J4[n][i].symbol[k]))
		    printf(" '%c'", J4[n][i].symbol[k]);
		else
		    printf(" %03d", J4[n][i].symbol[k]);
	    }
	    printf("\n");
	}
    }
#endif


    /*
     * Decoding - part 1
     * We're part way through a byte, so decode bit by bit up to the next
     * whole byte and then we start the fast decoding section.
     */
    cp = out->data;
    node_num = 0;
    while (in->bit != 0) {
	int b = get_bits(in, 1);
	if (t[node_num].l[b] != -1) {
	    if (t[node_num].l[b] != SYM_EOF) {
		*cp++ = t[node_num].l[b];
	    } else {
		out->byte = cp - out->data;
		return out;
	    }
	}
	node_num = t[node_num].c[b];
    }

    /*
     * Decoding - part 2
     *
     * We handle data nibble by nibble using the nibble to get an
     * h_jump4_t lookup from the J4[] table.
     * If top_bit is clear then we know we have no funny business (SYM_EOF)
     * so we use a fast decoding technique, otherwise we have to do a slower
     * loop with a check.
     */
    {
	int last_node = 0;
	unsigned char *last_cp = cp;
	h_jump4_t *x = &J4[node_num][in->data[in->byte] & 0x0f];
	int l = x->nsymbols;
	int b;

	/*
	 * This is the tight loop, so we over-optimise here by ignoring EOF
	 * and relying on knowing the length of the input data stream.
	 * This allows us to ignore the 9-bit data and only operate on
	 * the basic 0-255 symbols, glossing over the minor issue that EOF
	 * will look like an ordinary symbol.
	 */
	for (i = in->byte; i < in->alloc; i++) {
	    last_cp = cp;
	    last_node = node_num;

	    x = &J4[node_num][in->data[i] & 0x0f];
	    l = x->nsymbols;

	    for (j = 0; j < l; j++) {	
		*cp++ = x->symbol[j];
	    }
	    node_num = x->jump;

	    if (x->top_bit)
		break;

	    x = &J4[node_num][(in->data[i] >> 4) & 0x0f];
	    l = x->nsymbols;
	     
	    for (j = 0; j < l; j++) {	
		*cp++ = x->symbol[j];
	    }
	    node_num = x->jump;

	    if (x->top_bit)
		break;
	}

    
	/*
	 * Decoding - part 3
	 *
	 * The above optimisation has unfortunately added EOF to our data
	 * along with whatever else was packed in the last byte after the
	 * EOF symbol. So we rewind one byte and finish off decoding
	 * the slow way - walking the tree.
	 */
	cp = last_cp;
	node_num = last_node;
	in->byte = i;
	in->bit = 0;
	while (-1 != (b = get_bits(in, 1))) {
	    if (t[node_num].l[b] != -1) {
		if (t[node_num].l[b] != SYM_EOF) {
		    *cp++ = t[node_num].l[b];
		} else {
		    out->byte = cp - out->data;
		    return out;
		}
	    }
	    node_num = t[node_num].c[b];
	}
     }


    /* We shouldn't reach here */
    return NULL;
}

int init_decode_tables(huffman_codeset_t *cs) {
    int nnodes, i, j, n, nc;
    huffman_codes_t **c;
    int new_node, rec;
    h_jump4_t (*J4)[16] = NULL;
    htree_t *t;
    
    c = cs->codes;
    nc = cs->ncodes;

    /* Allocate memory for internal nodes (nsyms-1 for each code set) */
    for (nnodes = i = 0; i < nc; i++) {
	nnodes += c[i]->ncodes-1;
    }

    if (NULL == (t = (htree_t *)malloc(nnodes * sizeof(*t))))
	goto error;

    if (NULL == (J4 = (h_jump4_t (*)[16])malloc(nnodes * sizeof(*J4))))
	goto error;

    /*
     * Construct the tree from the codes.
     * We have one tree for all 'nc' huffman codes with each tree pointing
     * to the root of the next one (or first) tree whenever we emit a
     * symbol.
     *
     * This then effectively means the decoding step is identical to the
     * single huffman code function.
     */
    new_node = 0;
    for (rec = 0; rec < nc; rec++) {
	int root = new_node++;
	int next_root = rec == nc-1
	    ? 0
	    : root + c[rec]->ncodes-1;
	    
	t[root].l[0] = t[root].l[1] = -1;
	t[root].c[0] = t[root].c[1] = next_root;
	for (i = 0; i < c[rec]->ncodes; i++) {
	    int n = root;
	    unsigned int v = c[rec]->codes[i].code;

	    for (j = 0; j < c[rec]->codes[i].nbits-1; j++) {
		int b = v & 1;
		if (t[n].c[b] != next_root) {
		    n = t[n].c[b];
		} else {
		    n = (t[n].c[b] = new_node++);
		    t[n].c[0] = t[n].c[1] = next_root;
		    t[n].l[0] = t[n].l[1] = -1;
		}
		v >>= 1;
	    }
	    /* last bit */
	    t[n].l[v & 1] = c[rec]->codes[i].symbol;
	}
    }

    /*
    for (i = 0; i < new_node; i++) {
	printf("t[%d] = {(%d,%d), (%d,%d)}\n",
	       i,
	       t[i].l[0], t[i].l[1],
	       t[i].c[0], t[i].c[1]);
    }
    */

    /* Build the 16 wide lookup table per node */
    for (n = 0; n < new_node; n++) {
	for (j = 0; j < 16; j++) {
	    unsigned int v = j;
	    int n2 = n;
	    h_jump4_t *hj = &J4[n][j];
	    hj->nsymbols = 0;
	    hj->top_bit = 0;
	    for (i = 0; i < 4; i++) {
		int b = v & 1;
		if (t[n2].l[b] >= 0) {
		    hj->symbol[hj->nsymbols++] = t[n2].l[b];
		    if (t[n2].l[b] == SYM_EOF)
			if (!hj->top_bit)
			    hj->top_bit |= 1 << (hj->nsymbols-1);
		}
		n2 = t[n2].c[b];
		v >>= 1;
	    }
	    hj->jump = n2;
	    /*
	    printf("J4[%d][%d] = {'%.*s', %d}\n",
		   n, j, hj->nsymbols, hj->symbol, n2);
	    */
	}
    }

    cs->decode_t = t;
    cs->decode_J4 = J4;

    return 0;

 error:
    if (t)
	free(t);

    if (J4)
	free(J4);

    cs->decode_t = NULL;
    cs->decode_J4 = NULL;

    return -1;
}

/*
 * The opposite of huffman_encode().
 * Decode a huffman stream from 'block' using huffman codes 'c'.
 *
 * Returns: allocated block_t pointer on success
 *          NULL on failure.
 *
 * Method 1
 * --------
 *
 * At any node in our tree we can precompute a lookup table so that upon
 * reading the next 'k' bits we know the new node we'd end up in and what
 * symbols to export.
 * Then decoding simply works in fixed sets of k bits at a time.
 *
 * We use k=4 for efficient table space (they fit neatly in cache) and ease
 * of decoding 4-bits at a time. k=8 is about 20% faster as reading the input
 * byte by byte is easy, but the setup time is substantially longer
 * (16x at a guess) and the lookup tables no longer fit in the L1 cache.
 *
 * NB: This version also handles multiple interleaved huffman codes as
 * this support doesn't really slow down the decoding process.
 */
block_t *huffman_multi_decode(block_t *in, huffman_codeset_t *cs) {
    block_t *out = NULL;
    int i, j;
    int node_num;
    unsigned char *cp;
    huffman_codes_t **c;
    int nc;
    h_jump4_t (*J4)[16];
    htree_t *t;

    if (!cs)
	return NULL;
    c = cs->codes;
    nc = cs->ncodes;

    /* Ensure precomputed lookup tables exist */
    if (!cs->decode_t || !cs->decode_J4)
	if (-1 == init_decode_tables(cs))
	    return NULL;

    t  = cs->decode_t;
    J4 = cs->decode_J4;

    if (NULL == (out = block_create(NULL, 9*(in->alloc+1)))) {
	goto error;
    }

    /*
     * Decoding - part 1
     * We're part way through a byte, so decode bit by bit up to the next
     * whole byte and then we start the fast decoding section.
     */
    cp = out->data;
    node_num = 0;
    while (in->bit != 0) {
	int b = get_bits(in, 1);
	htree_t *t2 = &t[node_num];
	if (t2->l[b] != -1) {
	    if (t2->l[b] != SYM_EOF) {
		*cp++ = t2->l[b];
	    } else {
		out->byte = cp - out->data;
		goto success;
	    }
	}
	node_num = t2->c[b];
    }

    /*
     * Decoding - part 2
     *
     * We now handle data nibble by nibble using the nibble to get an
     * h_jump4_t lookup from the J4[] table.
     * If top_bit is clear then we know we have no funny business (SYM_EOF)
     * so we use a fast decoding technique, otherwise we have to do a slower
     * loop with a check.
     */
    {
	int last_node = node_num;
	unsigned char *last_cp = cp;
	h_jump4_t *x = &J4[node_num][in->data[in->byte] & 0x0f];
	int l = x->nsymbols;
	int b;

	/*
	 * This is the tight loop, so we over-optimise here by ignoring EOF
	 * and relying on knowing the length of the input data stream.
	 * This allows us to ignore the 9-bit data and only operate on
	 * the basic 0-255 symbols, glossing over the minor issue that EOF
	 * will look like an ordinary symbol.
	 */
	for (i = in->byte; i < in->alloc; i++) {
	    last_cp = cp;
	    last_node = node_num;

	    x = &J4[node_num][in->data[i] & 0x0f];
	    l = x->nsymbols;

	    /* printf("val=%d\n", in->data[i] & 0x0f); */
	    for (j = 0; j < l; j++) {
		*cp++ = x->symbol[j];
	    }
	    node_num = x->jump;

	    if (x->top_bit)
		break;

	    x = &J4[node_num][(in->data[i] >> 4) & 0x0f];
	    l = x->nsymbols;
	     
	    for (j = 0; j < l; j++) {	
		*cp++ = x->symbol[j];
	    }
	    node_num = x->jump;

	    if (x->top_bit)
		break;
	}

    
	/*
	 * Decoding - part 3
	 *
	 * The above optimisation has unfortunately added EOF to our data
	 * along with whatever else was packed in the last byte after the
	 * EOF symbol. So we rewind one byte and finish off decoding
	 * the slow way - walking the tree.
	 */
	cp = last_cp;
	node_num = last_node;
	in->byte = i;
	in->bit = 0;
	while (-1 != (b = get_bits(in, 1))) {
	    htree_t *t2 = &t[node_num];
	    if (t2->l[b] != -1) {
		if (t2->l[b] != SYM_EOF) {
		    *cp++ = t2->l[b];
		} else {
		    out->byte = cp - out->data;
		    goto success;
		}
	    }
	    node_num = t2->c[b];
	}
     }

 success:
    return out;

 error:
    if (out)
	block_destroy(out, 0);

    return NULL;
}

#if 0
/* A simple to understand (but slow) version of the above function */
block_t *huffman_multi_decode(block_t *in, huffman_codes_t **c, int nc) {
    block_t *out;
    htree_t (*t)[513]; /* number of internal nodes */
    int i, j, n, rec;
    int new_node, node_num;
    unsigned char *cp;
    int bC; /* byte count */

    if (NULL == (t = (htree_t (*)[513])malloc(nc * sizeof(*t))))
	return NULL;

    if (NULL == (out = block_create(NULL, 8*in->alloc+8))) {
	block_destroy(in, 0);
	free(t);
	return NULL;
    }

    /* Construct the tree from the codes */
    for (rec = 0; rec < nc; rec++) {
	new_node = 1;
	t[rec][0].l[0] = t[rec][0].l[1] = -1;
	t[rec][0].c[0] = t[rec][0].c[1] = 0;
	for (i = 0; i < c[rec]->ncodes; i++) {
	    int n = 0;
	    unsigned int v = c[rec]->codes[i].code;

	    for (j = 0; j < c[rec]->codes[i].nbits-1; j++) {
		int b = v & 1;
		if (t[rec][n].c[b]) {
		    n = t[rec][n].c[b];
		} else {
		    n = (t[rec][n].c[b] = new_node++);
		    t[rec][n].c[0] = t[rec][n].c[1] = 0;
		    t[rec][n].l[0] = t[rec][n].l[1] = -1;
		}
		v >>= 1;
	    }
	    /* last bit */
	    t[rec][n].l[v & 1] = c[rec]->codes[i].symbol;
	}
    }

    /*
     * Decoding - the slow way. How to speed up multi-code decoding?
     */
    cp = out->data;
    node_num = 0;
    bC = 0;
    while (in->byte < in->alloc) {
	htree_t *t2;
	int b = get_bits(in, 1);
	t2 = &t[bC%nc][node_num];
	if (t2->l[b] != -1) {
	    if (t2->l[b] != SYM_EOF) {
		*cp++ = t2->l[b];
		bC++;
		node_num = 0;
	    } else {
		out->byte = cp - out->data;
		free(t);
		return out;
	    }
	} else {
	    node_num = t2->c[b];
	}
    }

    /* We shouldn't reach here */
    free(t);
    return NULL;
}
#endif

/*
 * A slow version of the above huffman_decode function. This is designed as
 * a piecemeal decoder for purposes of restoring the huffman codes themselves.
 * NB: this only works for code lengths small enough to keep inside the
 * htab[] dimensions - IT DOES NOT CHECK THIS.
 *
 * Returns the next symbol
 *         -1 for failure
 */
int next_symbol(block_t *in, int *htab) {
    int b, v = 0, c = 1;
    while ((b = get_bits(in, 1)) != -1) {
	v = (v<<1) | b | (c <<= 1);

	if (htab[v] != -1)
	    return htab[v];
    }

    return -1;
}


/*
 * ---------------------------------------------------------------------------
 * Debug code. This turns the library into a stand-alone program for
 * easy debugging.x
 */
static void output_code_set(FILE *fp, huffman_codes_t *cds) {
    int i, j;
    int nbits_in = 0, nbits_out = 0;
    huffman_code_t *codes = cds->codes;
    int ncodes = cds->ncodes;

    fprintf(fp, "static huffman_code_t codes_FIXME[] = {\n");
    for (i = j = 0; i < ncodes; i++) {
	nbits_out += codes[i].nbits * codes[i].freq;
	nbits_in  += 8*codes[i].freq;
	if (j == 0)
	    fprintf(fp, "    ");
	if (codes[i].symbol == SYM_EOF) {
	    fprintf(fp, "{SYM_EOF,%3d}, ", codes[i].nbits);
	    j = 10;
	} else {
	    if (isalnum(codes[i].symbol)) {
		fprintf(fp, "{'%c',%3d}, ", codes[i].symbol, codes[i].nbits);
	    } else {
		fprintf(fp, "{%3d,%3d}, ", codes[i].symbol, codes[i].nbits);
	    }
	}
	j++;
	
	if (j >= 6) {
	    fputc('\n', fp);
	    j = 0;
	}
    }
    if (j)
	fputc('\n', fp);
    fprintf(fp, "};\n");
    fprintf(fp, "/* Expected compression to %f of input */\n",
	    (double)nbits_out/nbits_in);
}

static void output_code_set2(FILE *fp, huffman_codes_t *cds) {
    int i;
    huffman_code_t *codes = cds->codes;
    int ncodes = cds->ncodes;

    fprintf(fp, "huffman_code_t = {\n");
    for (i = 0; i < ncodes; i++) {
	fprintf(fp, "\t%d:\t%3d %c %2d %04x %d\n",
		i,codes[i].symbol, codes[i].symbol,
		codes[i].nbits, codes[i].code,
		codes[i].freq);
    }
    fprintf(fp, "};\n");
}


/*
 * --------------------------------------------------------------------------
 * A test main() to create an application capable of compressing and
 * uncompressing stdin.
 */

#ifdef TEST_MAIN
#include <fcntl.h>
/* #include <zlib.h> */

/*
 * Slurps the entirety of stdin into a malloced buffer and returns a pointer
 * to it.
 *
 * Returns: malloced buffer on success, *lenp equal to length
 *          NULL on failure
 */
static unsigned char *load(int *lenp, char *fn) {
    unsigned char *data = NULL;
    int dsize = 0;
    int dcurr = 0, len;
    int fd = 0;

    if (fn) {
	if (-1 == (fd = open(fn, O_RDONLY, 0))) {
	    perror(fn);
	    return NULL;
	}
    }

    do {
	if (dsize - dcurr < 8192) {
	    dsize = dsize ? dsize * 2 : 8192;
	    if (NULL == (data = realloc(data, dsize))) {
		if (fd)
		    close(fd);
		return NULL;
	    }
	}

	len = read(fd, data + dcurr, 8192);
	if (len > 0)
	    dcurr += len;
    } while (len > 0);

    if (len == -1) {
	perror("read");
	if (fd)
	    close(fd);
	return NULL;
    }

    if (fd)
	close(fd);

    *lenp = dcurr;
    return data;
}

/*
 * Returns 0 for success
 *        -1 for failure.
 */
int decode_main(unsigned char *data, int len, int code_set) {
    huffman_codeset_t *cs = NULL;
    block_t *blk_in;
    unsigned char *out;
    int err, out_len, i;
    int bfinal;

    blk_in = block_create(NULL, 1000 + len);

    if (code_set != 0) {
	cs = generate_code_set(code_set, 1,  /* no. codes */
			       NULL, 0,      /* data + size */
			       1,	     /* eof */
			       MAX_CODE_LEN,
			       0);	     /* all_codes */
	store_codes(blk_in, cs, 1);
    }

    if (blk_in->bit != 0) {
	blk_in->data[blk_in->byte] |= data[0];
	memcpy(&blk_in->data[blk_in->byte+1], data+1, len-1);
    } else {
	memcpy(&blk_in->data[blk_in->byte], data, len);
    }

    /* Do the decoding */
    do {
	block_t *out;
	if (!cs)
	    cs = restore_codes(blk_in, &bfinal);
	out = huffman_multi_decode(blk_in, cs);
	write(1, out->data, out->byte);
	block_destroy(out, 0);
	huffman_codeset_destroy(cs);
	cs = NULL;
    } while (!bfinal);

    block_destroy(blk_in, 0);

    return 0;
}

/*
 * Returns 0 for success
 *        -1 for failure.
 */
int encode_main(unsigned char *data, int len, int code_set, int rec_size,
		int blk_size, int dump_tree, int exit_after_tree) {
    /* Encoding */
    unsigned char *d2 = data;
    block_t *blk;
    huffman_codeset_t *cs;

    blk = block_create(NULL, 8192);
    fprintf(stderr, "Input %d bytes\n", len);

    do {
	int l2 = len > blk_size ? blk_size : len;
	int i;

	if (code_set != 0)
	    l2 = len; /* predefined code-sets have final-block bit set */
	cs = generate_code_set(code_set, rec_size,
			       d2, l2, /* Data and length */
			       1,      /* eof */
			       MAX_CODE_LEN,
			       0);     /* all codes */
	if (!cs)
	    return -1;

	if (dump_tree) {
	    for (i = 0; i < rec_size; i++) {
		printf("==Sub-set %d==\n", i);
		output_code_set(stdout, cs->codes[i]);
		/* output_code_set2(stdout, cs->codes[i]); */
	    }
	    if (exit_after_tree)
		return 0;
	}

	store_codes(blk, cs, l2 == len);
	if (code_set != 0) {
	    blk->data[blk->byte = 0] = 0;  /* starting bit no. preseved */
	} else {
	    /*
	      fprintf(stderr, "codes finished at %d bytes, %d bits\n",
		    blk->byte, blk->bit);
	    */
	}
	if (exit_after_tree) {
	    write(1, blk->data, blk->byte + (blk->bit != 0));
	    return 0;
	}

	huffman_multi_encode(blk, cs, code_set, d2, l2);

	huffman_codeset_destroy(cs);
	len -= l2;
	d2  += l2;

    } while (len > 0);
    write(1, blk->data, blk->byte + (blk->bit != 0));

    fprintf(stderr, "Output %ld bytes\n", blk->byte + (blk->bit != 0));
    block_destroy(blk, 0);

    return 0;
}

int main(int argc, char **argv) {
    unsigned char *data;
    int decode = 0;
    int dump_tree = 0;
    int exit_after_tree = 0;
    int code_set = CODE_INLINE;
    int blk_size = 0x7fff;
    int rec_size = 1;
    int c, r, len;
    char *fn = NULL;

    while ((c = getopt(argc, argv, "c:detxl:b:hr:i:")) != -1) {
	switch (c) {
	case 'b':
	    blk_size = atoi(optarg);
	    break;

	case 'c':
	    code_set = atoi(optarg);
	    break;

	case 'r':
	    rec_size = atoi(optarg);
	    break;

	case 'd':
	    decode = 1;
	    break;

	case 'e':
	    decode = 0;
	    break;

	case 't':
	    dump_tree = 1;
	    break;

	case 'x':
	    exit_after_tree = 1;
	    break;
	    
	case 'i':
	    fn = optarg;
	    break;

	default:
	    fprintf(stderr, "Usage: huffman_static [options] < stdin > stdout\n");
	    fprintf(stderr, "    Decoding options\n");
	    fprintf(stderr, "        -d\tdecode flag\n");
	    fprintf(stderr, "    Encoding options\n");
	    fprintf(stderr, "        -e\tencode flag\n");
	    fprintf(stderr, "        -b size\tspecify the block-size\n");
	    fprintf(stderr, "        -c code\tspecify code-set. 0 => inline\n");
	    fprintf(stderr, "        -t\tpretty-print the code-set used\n");
	    fprintf(stderr, "        -x\texit after outputting code-set\n");
	    exit(c != 'h');
	}
    }

    data = load(&len, fn);

    r = decode
	? decode_main(data, len, code_set)
	: encode_main(data, len, code_set, rec_size,
		      blk_size, dump_tree, exit_after_tree);

    free(data);

    return r == 0 ? 0 : 1;
}
#endif