! For a normalized diyfp w = f * 2^e, this function returns a (normalized) cached power-of-ten c = f_c * 2^e_c, such that the exponent of the product w * c satisfies (Definition 3.2 from [1]) alpha <= e_c + e + q <= gamma. */
| 14657 | alpha <= e_c + e + q <= gamma. |
| 14658 | */ |
| 14659 | inline cached_power get_cached_power_for_binary_exponent(int e) |
| 14660 | { |
| 14661 | // Now |
| 14662 | // |
| 14663 | // alpha <= e_c + e + q <= gamma (1) |
| 14664 | // ==> f_c * 2^alpha <= c * 2^e * 2^q |
| 14665 | // |
| 14666 | // and since the c's are normalized, 2^(q-1) <= f_c, |
| 14667 | // |
| 14668 | // ==> 2^(q - 1 + alpha) <= c * 2^(e + q) |
| 14669 | // ==> 2^(alpha - e - 1) <= c |
| 14670 | // |
| 14671 | // If c were an exact power of ten, i.e. c = 10^k, one may determine k as |
| 14672 | // |
| 14673 | // k = ceil( log_10( 2^(alpha - e - 1) ) ) |
| 14674 | // = ceil( (alpha - e - 1) * log_10(2) ) |
| 14675 | // |
| 14676 | // From the paper: |
| 14677 | // "In theory the result of the procedure could be wrong since c is rounded, |
| 14678 | // and the computation itself is approximated [...]. In practice, however, |
| 14679 | // this simple function is sufficient." |
| 14680 | // |
| 14681 | // For IEEE double precision floating-point numbers converted into |
| 14682 | // normalized diyfp's w = f * 2^e, with q = 64, |
| 14683 | // |
| 14684 | // e >= -1022 (min IEEE exponent) |
| 14685 | // -52 (p - 1) |
| 14686 | // -52 (p - 1, possibly normalize denormal IEEE numbers) |
| 14687 | // -11 (normalize the diyfp) |
| 14688 | // = -1137 |
| 14689 | // |
| 14690 | // and |
| 14691 | // |
| 14692 | // e <= +1023 (max IEEE exponent) |
| 14693 | // -52 (p - 1) |
| 14694 | // -11 (normalize the diyfp) |
| 14695 | // = 960 |
| 14696 | // |
| 14697 | // This binary exponent range [-1137,960] results in a decimal exponent |
| 14698 | // range [-307,324]. One does not need to store a cached power for each |
| 14699 | // k in this range. For each such k it suffices to find a cached power |
| 14700 | // such that the exponent of the product lies in [alpha,gamma]. |
| 14701 | // This implies that the difference of the decimal exponents of adjacent |
| 14702 | // table entries must be less than or equal to |
| 14703 | // |
| 14704 | // floor( (gamma - alpha) * log_10(2) ) = 8. |
| 14705 | // |
| 14706 | // (A smaller distance gamma-alpha would require a larger table.) |
| 14707 | |
| 14708 | // NB: |
| 14709 | // Actually this function returns c, such that -60 <= e_c + e + 64 <= -34. |
| 14710 | |
| 14711 | constexpr int kCachedPowersMinDecExp = -300; |
| 14712 | constexpr int kCachedPowersDecStep = 8; |
| 14713 | |
| 14714 | static constexpr std::array<cached_power, 79> kCachedPowers = |
| 14715 | { |
| 14716 | { |