enough.c 23.8 KB
Newer Older
M
Mark Adler 已提交
1 2
/* enough.c -- determine the maximum size of inflate's Huffman code tables over
 * all possible valid and complete Huffman codes, subject to a length limit.
3 4
 * Copyright (C) 2007, 2008, 2012 Mark Adler
 * Version 1.4  18 August 2012  Mark Adler
M
Mark Adler 已提交
5 6 7 8 9 10 11 12 13 14 15 16
 */

/* Version history:
   1.0   3 Jan 2007  First version (derived from codecount.c version 1.4)
   1.1   4 Jan 2007  Use faster incremental table usage computation
                     Prune examine() search on previously visited states
   1.2   5 Jan 2007  Comments clean up
                     As inflate does, decrease root for short codes
                     Refuse cases where inflate would increase root
   1.3  17 Feb 2008  Add argument for initial root table size
                     Fix bug for initial root table size == max - 1
                     Use a macro to compute the history index
17 18 19
   1.4  18 Aug 2012  Avoid shifts more than bits in type (caused endless loop!)
                     Clean up comparisons of different types
                     Clean up code indentation
M
Mark Adler 已提交
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
 */

/*
   Examine all possible Huffman codes for a given number of symbols and a
   maximum code length in bits to determine the maximum table size for zilb's
   inflate.  Only complete Huffman codes are counted.

   Two codes are considered distinct if the vectors of the number of codes per
   length are not identical.  So permutations of the symbol assignments result
   in the same code for the counting, as do permutations of the assignments of
   the bit values to the codes (i.e. only canonical codes are counted).

   We build a code from shorter to longer lengths, determining how many symbols
   are coded at each length.  At each step, we have how many symbols remain to
   be coded, what the last code length used was, and how many bit patterns of
   that length remain unused. Then we add one to the code length and double the
   number of unused patterns to graduate to the next code length.  We then
   assign all portions of the remaining symbols to that code length that
   preserve the properties of a correct and eventually complete code.  Those
   properties are: we cannot use more bit patterns than are available; and when
   all the symbols are used, there are exactly zero possible bit patterns
   remaining.

   The inflate Huffman decoding algorithm uses two-level lookup tables for
   speed.  There is a single first-level table to decode codes up to root bits
   in length (root == 9 in the current inflate implementation).  The table
   has 1 << root entries and is indexed by the next root bits of input.  Codes
   shorter than root bits have replicated table entries, so that the correct
   entry is pointed to regardless of the bits that follow the short code.  If
   the code is longer than root bits, then the table entry points to a second-
   level table.  The size of that table is determined by the longest code with
   that root-bit prefix.  If that longest code has length len, then the table
   has size 1 << (len - root), to index the remaining bits in that set of
   codes.  Each subsequent root-bit prefix then has its own sub-table.  The
   total number of table entries required by the code is calculated
   incrementally as the number of codes at each bit length is populated.  When
   all of the codes are shorter than root bits, then root is reduced to the
   longest code length, resulting in a single, smaller, one-level table.

   The inflate algorithm also provides for small values of root (relative to
   the log2 of the number of symbols), where the shortest code has more bits
   than root.  In that case, root is increased to the length of the shortest
   code.  This program, by design, does not handle that case, so it is verified
   that the number of symbols is less than 2^(root + 1).

   In order to speed up the examination (by about ten orders of magnitude for
   the default arguments), the intermediate states in the build-up of a code
   are remembered and previously visited branches are pruned.  The memory
   required for this will increase rapidly with the total number of symbols and
   the maximum code length in bits.  However this is a very small price to pay
   for the vast speedup.

   First, all of the possible Huffman codes are counted, and reachable
   intermediate states are noted by a non-zero count in a saved-results array.
   Second, the intermediate states that lead to (root + 1) bit or longer codes
   are used to look at all sub-codes from those junctures for their inflate
   memory usage.  (The amount of memory used is not affected by the number of
   codes of root bits or less in length.)  Third, the visited states in the
   construction of those sub-codes and the associated calculation of the table
   size is recalled in order to avoid recalculating from the same juncture.
   Beginning the code examination at (root + 1) bit codes, which is enabled by
   identifying the reachable nodes, accounts for about six of the orders of
   magnitude of improvement for the default arguments.  About another four
   orders of magnitude come from not revisiting previous states.  Out of
   approximately 2x10^16 possible Huffman codes, only about 2x10^6 sub-codes
   need to be examined to cover all of the possible table memory usage cases
   for the default arguments of 286 symbols limited to 15-bit codes.

   Note that an unsigned long long type is used for counting.  It is quite easy
   to exceed the capacity of an eight-byte integer with a large number of
   symbols and a large maximum code length, so multiple-precision arithmetic
   would need to replace the unsigned long long arithmetic in that case.  This
   program will abort if an overflow occurs.  The big_t type identifies where
   the counting takes place.

   An unsigned long long type is also used for calculating the number of
   possible codes remaining at the maximum length.  This limits the maximum
   code length to the number of bits in a long long minus the number of bits
   needed to represent the symbols in a flat code.  The code_t type identifies
   where the bit pattern counting takes place.
 */

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

