Home | History | Annotate | Download | only in examples
      1 /* enough.c -- determine the maximum size of inflate's Huffman code tables over
      2  * all possible valid and complete Huffman codes, subject to a length limit.
      3  * Copyright (C) 2007, 2008, 2012 Mark Adler
      4  * Version 1.4  18 August 2012  Mark Adler
      5  */
      6 
      7 /* Version history:
      8    1.0   3 Jan 2007  First version (derived from codecount.c version 1.4)
      9    1.1   4 Jan 2007  Use faster incremental table usage computation
     10                      Prune examine() search on previously visited states
     11    1.2   5 Jan 2007  Comments clean up
     12                      As inflate does, decrease root for short codes
     13                      Refuse cases where inflate would increase root
     14    1.3  17 Feb 2008  Add argument for initial root table size
     15                      Fix bug for initial root table size == max - 1
     16                      Use a macro to compute the history index
     17    1.4  18 Aug 2012  Avoid shifts more than bits in type (caused endless loop!)
     18                      Clean up comparisons of different types
     19                      Clean up code indentation
     20  */
     21 
     22 /*
     23    Examine all possible Huffman codes for a given number of symbols and a
     24    maximum code length in bits to determine the maximum table size for zilb's
     25    inflate.  Only complete Huffman codes are counted.
     26 
     27    Two codes are considered distinct if the vectors of the number of codes per
     28    length are not identical.  So permutations of the symbol assignments result
     29    in the same code for the counting, as do permutations of the assignments of
     30    the bit values to the codes (i.e. only canonical codes are counted).
     31 
     32    We build a code from shorter to longer lengths, determining how many symbols
     33    are coded at each length.  At each step, we have how many symbols remain to
     34    be coded, what the last code length used was, and how many bit patterns of
     35    that length remain unused. Then we add one to the code length and double the
     36    number of unused patterns to graduate to the next code length.  We then
     37    assign all portions of the remaining symbols to that code length that
     38    preserve the properties of a correct and eventually complete code.  Those
     39    properties are: we cannot use more bit patterns than are available; and when
     40    all the symbols are used, there are exactly zero possible bit patterns
     41    remaining.
     42 
     43    The inflate Huffman decoding algorithm uses two-level lookup tables for
     44    speed.  There is a single first-level table to decode codes up to root bits
     45    in length (root == 9 in the current inflate implementation).  The table
     46    has 1 << root entries and is indexed by the next root bits of input.  Codes
     47    shorter than root bits have replicated table entries, so that the correct
     48    entry is pointed to regardless of the bits that follow the short code.  If
     49    the code is longer than root bits, then the table entry points to a second-
     50    level table.  The size of that table is determined by the longest code with
     51    that root-bit prefix.  If that longest code has length len, then the table
     52    has size 1 << (len - root), to index the remaining bits in that set of
     53    codes.  Each subsequent root-bit prefix then has its own sub-table.  The
     54    total number of table entries required by the code is calculated
     55    incrementally as the number of codes at each bit length is populated.  When
     56    all of the codes are shorter than root bits, then root is reduced to the
     57    longest code length, resulting in a single, smaller, one-level table.
     58 
     59    The inflate algorithm also provides for small values of root (relative to
     60    the log2 of the number of symbols), where the shortest code has more bits
     61    than root.  In that case, root is increased to the length of the shortest
     62    code.  This program, by design, does not handle that case, so it is verified
     63    that the number of symbols is less than 2^(root + 1).
     64 
     65    In order to speed up the examination (by about ten orders of magnitude for
     66    the default arguments), the intermediate states in the build-up of a code
     67    are remembered and previously visited branches are pruned.  The memory
     68    required for this will increase rapidly with the total number of symbols and
     69    the maximum code length in bits.  However this is a very small price to pay
     70    for the vast speedup.
     71 
     72    First, all of the possible Huffman codes are counted, and reachable
     73    intermediate states are noted by a non-zero count in a saved-results array.
     74    Second, the intermediate states that lead to (root + 1) bit or longer codes
     75    are used to look at all sub-codes from those junctures for their inflate
     76    memory usage.  (The amount of memory used is not affected by the number of
     77    codes of root bits or less in length.)  Third, the visited states in the
     78    construction of those sub-codes and the associated calculation of the table
     79    size is recalled in order to avoid recalculating from the same juncture.
     80    Beginning the code examination at (root + 1) bit codes, which is enabled by
     81    identifying the reachable nodes, accounts for about six of the orders of
     82    magnitude of improvement for the default arguments.  About another four
     83    orders of magnitude come from not revisiting previous states.  Out of
     84    approximately 2x10^16 possible Huffman codes, only about 2x10^6 sub-codes
     85    need to be examined to cover all of the possible table memory usage cases
     86    for the default arguments of 286 symbols limited to 15-bit codes.
     87 
     88    Note that an unsigned long long type is used for counting.  It is quite easy
     89    to exceed the capacity of an eight-byte integer with a large number of
     90    symbols and a large maximum code length, so multiple-precision arithmetic
     91    would need to replace the unsigned long long arithmetic in that case.  This
     92    program will abort if an overflow occurs.  The big_t type identifies where
     93    the counting takes place.
     94 
     95    An unsigned long long type is also used for calculating the number of
     96    possible codes remaining at the maximum length.  This limits the maximum
     97    code length to the number of bits in a long long minus the number of bits
     98    needed to represent the symbols in a flat code.  The code_t type identifies
     99    where the bit pattern counting takes place.
    100  */
    101 
    102 #include <stdio.h>
    103 #include <stdlib.h>
    104 #include <string.h>
    105 #include <assert.h>
    106 
    107 #define local static
    108 
    109 /* special data types */
    110 typedef unsigned long long big_t;   /* type for code counting */
    111 typedef unsigned long long code_t;  /* type for bit pattern counting */
    112 struct tab {                        /* type for been here check */
    113     size_t len;         /* length of bit vector in char's */
    114     char *vec;          /* allocated bit vector */
    115 };
    116 
    117 /* The array for saving results, num[], is indexed with this triplet:
    118 
    119       syms: number of symbols remaining to code
    120       left: number of available bit patterns at length len
    121       len: number of bits in the codes currently being assigned
    122 
    123    Those indices are constrained thusly when saving results:
    124 
    125       syms: 3..totsym (totsym == total symbols to code)
    126       left: 2..syms - 1, but only the evens (so syms == 8 -> 2, 4, 6)
    127       len: 1..max - 1 (max == maximum code length in bits)
    128 
    129    syms == 2 is not saved since that immediately leads to a single code.  left
    130    must be even, since it represents the number of available bit patterns at
    131    the current length, which is double the number at the previous length.
    132    left ends at syms-1 since left == syms immediately results in a single code.
    133    (left > sym is not allowed since that would result in an incomplete code.)
    134    len is less than max, since the code completes immediately when len == max.
    135 
    136    The offset into the array is calculated for the three indices with the
    137    first one (syms) being outermost, and the last one (len) being innermost.
    138    We build the array with length max-1 lists for the len index, with syms-3
    139    of those for each symbol.  There are totsym-2 of those, with each one
    140    varying in length as a function of sym.  See the calculation of index in
    141    count() for the index, and the calculation of size in main() for the size
    142    of the array.
    143 
    144    For the deflate example of 286 symbols limited to 15-bit codes, the array
    145    has 284,284 entries, taking up 2.17 MB for an 8-byte big_t.  More than
    146    half of the space allocated for saved results is actually used -- not all
    147    possible triplets are reached in the generation of valid Huffman codes.
    148  */
    149 
    150 /* The array for tracking visited states, done[], is itself indexed identically
    151    to the num[] array as described above for the (syms, left, len) triplet.
    152    Each element in the array is further indexed by the (mem, rem) doublet,
    153    where mem is the amount of inflate table space used so far, and rem is the
    154    remaining unused entries in the current inflate sub-table.  Each indexed
    155    element is simply one bit indicating whether the state has been visited or
    156    not.  Since the ranges for mem and rem are not known a priori, each bit
    157    vector is of a variable size, and grows as needed to accommodate the visited
    158    states.  mem and rem are used to calculate a single index in a triangular
    159    array.  Since the range of mem is expected in the default case to be about
    160    ten times larger than the range of rem, the array is skewed to reduce the
    161    memory usage, with eight times the range for mem than for rem.  See the
    162    calculations for offset and bit in beenhere() for the details.
    163 
    164    For the deflate example of 286 symbols limited to 15-bit codes, the bit
    165    vectors grow to total approximately 21 MB, in addition to the 4.3 MB done[]
    166    array itself.
    167  */
    168 
    169 /* Globals to avoid propagating constants or constant pointers recursively */
    170 local int max;          /* maximum allowed bit length for the codes */
    171 local int root;         /* size of base code table in bits */
    172 local int large;        /* largest code table so far */
    173 local size_t size;      /* number of elements in num and done */
    174 local int *code;        /* number of symbols assigned to each bit length */
    175 local big_t *num;       /* saved results array for code counting */
    176 local struct tab *done; /* states already evaluated array */
    177 
    178 /* Index function for num[] and done[] */
    179 #define INDEX(i,j,k) (((size_t)((i-1)>>1)*((i-2)>>1)+(j>>1)-1)*(max-1)+k-1)
    180 
    181 /* Free allocated space.  Uses globals code, num, and done. */
    182 local void cleanup(void)
    183 {
    184     size_t n;
    185 
    186     if (done != NULL) {
    187         for (n = 0; n < size; n++)
    188             if (done[n].len)
    189                 free(done[n].vec);
    190         free(done);
    191     }
    192     if (num != NULL)
    193         free(num);
    194     if (code != NULL)
    195         free(code);
    196 }
    197 
    198 /* Return the number of possible Huffman codes using bit patterns of lengths
    199    len through max inclusive, coding syms symbols, with left bit patterns of
    200    length len unused -- return -1 if there is an overflow in the counting.
    201    Keep a record of previous results in num to prevent repeating the same
    202    calculation.  Uses the globals max and num. */
    203 local big_t count(int syms, int len, int left)
    204 {
    205     big_t sum;          /* number of possible codes from this juncture */
    206     big_t got;          /* value returned from count() */
    207     int least;          /* least number of syms to use at this juncture */
    208     int most;           /* most number of syms to use at this juncture */
    209     int use;            /* number of bit patterns to use in next call */
    210     size_t index;       /* index of this case in *num */
    211 
    212     /* see if only one possible code */
    213     if (syms == left)
    214         return 1;
    215 
    216     /* note and verify the expected state */
    217     assert(syms > left && left > 0 && len < max);
    218 
    219     /* see if we've done this one already */
    220     index = INDEX(syms, left, len);
    221     got = num[index];
    222     if (got)
    223         return got;         /* we have -- return the saved result */
    224 
    225     /* we need to use at least this many bit patterns so that the code won't be
    226        incomplete at the next length (more bit patterns than symbols) */
    227     least = (left << 1) - syms;
    228     if (least < 0)
    229         least = 0;
    230 
    231     /* we can use at most this many bit patterns, lest there not be enough
    232        available for the remaining symbols at the maximum length (if there were
    233        no limit to the code length, this would become: most = left - 1) */
    234     most = (((code_t)left << (max - len)) - syms) /
    235             (((code_t)1 << (max - len)) - 1);
    236 
    237     /* count all possible codes from this juncture and add them up */
    238     sum = 0;
    239     for (use = least; use <= most; use++) {
    240         got = count(syms - use, len + 1, (left - use) << 1);
    241         sum += got;
    242         if (got == (big_t)0 - 1 || sum < got)   /* overflow */
    243             return (big_t)0 - 1;
    244     }
    245 
    246     /* verify that all recursive calls are productive */
    247     assert(sum != 0);
    248 
    249     /* save the result and return it */
    250     num[index] = sum;
    251     return sum;
    252 }
    253 
    254 /* Return true if we've been here before, set to true if not.  Set a bit in a
    255    bit vector to indicate visiting this state.  Each (syms,len,left) state
    256    has a variable size bit vector indexed by (mem,rem).  The bit vector is
    257    lengthened if needed to allow setting the (mem,rem) bit. */
    258 local int beenhere(int syms, int len, int left, int mem, int rem)
    259 {
    260     size_t index;       /* index for this state's bit vector */
    261     size_t offset;      /* offset in this state's bit vector */
    262     int bit;            /* mask for this state's bit */
    263     size_t length;      /* length of the bit vector in bytes */
    264     char *vector;       /* new or enlarged bit vector */
    265 
    266     /* point to vector for (syms,left,len), bit in vector for (mem,rem) */
    267     index = INDEX(syms, left, len);
    268     mem -= 1 << root;
    269     offset = (mem >> 3) + rem;
    270     offset = ((offset * (offset + 1)) >> 1) + rem;
    271     bit = 1 << (mem & 7);
    272 
    273     /* see if we've been here */
    274     length = done[index].len;
    275     if (offset < length && (done[index].vec[offset] & bit) != 0)
    276         return 1;       /* done this! */
    277 
    278     /* we haven't been here before -- set the bit to show we have now */
    279 
    280     /* see if we need to lengthen the vector in order to set the bit */
    281     if (length <= offset) {
    282         /* if we have one already, enlarge it, zero out the appended space */
    283         if (length) {
    284             do {
    285                 length <<= 1;
    286             } while (length <= offset);
    287             vector = realloc(done[index].vec, length);
    288             if (vector != NULL)
    289                 memset(vector + done[index].len, 0, length - done[index].len);
    290         }
    291 
    292         /* otherwise we need to make a new vector and zero it out */
    293         else {
    294             length = 1 << (len - root);
    295             while (length <= offset)
    296                 length <<= 1;
    297             vector = calloc(length, sizeof(char));
    298         }
    299 
    300         /* in either case, bail if we can't get the memory */
    301         if (vector == NULL) {
    302             fputs("abort: unable to allocate enough memory\n", stderr);
    303             cleanup();
    304             exit(1);
    305         }
    306 
    307         /* install the new vector */
    308         done[index].len = length;
    309         done[index].vec = vector;
    310     }
    311 
    312     /* set the bit */
    313     done[index].vec[offset] |= bit;
    314     return 0;
    315 }
    316 
    317 /* Examine all possible codes from the given node (syms, len, left).  Compute
    318    the amount of memory required to build inflate's decoding tables, where the
    319    number of code structures used so far is mem, and the number remaining in
    320    the current sub-table is rem.  Uses the globals max, code, root, large, and
    321    done. */
    322 local void examine(int syms, int len, int left, int mem, int rem)
    323 {
    324     int least;          /* least number of syms to use at this juncture */
    325     int most;           /* most number of syms to use at this juncture */
    326     int use;            /* number of bit patterns to use in next call */
    327 
    328     /* see if we have a complete code */
    329     if (syms == left) {
    330         /* set the last code entry */
    331         code[len] = left;
    332 
    333         /* complete computation of memory used by this code */
    334         while (rem < left) {
    335             left -= rem;
    336             rem = 1 << (len - root);
    337             mem += rem;
    338         }
    339         assert(rem == left);
    340 
    341         /* if this is a new maximum, show the entries used and the sub-code */
    342         if (mem > large) {
    343             large = mem;
    344             printf("max %d: ", mem);
    345             for (use = root + 1; use <= max; use++)
    346                 if (code[use])
    347                     printf("%d[%d] ", code[use], use);
    348             putchar('\n');
    349             fflush(stdout);
    350         }
    351 
    352         /* remove entries as we drop back down in the recursion */
    353         code[len] = 0;
    354         return;
    355     }
    356 
    357     /* prune the tree if we can */
    358     if (beenhere(syms, len, left, mem, rem))
    359         return;
    360 
    361     /* we need to use at least this many bit patterns so that the code won't be
    362        incomplete at the next length (more bit patterns than symbols) */
    363     least = (left << 1) - syms;
    364     if (least < 0)
    365         least = 0;
    366 
    367     /* we can use at most this many bit patterns, lest there not be enough
    368        available for the remaining symbols at the maximum length (if there were
    369        no limit to the code length, this would become: most = left - 1) */
    370     most = (((code_t)left << (max - len)) - syms) /
    371             (((code_t)1 << (max - len)) - 1);
    372 
    373     /* occupy least table spaces, creating new sub-tables as needed */
    374     use = least;
    375     while (rem < use) {
    376         use -= rem;
    377         rem = 1 << (len - root);
    378         mem += rem;
    379     }
    380     rem -= use;
    381 
    382     /* examine codes from here, updating table space as we go */
    383     for (use = least; use <= most; use++) {
    384         code[len] = use;
    385         examine(syms - use, len + 1, (left - use) << 1,
    386                 mem + (rem ? 1 << (len - root) : 0), rem << 1);
    387         if (rem == 0) {
    388             rem = 1 << (len - root);
    389             mem += rem;
    390         }
    391         rem--;
    392     }
    393 
    394     /* remove entries as we drop back down in the recursion */
    395     code[len] = 0;
    396 }
    397 
    398 /* Look at all sub-codes starting with root + 1 bits.  Look at only the valid
    399    intermediate code states (syms, left, len).  For each completed code,
    400    calculate the amount of memory required by inflate to build the decoding
    401    tables. Find the maximum amount of memory required and show the code that
    402    requires that maximum.  Uses the globals max, root, and num. */
    403 local void enough(int syms)
    404 {
    405     int n;              /* number of remaing symbols for this node */
    406     int left;           /* number of unused bit patterns at this length */
    407     size_t index;       /* index of this case in *num */
    408 
    409     /* clear code */
    410     for (n = 0; n <= max; n++)
    411         code[n] = 0;
    412 
    413     /* look at all (root + 1) bit and longer codes */
    414     large = 1 << root;              /* base table */
    415     if (root < max)                 /* otherwise, there's only a base table */
    416         for (n = 3; n <= syms; n++)
    417             for (left = 2; left < n; left += 2)
    418             {
    419                 /* look at all reachable (root + 1) bit nodes, and the
    420                    resulting codes (complete at root + 2 or more) */
    421                 index = INDEX(n, left, root + 1);
    422                 if (root + 1 < max && num[index])       /* reachable node */
    423                     examine(n, root + 1, left, 1 << root, 0);
    424 
    425                 /* also look at root bit codes with completions at root + 1
    426                    bits (not saved in num, since complete), just in case */
    427                 if (num[index - 1] && n <= left << 1)
    428                     examine((n - left) << 1, root + 1, (n - left) << 1,
    429                             1 << root, 0);
    430             }
    431 
    432     /* done */
    433     printf("done: maximum of %d table entries\n", large);
    434 }
    435 
    436 /*
    437    Examine and show the total number of possible Huffman codes for a given
    438    maximum number of symbols, initial root table size, and maximum code length
    439    in bits -- those are the command arguments in that order.  The default
    440    values are 286, 9, and 15 respectively, for the deflate literal/length code.
    441    The possible codes are counted for each number of coded symbols from two to
    442    the maximum.  The counts for each of those and the total number of codes are
    443    shown.  The maximum number of inflate table entires is then calculated
    444    across all possible codes.  Each new maximum number of table entries and the
    445    associated sub-code (starting at root + 1 == 10 bits) is shown.
    446 
    447    To count and examine Huffman codes that are not length-limited, provide a
    448    maximum length equal to the number of symbols minus one.
    449 
    450    For the deflate literal/length code, use "enough".  For the deflate distance
    451    code, use "enough 30 6".
    452 
    453    This uses the %llu printf format to print big_t numbers, which assumes that
    454    big_t is an unsigned long long.  If the big_t type is changed (for example
    455    to a multiple precision type), the method of printing will also need to be
    456    updated.
    457  */
    458 int main(int argc, char **argv)
    459 {
    460     int syms;           /* total number of symbols to code */
    461     int n;              /* number of symbols to code for this run */
    462     big_t got;          /* return value of count() */
    463     big_t sum;          /* accumulated number of codes over n */
    464     code_t word;        /* for counting bits in code_t */
    465 
    466     /* set up globals for cleanup() */
    467     code = NULL;
    468     num = NULL;
    469     done = NULL;
    470 
    471     /* get arguments -- default to the deflate literal/length code */
    472     syms = 286;
    473     root = 9;
    474     max = 15;
    475     if (argc > 1) {
    476         syms = atoi(argv[1]);
    477         if (argc > 2) {
    478             root = atoi(argv[2]);
    479             if (argc > 3)
    480                 max = atoi(argv[3]);
    481         }
    482     }
    483     if (argc > 4 || syms < 2 || root < 1 || max < 1) {
    484         fputs("invalid arguments, need: [sym >= 2 [root >= 1 [max >= 1]]]\n",
    485               stderr);
    486         return 1;
    487     }
    488 
    489     /* if not restricting the code length, the longest is syms - 1 */
    490     if (max > syms - 1)
    491         max = syms - 1;
    492 
    493     /* determine the number of bits in a code_t */
    494     for (n = 0, word = 1; word; n++, word <<= 1)
    495         ;
    496 
    497     /* make sure that the calculation of most will not overflow */
    498     if (max > n || (code_t)(syms - 2) >= (((code_t)0 - 1) >> (max - 1))) {
    499         fputs("abort: code length too long for internal types\n", stderr);
    500         return 1;
    501     }
    502 
    503     /* reject impossible code requests */
    504     if ((code_t)(syms - 1) > ((code_t)1 << max) - 1) {
    505         fprintf(stderr, "%d symbols cannot be coded in %d bits\n",
    506                 syms, max);
    507         return 1;
    508     }
    509 
    510     /* allocate code vector */
    511     code = calloc(max + 1, sizeof(int));
    512     if (code == NULL) {
    513         fputs("abort: unable to allocate enough memory\n", stderr);
    514         return 1;
    515     }
    516 
    517     /* determine size of saved results array, checking for overflows,
    518        allocate and clear the array (set all to zero with calloc()) */
    519     if (syms == 2)              /* iff max == 1 */
    520         num = NULL;             /* won't be saving any results */
    521     else {
    522         size = syms >> 1;
    523         if (size > ((size_t)0 - 1) / (n = (syms - 1) >> 1) ||
    524                 (size *= n, size > ((size_t)0 - 1) / (n = max - 1)) ||
    525                 (size *= n, size > ((size_t)0 - 1) / sizeof(big_t)) ||
    526                 (num = calloc(size, sizeof(big_t))) == NULL) {
    527             fputs("abort: unable to allocate enough memory\n", stderr);
    528             cleanup();
    529             return 1;
    530         }
    531     }
    532 
    533     /* count possible codes for all numbers of symbols, add up counts */
    534     sum = 0;
    535     for (n = 2; n <= syms; n++) {
    536         got = count(n, 1, 2);
    537         sum += got;
    538         if (got == (big_t)0 - 1 || sum < got) {     /* overflow */
    539             fputs("abort: can't count that high!\n", stderr);
    540             cleanup();
    541             return 1;
    542         }
    543         printf("%llu %d-codes\n", got, n);
    544     }
    545     printf("%llu total codes for 2 to %d symbols", sum, syms);
    546     if (max < syms - 1)
    547         printf(" (%d-bit length limit)\n", max);
    548     else
    549         puts(" (no length limit)");
    550 
    551     /* allocate and clear done array for beenhere() */
    552     if (syms == 2)
    553         done = NULL;
    554     else if (size > ((size_t)0 - 1) / sizeof(struct tab) ||
    555              (done = calloc(size, sizeof(struct tab))) == NULL) {
    556         fputs("abort: unable to allocate enough memory\n", stderr);
    557         cleanup();
    558         return 1;
    559     }
    560 
    561     /* find and show maximum inflate table usage */
    562     if (root > max)                 /* reduce root to max length */
    563         root = max;
    564     if ((code_t)syms < ((code_t)1 << (root + 1)))
    565         enough(syms);
    566     else
    567         puts("cannot handle minimum code lengths > root");
    568 
    569     /* done */
    570     cleanup();
    571     return 0;
    572 }
    573