Compute floor(log2(x!)), exactly up to x=57; an underestimate up to x=2^32-1. */
| 17 | |
| 18 | /** Compute floor(log2(x!)), exactly up to x=57; an underestimate up to x=2^32-1. */ |
| 19 | uint64_t Log2Factorial(uint32_t x) { |
| 20 | //! Values of floor(106*log2(1 + i/32)) for i=0..31 |
| 21 | static constexpr uint8_t T[32] = { |
| 22 | 0, 4, 9, 13, 18, 22, 26, 30, 34, 37, 41, 45, 48, 52, 55, 58, 62, 65, 68, |
| 23 | 71, 74, 77, 80, 82, 85, 88, 90, 93, 96, 98, 101, 103 |
| 24 | }; |
| 25 | int bits = CountBits(x, 32); |
| 26 | // Compute an (under)estimate of floor(106*log2(x)). |
| 27 | // This works by relying on floor(log2(x)) = countbits(x)-1, and adding |
| 28 | // precision using the top 6 bits of x (the highest one of which is always |
| 29 | // one). |
| 30 | unsigned l2_106 = 106 * (bits - 1) + T[((x << (32 - bits)) >> 26) & 31]; |
| 31 | // Based on Stirling approximation for log2(x!): |
| 32 | // log2(x!) = log(x!) / log(2) |
| 33 | // = ((x + 1/2) * log(x) - x + log(2*pi)/2 + ...) / log(2) |
| 34 | // = (x + 1/2) * log2(x) - x/log(2) + log2(2*pi)/2 + ... |
| 35 | // = 1/2*(2*x+1)*log2(x) - (1/log(2))*x + log2(2*pi)/2 + ... |
| 36 | // = 1/212*(2*x+1)*(106*log2(x)) + (-1/log(2))*x + log2(2*pi)/2 + ... |
| 37 | // where 418079/88632748 is exactly 1/212 |
| 38 | // -127870026/88632748 is slightly less than -1/log(2) |
| 39 | // 117504694/88632748 is less than log2(2*pi)/2 |
| 40 | // A correction term is only needed for x < 3. |
| 41 | // |
| 42 | // See doc/log2_factorial.sage for how these constants were obtained. |
| 43 | return (418079 * (2 * uint64_t{x} + 1) * l2_106 - 127870026 * uint64_t{x} + 117504694 + 88632748 * (x < 3)) / 88632748; |
| 44 | } |
| 45 | |
| 46 | /** Compute floor(log2(2^(bits * capacity) / sum((2^bits - 1) choose k, k=0..capacity))), for bits>1 |
| 47 | * |
no test coverage detected