examples/enough.c

   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 struct {
 171     int max;            /* maximum allowed bit length for the codes */
 172     int root;           /* size of base code table in bits */
 173     int large;          /* largest code table so far */
 174     size_t size;        /* number of elements in num and done */
 175     int *code;          /* number of symbols assigned to each bit length */
 176     big_t *num;         /* saved results array for code counting */
 177     struct tab *done;   /* states already evaluated array */
 178 } g;
 179
 180 /* Index function for num[] and done[] */
 181 #define INDEX(i,j,k) (((size_t)((i-1)>>1)*((i-2)>>1)+(j>>1)-1)*(g.max-1)+k-1)
 182
 183 /* Free allocated space.  Uses globals code, num, and done. */
 184 local void cleanup(void)
 185 {
 186     size_t n;
 187
 188     if (g.done != NULL) {
 189         for (n = 0; n < g.size; n++)
 190             if (g.done[n].len)
 191                 free(g.done[n].vec);
 192         free(g.done);
 193     }
 194     if (g.num != NULL)
 195         free(g.num);
 196     if (g.code != NULL)
 197         free(g.code);
 198 }
 199
 200 /* Return the number of possible Huffman codes using bit patterns of lengths
 201    len through max inclusive, coding syms symbols, with left bit patterns of
 202    length len unused -- return -1 if there is an overflow in the counting.
 203    Keep a record of previous results in num to prevent repeating the same
 204    calculation.  Uses the globals max and num. */
 205 local big_t count(int syms, int len, int left)
 206 {
 207     big_t sum;          /* number of possible codes from this juncture */
 208     big_t got;          /* value returned from count() */
 209     int least;          /* least number of syms to use at this juncture */
 210     int most;           /* most number of syms to use at this juncture */
 211     int use;            /* number of bit patterns to use in next call */
 212     size_t index;       /* index of this case in *num */
 213
 214     /* see if only one possible code */
 215     if (syms == left)
 216         return 1;
 217
 218     /* note and verify the expected state */
 219     assert(syms > left && left > 0 && len < g.max);
 220
 221     /* see if we've done this one already */
 222     index = INDEX(syms, left, len);
 223     got = g.num[index];
 224     if (got)
 225         return got;         /* we have -- return the saved result */
 226
 227     /* we need to use at least this many bit patterns so that the code won't be
 228        incomplete at the next length (more bit patterns than symbols) */
 229     least = (left << 1) - syms;
 230     if (least < 0)
 231         least = 0;
 232
 233     /* we can use at most this many bit patterns, lest there not be enough
 234        available for the remaining symbols at the maximum length (if there were
 235        no limit to the code length, this would become: most = left - 1) */
 236     most = (((code_t)left << (g.max - len)) - syms) /
 237             (((code_t)1 << (g.max - len)) - 1);
 238
 239     /* count all possible codes from this juncture and add them up */
 240     sum = 0;
 241     for (use = least; use <= most; use++) {
 242         got = count(syms - use, len + 1, (left - use) << 1);
 243         sum += got;
 244         if (got == (big_t)0 - 1 || sum < got)   /* overflow */
 245             return (big_t)0 - 1;
 246     }
 247
 248     /* verify that all recursive calls are productive */
 249     assert(sum != 0);
 250
 251     /* save the result and return it */
 252     g.num[index] = sum;
 253     return sum;
 254 }
 255
 256 /* Return true if we've been here before, set to true if not.  Set a bit in a
 257    bit vector to indicate visiting this state.  Each (syms,len,left) state
 258    has a variable size bit vector indexed by (mem,rem).  The bit vector is
 259    lengthened if needed to allow setting the (mem,rem) bit. */
 260 local int beenhere(int syms, int len, int left, int mem, int rem)
 261 {
 262     size_t index;       /* index for this state's bit vector */
 263     size_t offset;      /* offset in this state's bit vector */
 264     int bit;            /* mask for this state's bit */
 265     size_t length;      /* length of the bit vector in bytes */
 266     char *vector;       /* new or enlarged bit vector */
 267
 268     /* point to vector for (syms,left,len), bit in vector for (mem,rem) */
 269     index = INDEX(syms, left, len);
 270     mem -= 1 << g.root;
 271     offset = (mem >> 3) + rem;
 272     offset = ((offset * (offset + 1)) >> 1) + rem;
 273     bit = 1 << (mem & 7);
 274
 275     /* see if we've been here */
 276     length = g.done[index].len;
 277     if (offset < length && (g.done[index].vec[offset] & bit) != 0)
 278         return 1;       /* done this! */
 279
 280     /* we haven't been here before -- set the bit to show we have now */
 281
 282     /* see if we need to lengthen the vector in order to set the bit */
 283     if (length <= offset) {
 284         /* if we have one already, enlarge it, zero out the appended space */
 285         if (length) {
 286             do {
 287                 length <<= 1;
 288             } while (length <= offset);
 289             vector = realloc(g.done[index].vec, length);
 290             if (vector != NULL)
 291                 memset(vector + g.done[index].len, 0,
 292                        length - g.done[index].len);
 293         }
 294
 295         /* otherwise we need to make a new vector and zero it out */
 296         else {
 297             length = 1 << (len - g.root);
 298             while (length <= offset)
 299                 length <<= 1;
 300             vector = calloc(length, sizeof(char));
 301         }
 302
 303         /* in either case, bail if we can't get the memory */
 304         if (vector == NULL) {
 305             fputs("abort: unable to allocate enough memory\n", stderr);
 306             cleanup();
 307             exit(1);
 308         }
 309
 310         /* install the new vector */
 311         g.done[index].len = length;
 312         g.done[index].vec = vector;
 313     }
 314
 315     /* set the bit */
 316     g.done[index].vec[offset] |= bit;
 317     return 0;
 318 }
 319
 320 /* Examine all possible codes from the given node (syms, len, left).  Compute
 321    the amount of memory required to build inflate's decoding tables, where the
 322    number of code structures used so far is mem, and the number remaining in
 323    the current sub-table is rem.  Uses the globals max, code, root, large, and
 324    done. */
 325 local void examine(int syms, int len, int left, int mem, int rem)
 326 {
 327     int least;          /* least number of syms to use at this juncture */
 328     int most;           /* most number of syms to use at this juncture */
 329     int use;            /* number of bit patterns to use in next call */
 330
 331     /* see if we have a complete code */
 332     if (syms == left) {
 333         /* set the last code entry */
 334         g.code[len] = left;
 335
 336         /* complete computation of memory used by this code */
 337         while (rem < left) {
 338             left -= rem;
 339             rem = 1 << (len - g.root);
 340             mem += rem;
 341         }
 342         assert(rem == left);
 343
 344         /* if this is a new maximum, show the entries used and the sub-code */
 345         if (mem > g.large) {
 346             g.large = mem;
 347             printf("max %d: ", mem);
 348             for (use = g.root + 1; use <= g.max; use++)
 349                 if (g.code[use])
 350                     printf("%d[%d] ", g.code[use], use);
 351             putchar('\n');
 352             fflush(stdout);
 353         }
 354
 355         /* remove entries as we drop back down in the recursion */
 356         g.code[len] = 0;
 357         return;
 358     }
 359
 360     /* prune the tree if we can */
 361     if (beenhere(syms, len, left, mem, rem))
 362         return;
 363
 364     /* we need to use at least this many bit patterns so that the code won't be
 365        incomplete at the next length (more bit patterns than symbols) */
 366     least = (left << 1) - syms;
 367     if (least < 0)
 368         least = 0;
 369
 370     /* we can use at most this many bit patterns, lest there not be enough
 371        available for the remaining symbols at the maximum length (if there were
 372        no limit to the code length, this would become: most = left - 1) */
 373     most = (((code_t)left << (g.max - len)) - syms) /
 374             (((code_t)1 << (g.max - len)) - 1);
 375
 376     /* occupy least table spaces, creating new sub-tables as needed */
 377     use = least;
 378     while (rem < use) {
 379         use -= rem;
 380         rem = 1 << (len - g.root);
 381         mem += rem;
 382     }
 383     rem -= use;
 384
 385     /* examine codes from here, updating table space as we go */
 386     for (use = least; use <= most; use++) {
 387         g.code[len] = use;
 388         examine(syms - use, len + 1, (left - use) << 1,
 389                 mem + (rem ? 1 << (len - g.root) : 0), rem << 1);
 390         if (rem == 0) {
 391             rem = 1 << (len - g.root);
 392             mem += rem;
 393         }
 394         rem--;
 395     }
 396
 397     /* remove entries as we drop back down in the recursion */
 398     g.code[len] = 0;
 399 }
 400
 401 /* Look at all sub-codes starting with root + 1 bits.  Look at only the valid
 402    intermediate code states (syms, left, len).  For each completed code,
 403    calculate the amount of memory required by inflate to build the decoding
 404    tables. Find the maximum amount of memory required and show the code that
 405    requires that maximum.  Uses the globals max, root, and num. */
 406 local void enough(int syms)
 407 {
 408     int n;              /* number of remaing symbols for this node */
 409     int left;           /* number of unused bit patterns at this length */
 410     size_t index;       /* index of this case in *num */
 411
 412     /* clear code */
 413     for (n = 0; n <= g.max; n++)
 414         g.code[n] = 0;
 415
 416     /* look at all (root + 1) bit and longer codes */
 417     g.large = 1 << g.root;          /* base table */
 418     if (g.root < g.max)             /* otherwise, there's only a base table */
 419         for (n = 3; n <= syms; n++)
 420             for (left = 2; left < n; left += 2)
 421             {
 422                 /* look at all reachable (root + 1) bit nodes, and the
 423                    resulting codes (complete at root + 2 or more) */
 424                 index = INDEX(n, left, g.root + 1);
 425                 if (g.root + 1 < g.max && g.num[index]) /* reachable node */
 426                     examine(n, g.root + 1, left, 1 << g.root, 0);
 427
 428                 /* also look at root bit codes with completions at root + 1
 429                    bits (not saved in num, since complete), just in case */
 430                 if (g.num[index - 1] && n <= left << 1)
 431                     examine((n - left) << 1, g.root + 1, (n - left) << 1,
 432                             1 << g.root, 0);
 433             }
 434
 435     /* done */
 436     printf("done: maximum of %d table entries\n", g.large);
 437 }
 438
 439 /*
 440    Examine and show the total number of possible Huffman codes for a given
 441    maximum number of symbols, initial root table size, and maximum code length
 442    in bits -- those are the command arguments in that order.  The default
 443    values are 286, 9, and 15 respectively, for the deflate literal/length code.
 444    The possible codes are counted for each number of coded symbols from two to
 445    the maximum.  The counts for each of those and the total number of codes are
 446    shown.  The maximum number of inflate table entires is then calculated
 447    across all possible codes.  Each new maximum number of table entries and the
 448    associated sub-code (starting at root + 1 == 10 bits) is shown.
 449
 450    To count and examine Huffman codes that are not length-limited, provide a
 451    maximum length equal to the number of symbols minus one.
 452
 453    For the deflate literal/length code, use "enough".  For the deflate distance
 454    code, use "enough 30 6".
 455
 456    This uses the %llu printf format to print big_t numbers, which assumes that
 457    big_t is an unsigned long long.  If the big_t type is changed (for example
 458    to a multiple precision type), the method of printing will also need to be
 459    updated.
 460  */
 461 int main(int argc, char **argv)
 462 {
 463     int syms;           /* total number of symbols to code */
 464     int n;              /* number of symbols to code for this run */
 465     big_t got;          /* return value of count() */
 466     big_t sum;          /* accumulated number of codes over n */
 467     code_t word;        /* for counting bits in code_t */
 468
 469     /* set up globals for cleanup() */
 470     g.code = NULL;
 471     g.num = NULL;
 472     g.done = NULL;
 473
 474     /* get arguments -- default to the deflate literal/length code */
 475     syms = 286;
 476     g.root = 9;
 477     g.max = 15;
 478     if (argc > 1) {
 479         syms = atoi(argv[1]);
 480         if (argc > 2) {
 481             g.root = atoi(argv[2]);
 482             if (argc > 3)
 483                 g.max = atoi(argv[3]);
 484         }
 485     }
 486     if (argc > 4 || syms < 2 || g.root < 1 || g.max < 1) {
 487         fputs("invalid arguments, need: [sym >= 2 [root >= 1 [max >= 1]]]\n",
 488               stderr);
 489         return 1;
 490     }
 491
 492     /* if not restricting the code length, the longest is syms - 1 */
 493     if (g.max > syms - 1)
 494         g.max = syms - 1;
 495
 496     /* determine the number of bits in a code_t */
 497     for (n = 0, word = 1; word; n++, word <<= 1)
 498         ;
 499
 500     /* make sure that the calculation of most will not overflow */
 501     if (g.max > n || (code_t)(syms - 2) >= (((code_t)0 - 1) >> (g.max - 1))) {
 502         fputs("abort: code length too long for internal types\n", stderr);
 503         return 1;
 504     }
 505
 506     /* reject impossible code requests */
 507     if ((code_t)(syms - 1) > ((code_t)1 << g.max) - 1) {
 508         fprintf(stderr, "%d symbols cannot be coded in %d bits\n",
 509                 syms, g.max);
 510         return 1;
 511     }
 512
 513     /* allocate code vector */
 514     g.code = calloc(g.max + 1, sizeof(int));
 515     if (g.code == NULL) {
 516         fputs("abort: unable to allocate enough memory\n", stderr);
 517         return 1;
 518     }
 519
 520     /* determine size of saved results array, checking for overflows,
 521        allocate and clear the array (set all to zero with calloc()) */
 522     if (syms == 2)              /* iff max == 1 */
 523         g.num = NULL;           /* won't be saving any results */
 524     else {
 525         g.size = syms >> 1;
 526         if (g.size > ((size_t)0 - 1) / (n = (syms - 1) >> 1) ||
 527                 (g.size *= n, g.size > ((size_t)0 - 1) / (n = g.max - 1)) ||
 528                 (g.size *= n, g.size > ((size_t)0 - 1) / sizeof(big_t)) ||
 529                 (g.num = calloc(g.size, sizeof(big_t))) == NULL) {
 530             fputs("abort: unable to allocate enough memory\n", stderr);
 531             cleanup();
 532             return 1;
 533         }
 534     }
 535
 536     /* count possible codes for all numbers of symbols, add up counts */
 537     sum = 0;
 538     for (n = 2; n <= syms; n++) {
 539         got = count(n, 1, 2);
 540         sum += got;
 541         if (got == (big_t)0 - 1 || sum < got) {     /* overflow */
 542             fputs("abort: can't count that high!\n", stderr);
 543             cleanup();
 544             return 1;
 545         }
 546         printf("%llu %d-codes\n", got, n);
 547     }
 548     printf("%llu total codes for 2 to %d symbols", sum, syms);
 549     if (g.max < syms - 1)
 550         printf(" (%d-bit length limit)\n", g.max);
 551     else
 552         puts(" (no length limit)");
 553
 554     /* allocate and clear done array for beenhere() */
 555     if (syms == 2)
 556         g.done = NULL;
 557     else if (g.size > ((size_t)0 - 1) / sizeof(struct tab) ||
 558              (g.done = calloc(g.size, sizeof(struct tab))) == NULL) {
 559         fputs("abort: unable to allocate enough memory\n", stderr);
 560         cleanup();
 561         return 1;
 562     }
 563
 564     /* find and show maximum inflate table usage */
 565     if (g.root > g.max)             /* reduce root to max length */
 566         g.root = g.max;
 567     if ((code_t)syms < ((code_t)1 << (g.root + 1)))
 568         enough(syms);
 569     else
 570         puts("cannot handle minimum code lengths > root");
 571
 572     /* done */
 573     cleanup();
 574     return 0;
 575 }