! 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. */
| 17109 | alpha <= e_c + e + q <= gamma. |
| 17110 | */ |
| 17111 | inline cached_power get_cached_power_for_binary_exponent(int e) |
| 17112 | { |
| 17113 | // Now |
| 17114 | // |
| 17115 | // alpha <= e_c + e + q <= gamma (1) |
| 17116 | // ==> f_c * 2^alpha <= c * 2^e * 2^q |
| 17117 | // |
| 17118 | // and since the c's are normalized, 2^(q-1) <= f_c, |
| 17119 | // |
| 17120 | // ==> 2^(q - 1 + alpha) <= c * 2^(e + q) |
| 17121 | // ==> 2^(alpha - e - 1) <= c |
| 17122 | // |
| 17123 | // If c were an exact power of ten, i.e. c = 10^k, one may determine k as |
| 17124 | // |
| 17125 | // k = ceil( log_10( 2^(alpha - e - 1) ) ) |
| 17126 | // = ceil( (alpha - e - 1) * log_10(2) ) |
| 17127 | // |
| 17128 | // From the paper: |
| 17129 | // "In theory the result of the procedure could be wrong since c is rounded, |
| 17130 | // and the computation itself is approximated [...]. In practice, however, |
| 17131 | // this simple function is sufficient." |
| 17132 | // |
| 17133 | // For IEEE double precision floating-point numbers converted into |
| 17134 | // normalized diyfp's w = f * 2^e, with q = 64, |
| 17135 | // |
| 17136 | // e >= -1022 (min IEEE exponent) |
| 17137 | // -52 (p - 1) |
| 17138 | // -52 (p - 1, possibly normalize denormal IEEE numbers) |
| 17139 | // -11 (normalize the diyfp) |
| 17140 | // = -1137 |
| 17141 | // |
| 17142 | // and |
| 17143 | // |
| 17144 | // e <= +1023 (max IEEE exponent) |
| 17145 | // -52 (p - 1) |
| 17146 | // -11 (normalize the diyfp) |
| 17147 | // = 960 |
| 17148 | // |
| 17149 | // This binary exponent range [-1137,960] results in a decimal exponent |
| 17150 | // range [-307,324]. One does not need to store a cached power for each |
| 17151 | // k in this range. For each such k it suffices to find a cached power |
| 17152 | // such that the exponent of the product lies in [alpha,gamma]. |
| 17153 | // This implies that the difference of the decimal exponents of adjacent |
| 17154 | // table entries must be less than or equal to |
| 17155 | // |
| 17156 | // floor( (gamma - alpha) * log_10(2) ) = 8. |
| 17157 | // |
| 17158 | // (A smaller distance gamma-alpha would require a larger table.) |
| 17159 | |
| 17160 | // NB: |
| 17161 | // Actually this function returns c, such that -60 <= e_c + e + 64 <= -34. |
| 17162 | |
| 17163 | constexpr int kCachedPowersMinDecExp = -300; |
| 17164 | constexpr int kCachedPowersDecStep = 8; |
| 17165 | |
| 17166 | static constexpr std::array<cached_power, 79> kCachedPowers = |
| 17167 | { |
| 17168 | { |