#define local static

/* special data types */
typedef unsigned long long big_t;   /* type for code counting */
typedef unsigned long long code_t;  /* type for bit pattern counting */
struct tab {                        /* type for been here check */
    size_t len;         /* length of bit vector in char's */
    char *vec;          /* allocated bit vector */
};

/* The array for saving results, num[], is indexed with this triplet:

      syms: number of symbols remaining to code
      left: number of available bit patterns at length len
      len: number of bits in the codes currently being assigned

   Those indices are constrained thusly when saving results:

      syms: 3..totsym (totsym == total symbols to code)
      left: 2..syms - 1, but only the evens (so syms == 8 -> 2, 4, 6)
      len: 1..max - 1 (max == maximum code length in bits)

   syms == 2 is not saved since that immediately leads to a single code.  left
   must be even, since it represents the number of available bit patterns at
   the current length, which is double the number at the previous length.
   left ends at syms-1 since left == syms immediately results in a single code.
   (left > sym is not allowed since that would result in an incomplete code.)
   len is less than max, since the code completes immediately when len == max.

   The offset into the array is calculated for the three indices with the
   first one (syms) being outermost, and the last one (len) being innermost.
   We build the array with length max-1 lists for the len index, with syms-3
   of those for each symbol.  There are totsym-2 of those, with each one
   varying in length as a function of sym.  See the calculation of index in
   count() for the index, and the calculation of size in main() for the size
   of the array.

   For the deflate example of 286 symbols limited to 15-bit codes, the array
   has 284,284 entries, taking up 2.17 MB for an 8-byte big_t.  More than
   half of the space allocated for saved results is actually used -- not all
M
Mark Adler 已提交
147
   possible triplets are reached in the generation of valid Huffman codes.
M
Mark Adler 已提交
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241
 */

/* The array for tracking visited states, done[], is itself indexed identically
   to the num[] array as described above for the (syms, left, len) triplet.
   Each element in the array is further indexed by the (mem, rem) doublet,
   where mem is the amount of inflate table space used so far, and rem is the
   remaining unused entries in the current inflate sub-table.  Each indexed
   element is simply one bit indicating whether the state has been visited or
   not.  Since the ranges for mem and rem are not known a priori, each bit
   vector is of a variable size, and grows as needed to accommodate the visited
   states.  mem and rem are used to calculate a single index in a triangular
   array.  Since the range of mem is expected in the default case to be about
   ten times larger than the range of rem, the array is skewed to reduce the
   memory usage, with eight times the range for mem than for rem.  See the
   calculations for offset and bit in beenhere() for the details.

   For the deflate example of 286 symbols limited to 15-bit codes, the bit
   vectors grow to total approximately 21 MB, in addition to the 4.3 MB done[]
   array itself.
 */

/* Globals to avoid propagating constants or constant pointers recursively */
local int max;          /* maximum allowed bit length for the codes */
local int root;         /* size of base code table in bits */
local int large;        /* largest code table so far */
local size_t size;      /* number of elements in num and done */
local int *code;        /* number of symbols assigned to each bit length */
local big_t *num;       /* saved results array for code counting */
local struct tab *done; /* states already evaluated array */

/* Index function for num[] and done[] */
#define INDEX(i,j,k) (((size_t)((i-1)>>1)*((i-2)>>1)+(j>>1)-1)*(max-1)+k-1)

/* Free allocated space.  Uses globals code, num, and done. */
local void cleanup(void)
{
    size_t n;

    if (done != NULL) {
        for (n = 0; n < size; n++)
            if (done[n].len)
                free(done[n].vec);
        free(done);
    }
    if (num != NULL)
        free(num);
    if (code != NULL)
        free(code);
}

