! 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. */
| 14732 | alpha <= e_c + e + q <= gamma. |
| 14733 | */ |
| 14734 | inline cached_power get_cached_power_for_binary_exponent(int e) |
| 14735 | { |
| 14736 | // Now |
| 14737 | // |
| 14738 | // alpha <= e_c + e + q <= gamma (1) |
| 14739 | // ==> f_c * 2^alpha <= c * 2^e * 2^q |
| 14740 | // |
| 14741 | // and since the c's are normalized, 2^(q-1) <= f_c, |
| 14742 | // |
| 14743 | // ==> 2^(q - 1 + alpha) <= c * 2^(e + q) |
| 14744 | // ==> 2^(alpha - e - 1) <= c |
| 14745 | // |
| 14746 | // If c were an exact power of ten, i.e. c = 10^k, one may determine k as |
| 14747 | // |
| 14748 | // k = ceil( log_10( 2^(alpha - e - 1) ) ) |
| 14749 | // = ceil( (alpha - e - 1) * log_10(2) ) |
| 14750 | // |
| 14751 | // From the paper: |
| 14752 | // "In theory the result of the procedure could be wrong since c is rounded, |
| 14753 | // and the computation itself is approximated [...]. In practice, however, |
| 14754 | // this simple function is sufficient." |
| 14755 | // |
| 14756 | // For IEEE double precision floating-point numbers converted into |
| 14757 | // normalized diyfp's w = f * 2^e, with q = 64, |
| 14758 | // |
| 14759 | // e >= -1022 (min IEEE exponent) |
| 14760 | // -52 (p - 1) |
| 14761 | // -52 (p - 1, possibly normalize denormal IEEE numbers) |
| 14762 | // -11 (normalize the diyfp) |
| 14763 | // = -1137 |
| 14764 | // |
| 14765 | // and |
| 14766 | // |
| 14767 | // e <= +1023 (max IEEE exponent) |
| 14768 | // -52 (p - 1) |
| 14769 | // -11 (normalize the diyfp) |
| 14770 | // = 960 |
| 14771 | // |
| 14772 | // This binary exponent range [-1137,960] results in a decimal exponent |
| 14773 | // range [-307,324]. One does not need to store a cached power for each |
| 14774 | // k in this range. For each such k it suffices to find a cached power |
| 14775 | // such that the exponent of the product lies in [alpha,gamma]. |
| 14776 | // This implies that the difference of the decimal exponents of adjacent |
| 14777 | // table entries must be less than or equal to |
| 14778 | // |
| 14779 | // floor( (gamma - alpha) * log_10(2) ) = 8. |
| 14780 | // |
| 14781 | // (A smaller distance gamma-alpha would require a larger table.) |
| 14782 | |
| 14783 | // NB: |
| 14784 | // Actually this function returns c, such that -60 <= e_c + e + 64 <= -34. |
| 14785 | |
| 14786 | constexpr int kCachedPowersMinDecExp = -300; |
| 14787 | constexpr int kCachedPowersDecStep = 8; |
| 14788 | |
| 14789 | static constexpr std::array<cached_power, 79> kCachedPowers = |
| 14790 | { |
| 14791 | { |