/* Return the number of possible Huffman codes using bit patterns of lengths
   len through max inclusive, coding syms symbols, with left bit patterns of
   length len unused -- return -1 if there is an overflow in the counting.
   Keep a record of previous results in num to prevent repeating the same
   calculation.  Uses the globals max and num. */
local big_t count(int syms, int len, int left)
{
    big_t sum;          /* number of possible codes from this juncture */
    big_t got;          /* value returned from count() */
    int least;          /* least number of syms to use at this juncture */
    int most;           /* most number of syms to use at this juncture */
    int use;            /* number of bit patterns to use in next call */
    size_t index;       /* index of this case in *num */

    /* see if only one possible code */
    if (syms == left)
        return 1;

    /* note and verify the expected state */
    assert(syms > left && left > 0 && len < max);

    /* see if we've done this one already */
    index = INDEX(syms, left, len);
    got = num[index];
    if (got)
        return got;         /* we have -- return the saved result */

    /* we need to use at least this many bit patterns so that the code won't be
       incomplete at the next length (more bit patterns than symbols) */
    least = (left << 1) - syms;
    if (least < 0)
        least = 0;

    /* we can use at most this many bit patterns, lest there not be enough
       available for the remaining symbols at the maximum length (if there were
       no limit to the code length, this would become: most = left - 1) */
    most = (((code_t)left << (max - len)) - syms) /
            (((code_t)1 << (max - len)) - 1);

    /* count all possible codes from this juncture and add them up */
    sum = 0;
    for (use = least; use <= most; use++) {
        got = count(syms - use, len + 1, (left - use) << 1);
        sum += got;
242 243
        if (got == (big_t)0 - 1 || sum < got)   /* overflow */
            return (big_t)0 - 1;
M
Mark Adler 已提交
244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463
    }

    /* verify that all recursive calls are productive */
    assert(sum != 0);

    /* save the result and return it */
    num[index] = sum;
    return sum;
}

/* Return true if we've been here before, set to true if not.  Set a bit in a
   bit vector to indicate visiting this state.  Each (syms,len,left) state
   has a variable size bit vector indexed by (mem,rem).  The bit vector is
   lengthened if needed to allow setting the (mem,rem) bit. */
local int beenhere(int syms, int len, int left, int mem, int rem)
{
    size_t index;       /* index for this state's bit vector */
    size_t offset;      /* offset in this state's bit vector */
    int bit;            /* mask for this state's bit */
    size_t length;      /* length of the bit vector in bytes */
    char *vector;       /* new or enlarged bit vector */

    /* point to vector for (syms,left,len), bit in vector for (mem,rem) */
    index = INDEX(syms, left, len);
    mem -= 1 << root;
    offset = (mem >> 3) + rem;
    offset = ((offset * (offset + 1)) >> 1) + rem;
    bit = 1 << (mem & 7);

    /* see if we've been here */
    length = done[index].len;
    if (offset < length && (done[index].vec[offset] & bit) != 0)
        return 1;       /* done this! */

    /* we haven't been here before -- set the bit to show we have now */

    /* see if we need to lengthen the vector in order to set the bit */
    if (length <= offset) {
        /* if we have one already, enlarge it, zero out the appended space */
        if (length) {
            do {
                length <<= 1;
            } while (length <= offset);
            vector = realloc(done[index].vec, length);
            if (vector != NULL)
                memset(vector + done[index].len, 0, length - done[index].len);
        }

        /* otherwise we need to make a new vector and zero it out */
        else {
            length = 1 << (len - root);
            while (length <= offset)
                length <<= 1;
            vector = calloc(length, sizeof(char));
        }

        /* in either case, bail if we can't get the memory */
        if (vector == NULL) {
            fputs("abort: unable to allocate enough memory\n", stderr);
            cleanup();
            exit(1);
        }

        /* install the new vector */
        done[index].len = length;
        done[index].vec = vector;
    }

    /* set the bit */
    done[index].vec[offset] |= bit;
    return 0;
}

/* Examine all possible codes from the given node (syms, len, left).  Compute
   the amount of memory required to build inflate's decoding tables, where the
   number of code structures used so far is mem, and the number remaining in
   the current sub-table is rem.  Uses the globals max, code, root, large, and
   done. */
local void examine(int syms, int len, int left, int mem, int rem)
{
    int least;          /* least number of syms to use at this juncture */
    int most;           /* most number of syms to use at this juncture */
    int use;            /* number of bit patterns to use in next call */

    /* see if we have a complete code */
    if (syms == left) {
        /* set the last code entry */
        code[len] = left;

        /* complete computation of memory used by this code */
        while (rem < left) {
            left -= rem;
            rem = 1 << (len - root);
            mem += rem;
        }
        assert(rem == left);

        /* if this is a new maximum, show the entries used and the sub-code */
        if (mem > large) {
            large = mem;
            printf("max %d: ", mem);
            for (use = root + 1; use <= max; use++)
                if (code[use])
                    printf("%d[%d] ", code[use], use);
            putchar('\n');
            fflush(stdout);
        }

        /* remove entries as we drop back down in the recursion */
        code[len] = 0;
        return;
    }

    /* prune the tree if we can */
    if (beenhere(syms, len, left, mem, rem))
        return;

    /* we need to use at least this many bit patterns so that the code won't be
       incomplete at the next length (more bit patterns than symbols) */
    least = (left << 1) - syms;
    if (least < 0)
        least = 0;

    /* we can use at most this many bit patterns, lest there not be enough
       available for the remaining symbols at the maximum length (if there were
       no limit to the code length, this would become: most = left - 1) */
    most = (((code_t)left << (max - len)) - syms) /
            (((code_t)1 << (max - len)) - 1);

    /* occupy least table spaces, creating new sub-tables as needed */
    use = least;
    while (rem < use) {
        use -= rem;
        rem = 1 << (len - root);
        mem += rem;
    }
    rem -= use;

    /* examine codes from here, updating table space as we go */
    for (use = least; use <= most; use++) {
        code[len] = use;
        examine(syms - use, len + 1, (left - use) << 1,
                mem + (rem ? 1 << (len - root) : 0), rem << 1);
        if (rem == 0) {
            rem = 1 << (len - root);
            mem += rem;
        }
        rem--;
    }

    /* remove entries as we drop back down in the recursion */
    code[len] = 0;
}

/* Look at all sub-codes starting with root + 1 bits.  Look at only the valid
   intermediate code states (syms, left, len).  For each completed code,
   calculate the amount of memory required by inflate to build the decoding
   tables. Find the maximum amount of memory required and show the code that
   requires that maximum.  Uses the globals max, root, and num. */
local void enough(int syms)
{
    int n;              /* number of remaing symbols for this node */
    int left;           /* number of unused bit patterns at this length */
    size_t index;       /* index of this case in *num */

    /* clear code */
    for (n = 0; n <= max; n++)
        code[n] = 0;

    /* look at all (root + 1) bit and longer codes */
    large = 1 << root;              /* base table */
    if (root < max)                 /* otherwise, there's only a base table */
        for (n = 3; n <= syms; n++)
            for (left = 2; left < n; left += 2)
            {
                /* look at all reachable (root + 1) bit nodes, and the
                   resulting codes (complete at root + 2 or more) */
                index = INDEX(n, left, root + 1);
                if (root + 1 < max && num[index])       /* reachable node */
                    examine(n, root + 1, left, 1 << root, 0);

                /* also look at root bit codes with completions at root + 1
                   bits (not saved in num, since complete), just in case */
                if (num[index - 1] && n <= left << 1)
                    examine((n - left) << 1, root + 1, (n - left) << 1,
                            1 << root, 0);
            }

    /* done */
    printf("done: maximum of %d table entries\n", large);
}

/*
   Examine and show the total number of possible Huffman codes for a given
   maximum number of symbols, initial root table size, and maximum code length
   in bits -- those are the command arguments in that order.  The default
   values are 286, 9, and 15 respectively, for the deflate literal/length code.
   The possible codes are counted for each number of coded symbols from two to
   the maximum.  The counts for each of those and the total number of codes are
   shown.  The maximum number of inflate table entires is then calculated
   across all possible codes.  Each new maximum number of table entries and the
   associated sub-code (starting at root + 1 == 10 bits) is shown.

   To count and examine Huffman codes that are not length-limited, provide a
   maximum length equal to the number of symbols minus one.

   For the deflate literal/length code, use "enough".  For the deflate distance
   code, use "enough 30 6".

   This uses the %llu printf format to print big_t numbers, which assumes that
   big_t is an unsigned long long.  If the big_t type is changed (for example
   to a multiple precision type), the method of printing will also need to be
   updated.
 */
int main(int argc, char **argv)
{
    int syms;           /* total number of symbols to code */
    int n;              /* number of symbols to code for this run */
    big_t got;          /* return value of count() */
    big_t sum;          /* accumulated number of codes over n */
464
    code_t word;        /* for counting bits in code_t */
M
Mark Adler 已提交
465 466 467 468 469 470 471 472

    /* set up globals for cleanup() */
    code = NULL;
    num = NULL;
    done = NULL;

    /* get arguments -- default to the deflate literal/length code */
    syms = 286;
473
    root = 9;
M
Mark Adler 已提交
474 475 476 477 478
    max = 15;
    if (argc > 1) {
        syms = atoi(argv[1]);
        if (argc > 2) {
            root = atoi(argv[2]);
479 480 481
            if (argc > 3)
                max = atoi(argv[3]);
        }
M
Mark Adler 已提交
482 483 484
    }
    if (argc > 4 || syms < 2 || root < 1 || max < 1) {
        fputs("invalid arguments, need: [sym >= 2 [root >= 1 [max >= 1]]]\n",
485
              stderr);
M
Mark Adler 已提交
486 487 488 489 490 491 492 493
        return 1;
    }

    /* if not restricting the code length, the longest is syms - 1 */
    if (max > syms - 1)
        max = syms - 1;

    /* determine the number of bits in a code_t */
494 495
    for (n = 0, word = 1; word; n++, word <<= 1)
        ;
M
Mark Adler 已提交
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))) {
M
Mark Adler 已提交
499 500 501 502 503
        fputs("abort: code length too long for internal types\n", stderr);
        return 1;
    }

    /* reject impossible code requests */
504
    if ((code_t)(syms - 1) > ((code_t)1 << max) - 1) {
M
Mark Adler 已提交
505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537
        fprintf(stderr, "%d symbols cannot be coded in %d bits\n",
                syms, max);
        return 1;
    }

    /* allocate code vector */
    code = calloc(max + 1, sizeof(int));
    if (code == NULL) {
        fputs("abort: unable to allocate enough memory\n", stderr);
        return 1;
    }

    /* determine size of saved results array, checking for overflows,
       allocate and clear the array (set all to zero with calloc()) */
    if (syms == 2)              /* iff max == 1 */
        num = NULL;             /* won't be saving any results */
    else {
        size = syms >> 1;
        if (size > ((size_t)0 - 1) / (n = (syms - 1) >> 1) ||
                (size *= n, size > ((size_t)0 - 1) / (n = max - 1)) ||
                (size *= n, size > ((size_t)0 - 1) / sizeof(big_t)) ||
                (num = calloc(size, sizeof(big_t))) == NULL) {
            fputs("abort: unable to allocate enough memory\n", stderr);
            cleanup();
            return 1;
        }
    }

    /* count possible codes for all numbers of symbols, add up counts */
    sum = 0;
    for (n = 2; n <= syms; n++) {
        got = count(n, 1, 2);
        sum += got;
538
        if (got == (big_t)0 - 1 || sum < got) {     /* overflow */
M
Mark Adler 已提交
539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561
            fputs("abort: can't count that high!\n", stderr);
            cleanup();
            return 1;
        }
        printf("%llu %d-codes\n", got, n);
    }
    printf("%llu total codes for 2 to %d symbols", sum, syms);
    if (max < syms - 1)
        printf(" (%d-bit length limit)\n", max);
    else
        puts(" (no length limit)");

    /* allocate and clear done array for beenhere() */
    if (syms == 2)
        done = NULL;
    else if (size > ((size_t)0 - 1) / sizeof(struct tab) ||
             (done = calloc(size, sizeof(struct tab))) == NULL) {
        fputs("abort: unable to allocate enough memory\n", stderr);
        cleanup();
        return 1;
    }

    /* find and show maximum inflate table usage */
562 563 564
    if (root > max)                 /* reduce root to max length */
        root = max;
    if ((code_t)syms < ((code_t)1 << (root + 1)))
M
Mark Adler 已提交
565 566 567 568 569 570 571 572
        enough(syms);
    else
        puts("cannot handle minimum code lengths > root");

    /* done */
    cleanup();
    return 0;
}