| 1 | // Copyright 2012 the V8 project authors. All rights reserved. | 
|---|
| 2 | // Redistribution and use in source and binary forms, with or without | 
|---|
| 3 | // modification, are permitted provided that the following conditions are | 
|---|
| 4 | // met: | 
|---|
| 5 | // | 
|---|
| 6 | //     * Redistributions of source code must retain the above copyright | 
|---|
| 7 | //       notice, this list of conditions and the following disclaimer. | 
|---|
| 8 | //     * Redistributions in binary form must reproduce the above | 
|---|
| 9 | //       copyright notice, this list of conditions and the following | 
|---|
| 10 | //       disclaimer in the documentation and/or other materials provided | 
|---|
| 11 | //       with the distribution. | 
|---|
| 12 | //     * Neither the name of Google Inc. nor the names of its | 
|---|
| 13 | //       contributors may be used to endorse or promote products derived | 
|---|
| 14 | //       from this software without specific prior written permission. | 
|---|
| 15 | // | 
|---|
| 16 | // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS | 
|---|
| 17 | // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT | 
|---|
| 18 | // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR | 
|---|
| 19 | // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT | 
|---|
| 20 | // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, | 
|---|
| 21 | // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT | 
|---|
| 22 | // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, | 
|---|
| 23 | // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY | 
|---|
| 24 | // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT | 
|---|
| 25 | // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE | 
|---|
| 26 | // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | 
|---|
| 27 |  | 
|---|
| 28 | #include "fast-dtoa.h" | 
|---|
| 29 |  | 
|---|
| 30 | #include "cached-powers.h" | 
|---|
| 31 | #include "diy-fp.h" | 
|---|
| 32 | #include "ieee.h" | 
|---|
| 33 |  | 
|---|
| 34 | namespace double_conversion { | 
|---|
| 35 |  | 
|---|
| 36 | // The minimal and maximal target exponent define the range of w's binary | 
|---|
| 37 | // exponent, where 'w' is the result of multiplying the input by a cached power | 
|---|
| 38 | // of ten. | 
|---|
| 39 | // | 
|---|
| 40 | // A different range might be chosen on a different platform, to optimize digit | 
|---|
| 41 | // generation, but a smaller range requires more powers of ten to be cached. | 
|---|
| 42 | static const int kMinimalTargetExponent = -60; | 
|---|
| 43 | static const int kMaximalTargetExponent = -32; | 
|---|
| 44 |  | 
|---|
| 45 |  | 
|---|
| 46 | // Adjusts the last digit of the generated number, and screens out generated | 
|---|
| 47 | // solutions that may be inaccurate. A solution may be inaccurate if it is | 
|---|
| 48 | // outside the safe interval, or if we cannot prove that it is closer to the | 
|---|
| 49 | // input than a neighboring representation of the same length. | 
|---|
| 50 | // | 
|---|
| 51 | // Input: * buffer containing the digits of too_high / 10^kappa | 
|---|
| 52 | //        * the buffer's length | 
|---|
| 53 | //        * distance_too_high_w == (too_high - w).f() * unit | 
|---|
| 54 | //        * unsafe_interval == (too_high - too_low).f() * unit | 
|---|
| 55 | //        * rest = (too_high - buffer * 10^kappa).f() * unit | 
|---|
| 56 | //        * ten_kappa = 10^kappa * unit | 
|---|
| 57 | //        * unit = the common multiplier | 
|---|
| 58 | // Output: returns true if the buffer is guaranteed to contain the closest | 
|---|
| 59 | //    representable number to the input. | 
|---|
| 60 | //  Modifies the generated digits in the buffer to approach (round towards) w. | 
|---|
| 61 | static bool RoundWeed(Vector<char> buffer, | 
|---|
| 62 | int length, | 
|---|
| 63 | uint64_t distance_too_high_w, | 
|---|
| 64 | uint64_t unsafe_interval, | 
|---|
| 65 | uint64_t rest, | 
|---|
| 66 | uint64_t ten_kappa, | 
|---|
| 67 | uint64_t unit) { | 
|---|
| 68 | uint64_t small_distance = distance_too_high_w - unit; | 
|---|
| 69 | uint64_t big_distance = distance_too_high_w + unit; | 
|---|
| 70 | // Let w_low  = too_high - big_distance, and | 
|---|
| 71 | //     w_high = too_high - small_distance. | 
|---|
| 72 | // Note: w_low < w < w_high | 
|---|
| 73 | // | 
|---|
| 74 | // The real w (* unit) must lie somewhere inside the interval | 
|---|
| 75 | // ]w_low; w_high[ (often written as "(w_low; w_high)") | 
|---|
| 76 |  | 
|---|
| 77 | // Basically the buffer currently contains a number in the unsafe interval | 
|---|
| 78 | // ]too_low; too_high[ with too_low < w < too_high | 
|---|
| 79 | // | 
|---|
| 80 | //  too_high - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | 
|---|
| 81 | //                     ^v 1 unit            ^      ^                 ^      ^ | 
|---|
| 82 | //  boundary_high ---------------------     .      .                 .      . | 
|---|
| 83 | //                     ^v 1 unit            .      .                 .      . | 
|---|
| 84 | //   - - - - - - - - - - - - - - - - - - -  +  - - + - - - - - -     .      . | 
|---|
| 85 | //                                          .      .         ^       .      . | 
|---|
| 86 | //                                          .  big_distance  .       .      . | 
|---|
| 87 | //                                          .      .         .       .    rest | 
|---|
| 88 | //                              small_distance     .         .       .      . | 
|---|
| 89 | //                                          v      .         .       .      . | 
|---|
| 90 | //  w_high - - - - - - - - - - - - - - - - - -     .         .       .      . | 
|---|
| 91 | //                     ^v 1 unit                   .         .       .      . | 
|---|
| 92 | //  w ----------------------------------------     .         .       .      . | 
|---|
| 93 | //                     ^v 1 unit                   v         .       .      . | 
|---|
| 94 | //  w_low  - - - - - - - - - - - - - - - - - - - - -         .       .      . | 
|---|
| 95 | //                                                           .       .      v | 
|---|
| 96 | //  buffer --------------------------------------------------+-------+-------- | 
|---|
| 97 | //                                                           .       . | 
|---|
| 98 | //                                                  safe_interval    . | 
|---|
| 99 | //                                                           v       . | 
|---|
| 100 | //   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     . | 
|---|
| 101 | //                     ^v 1 unit                                     . | 
|---|
| 102 | //  boundary_low -------------------------                     unsafe_interval | 
|---|
| 103 | //                     ^v 1 unit                                     v | 
|---|
| 104 | //  too_low  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | 
|---|
| 105 | // | 
|---|
| 106 | // | 
|---|
| 107 | // Note that the value of buffer could lie anywhere inside the range too_low | 
|---|
| 108 | // to too_high. | 
|---|
| 109 | // | 
|---|
| 110 | // boundary_low, boundary_high and w are approximations of the real boundaries | 
|---|
| 111 | // and v (the input number). They are guaranteed to be precise up to one unit. | 
|---|
| 112 | // In fact the error is guaranteed to be strictly less than one unit. | 
|---|
| 113 | // | 
|---|
| 114 | // Anything that lies outside the unsafe interval is guaranteed not to round | 
|---|
| 115 | // to v when read again. | 
|---|
| 116 | // Anything that lies inside the safe interval is guaranteed to round to v | 
|---|
| 117 | // when read again. | 
|---|
| 118 | // If the number inside the buffer lies inside the unsafe interval but not | 
|---|
| 119 | // inside the safe interval then we simply do not know and bail out (returning | 
|---|
| 120 | // false). | 
|---|
| 121 | // | 
|---|
| 122 | // Similarly we have to take into account the imprecision of 'w' when finding | 
|---|
| 123 | // the closest representation of 'w'. If we have two potential | 
|---|
| 124 | // representations, and one is closer to both w_low and w_high, then we know | 
|---|
| 125 | // it is closer to the actual value v. | 
|---|
| 126 | // | 
|---|
| 127 | // By generating the digits of too_high we got the largest (closest to | 
|---|
| 128 | // too_high) buffer that is still in the unsafe interval. In the case where | 
|---|
| 129 | // w_high < buffer < too_high we try to decrement the buffer. | 
|---|
| 130 | // This way the buffer approaches (rounds towards) w. | 
|---|
| 131 | // There are 3 conditions that stop the decrementation process: | 
|---|
| 132 | //   1) the buffer is already below w_high | 
|---|
| 133 | //   2) decrementing the buffer would make it leave the unsafe interval | 
|---|
| 134 | //   3) decrementing the buffer would yield a number below w_high and farther | 
|---|
| 135 | //      away than the current number. In other words: | 
|---|
| 136 | //              (buffer{-1} < w_high) && w_high - buffer{-1} > buffer - w_high | 
|---|
| 137 | // Instead of using the buffer directly we use its distance to too_high. | 
|---|
| 138 | // Conceptually rest ~= too_high - buffer | 
|---|
| 139 | // We need to do the following tests in this order to avoid over- and | 
|---|
| 140 | // underflows. | 
|---|
| 141 | ASSERT(rest <= unsafe_interval); | 
|---|
| 142 | while (rest < small_distance &&  // Negated condition 1 | 
|---|
| 143 | unsafe_interval - rest >= ten_kappa &&  // Negated condition 2 | 
|---|
| 144 | (rest + ten_kappa < small_distance ||  // buffer{-1} > w_high | 
|---|
| 145 | small_distance - rest >= rest + ten_kappa - small_distance)) { | 
|---|
| 146 | buffer[length - 1]--; | 
|---|
| 147 | rest += ten_kappa; | 
|---|
| 148 | } | 
|---|
| 149 |  | 
|---|
| 150 | // We have approached w+ as much as possible. We now test if approaching w- | 
|---|
| 151 | // would require changing the buffer. If yes, then we have two possible | 
|---|
| 152 | // representations close to w, but we cannot decide which one is closer. | 
|---|
| 153 | if (rest < big_distance && | 
|---|
| 154 | unsafe_interval - rest >= ten_kappa && | 
|---|
| 155 | (rest + ten_kappa < big_distance || | 
|---|
| 156 | big_distance - rest > rest + ten_kappa - big_distance)) { | 
|---|
| 157 | return false; | 
|---|
| 158 | } | 
|---|
| 159 |  | 
|---|
| 160 | // Weeding test. | 
|---|
| 161 | //   The safe interval is [too_low + 2 ulp; too_high - 2 ulp] | 
|---|
| 162 | //   Since too_low = too_high - unsafe_interval this is equivalent to | 
|---|
| 163 | //      [too_high - unsafe_interval + 4 ulp; too_high - 2 ulp] | 
|---|
| 164 | //   Conceptually we have: rest ~= too_high - buffer | 
|---|
| 165 | return (2 * unit <= rest) && (rest <= unsafe_interval - 4 * unit); | 
|---|
| 166 | } | 
|---|
| 167 |  | 
|---|
| 168 |  | 
|---|
| 169 | // Rounds the buffer upwards if the result is closer to v by possibly adding | 
|---|
| 170 | // 1 to the buffer. If the precision of the calculation is not sufficient to | 
|---|
| 171 | // round correctly, return false. | 
|---|
| 172 | // The rounding might shift the whole buffer in which case the kappa is | 
|---|
| 173 | // adjusted. For example "99", kappa = 3 might become "10", kappa = 4. | 
|---|
| 174 | // | 
|---|
| 175 | // If 2*rest > ten_kappa then the buffer needs to be round up. | 
|---|
| 176 | // rest can have an error of +/- 1 unit. This function accounts for the | 
|---|
| 177 | // imprecision and returns false, if the rounding direction cannot be | 
|---|
| 178 | // unambiguously determined. | 
|---|
| 179 | // | 
|---|
| 180 | // Precondition: rest < ten_kappa. | 
|---|
| 181 | static bool RoundWeedCounted(Vector<char> buffer, | 
|---|
| 182 | int length, | 
|---|
| 183 | uint64_t rest, | 
|---|
| 184 | uint64_t ten_kappa, | 
|---|
| 185 | uint64_t unit, | 
|---|
| 186 | int* kappa) { | 
|---|
| 187 | ASSERT(rest < ten_kappa); | 
|---|
| 188 | // The following tests are done in a specific order to avoid overflows. They | 
|---|
| 189 | // will work correctly with any uint64 values of rest < ten_kappa and unit. | 
|---|
| 190 | // | 
|---|
| 191 | // If the unit is too big, then we don't know which way to round. For example | 
|---|
| 192 | // a unit of 50 means that the real number lies within rest +/- 50. If | 
|---|
| 193 | // 10^kappa == 40 then there is no way to tell which way to round. | 
|---|
| 194 | if (unit >= ten_kappa) return false; | 
|---|
| 195 | // Even if unit is just half the size of 10^kappa we are already completely | 
|---|
| 196 | // lost. (And after the previous test we know that the expression will not | 
|---|
| 197 | // over/underflow.) | 
|---|
| 198 | if (ten_kappa - unit <= unit) return false; | 
|---|
| 199 | // If 2 * (rest + unit) <= 10^kappa we can safely round down. | 
|---|
| 200 | if ((ten_kappa - rest > rest) && (ten_kappa - 2 * rest >= 2 * unit)) { | 
|---|
| 201 | return true; | 
|---|
| 202 | } | 
|---|
| 203 | // If 2 * (rest - unit) >= 10^kappa, then we can safely round up. | 
|---|
| 204 | if ((rest > unit) && (ten_kappa - (rest - unit) <= (rest - unit))) { | 
|---|
| 205 | // Increment the last digit recursively until we find a non '9' digit. | 
|---|
| 206 | buffer[length - 1]++; | 
|---|
| 207 | for (int i = length - 1; i > 0; --i) { | 
|---|
| 208 | if (buffer[i] != '0' + 10) break; | 
|---|
| 209 | buffer[i] = '0'; | 
|---|
| 210 | buffer[i - 1]++; | 
|---|
| 211 | } | 
|---|
| 212 | // If the first digit is now '0'+ 10 we had a buffer with all '9's. With the | 
|---|
| 213 | // exception of the first digit all digits are now '0'. Simply switch the | 
|---|
| 214 | // first digit to '1' and adjust the kappa. Example: "99" becomes "10" and | 
|---|
| 215 | // the power (the kappa) is increased. | 
|---|
| 216 | if (buffer[0] == '0' + 10) { | 
|---|
| 217 | buffer[0] = '1'; | 
|---|
| 218 | (*kappa) += 1; | 
|---|
| 219 | } | 
|---|
| 220 | return true; | 
|---|
| 221 | } | 
|---|
| 222 | return false; | 
|---|
| 223 | } | 
|---|
| 224 |  | 
|---|
| 225 | // Returns the biggest power of ten that is less than or equal to the given | 
|---|
| 226 | // number. We furthermore receive the maximum number of bits 'number' has. | 
|---|
| 227 | // | 
|---|
| 228 | // Returns power == 10^(exponent_plus_one-1) such that | 
|---|
| 229 | //    power <= number < power * 10. | 
|---|
| 230 | // If number_bits == 0 then 0^(0-1) is returned. | 
|---|
| 231 | // The number of bits must be <= 32. | 
|---|
| 232 | // Precondition: number < (1 << (number_bits + 1)). | 
|---|
| 233 |  | 
|---|
| 234 | // Inspired by the method for finding an integer log base 10 from here: | 
|---|
| 235 | // http://graphics.stanford.edu/~seander/bithacks.html#IntegerLog10 | 
|---|
| 236 | static unsigned int const kSmallPowersOfTen[] = | 
|---|
| 237 | {0, 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, | 
|---|
| 238 | 1000000000}; | 
|---|
| 239 |  | 
|---|
| 240 | static void BiggestPowerTen(uint32_t number, | 
|---|
| 241 | int number_bits, | 
|---|
| 242 | uint32_t* power, | 
|---|
| 243 | int* exponent_plus_one) { | 
|---|
| 244 | ASSERT(number < (1u << (number_bits + 1))); | 
|---|
| 245 | // 1233/4096 is approximately 1/lg(10). | 
|---|
| 246 | int exponent_plus_one_guess = ((number_bits + 1) * 1233 >> 12); | 
|---|
| 247 | // We increment to skip over the first entry in the kPowersOf10 table. | 
|---|
| 248 | // Note: kPowersOf10[i] == 10^(i-1). | 
|---|
| 249 | exponent_plus_one_guess++; | 
|---|
| 250 | // We don't have any guarantees that 2^number_bits <= number. | 
|---|
| 251 | if (number < kSmallPowersOfTen[exponent_plus_one_guess] && exponent_plus_one_guess > 0) { | 
|---|
| 252 | exponent_plus_one_guess--; | 
|---|
| 253 | } | 
|---|
| 254 | *power = kSmallPowersOfTen[exponent_plus_one_guess]; | 
|---|
| 255 | *exponent_plus_one = exponent_plus_one_guess; | 
|---|
| 256 | } | 
|---|
| 257 |  | 
|---|
| 258 | // Generates the digits of input number w. | 
|---|
| 259 | // w is a floating-point number (DiyFp), consisting of a significand and an | 
|---|
| 260 | // exponent. Its exponent is bounded by kMinimalTargetExponent and | 
|---|
| 261 | // kMaximalTargetExponent. | 
|---|
| 262 | //       Hence -60 <= w.e() <= -32. | 
|---|
| 263 | // | 
|---|
| 264 | // Returns false if it fails, in which case the generated digits in the buffer | 
|---|
| 265 | // should not be used. | 
|---|
| 266 | // Preconditions: | 
|---|
| 267 | //  * low, w and high are correct up to 1 ulp (unit in the last place). That | 
|---|
| 268 | //    is, their error must be less than a unit of their last digits. | 
|---|
| 269 | //  * low.e() == w.e() == high.e() | 
|---|
| 270 | //  * low < w < high, and taking into account their error: low~ <= high~ | 
|---|
| 271 | //  * kMinimalTargetExponent <= w.e() <= kMaximalTargetExponent | 
|---|
| 272 | // Postconditions: returns false if procedure fails. | 
|---|
| 273 | //   otherwise: | 
|---|
| 274 | //     * buffer is not null-terminated, but len contains the number of digits. | 
|---|
| 275 | //     * buffer contains the shortest possible decimal digit-sequence | 
|---|
| 276 | //       such that LOW < buffer * 10^kappa < HIGH, where LOW and HIGH are the | 
|---|
| 277 | //       correct values of low and high (without their error). | 
|---|
| 278 | //     * if more than one decimal representation gives the minimal number of | 
|---|
| 279 | //       decimal digits then the one closest to W (where W is the correct value | 
|---|
| 280 | //       of w) is chosen. | 
|---|
| 281 | // Remark: this procedure takes into account the imprecision of its input | 
|---|
| 282 | //   numbers. If the precision is not enough to guarantee all the postconditions | 
|---|
| 283 | //   then false is returned. This usually happens rarely (~0.5%). | 
|---|
| 284 | // | 
|---|
| 285 | // Say, for the sake of example, that | 
|---|
| 286 | //   w.e() == -48, and w.f() == 0x1234567890abcdef | 
|---|
| 287 | // w's value can be computed by w.f() * 2^w.e() | 
|---|
| 288 | // We can obtain w's integral digits by simply shifting w.f() by -w.e(). | 
|---|
| 289 | //  -> w's integral part is 0x1234 | 
|---|
| 290 | //  w's fractional part is therefore 0x567890abcdef. | 
|---|
| 291 | // Printing w's integral part is easy (simply print 0x1234 in decimal). | 
|---|
| 292 | // In order to print its fraction we repeatedly multiply the fraction by 10 and | 
|---|
| 293 | // get each digit. Example the first digit after the point would be computed by | 
|---|
| 294 | //   (0x567890abcdef * 10) >> 48. -> 3 | 
|---|
| 295 | // The whole thing becomes slightly more complicated because we want to stop | 
|---|
| 296 | // once we have enough digits. That is, once the digits inside the buffer | 
|---|
| 297 | // represent 'w' we can stop. Everything inside the interval low - high | 
|---|
| 298 | // represents w. However we have to pay attention to low, high and w's | 
|---|
| 299 | // imprecision. | 
|---|
| 300 | static bool DigitGen(DiyFp low, | 
|---|
| 301 | DiyFp w, | 
|---|
| 302 | DiyFp high, | 
|---|
| 303 | Vector<char> buffer, | 
|---|
| 304 | int* length, | 
|---|
| 305 | int* kappa) { | 
|---|
| 306 | ASSERT(low.e() == w.e() && w.e() == high.e()); | 
|---|
| 307 | ASSERT(low.f() + 1 <= high.f() - 1); | 
|---|
| 308 | ASSERT(kMinimalTargetExponent <= w.e() && w.e() <= kMaximalTargetExponent); | 
|---|
| 309 | // low, w and high are imprecise, but by less than one ulp (unit in the last | 
|---|
| 310 | // place). | 
|---|
| 311 | // If we remove (resp. add) 1 ulp from low (resp. high) we are certain that | 
|---|
| 312 | // the new numbers are outside of the interval we want the final | 
|---|
| 313 | // representation to lie in. | 
|---|
| 314 | // Inversely adding (resp. removing) 1 ulp from low (resp. high) would yield | 
|---|
| 315 | // numbers that are certain to lie in the interval. We will use this fact | 
|---|
| 316 | // later on. | 
|---|
| 317 | // We will now start by generating the digits within the uncertain | 
|---|
| 318 | // interval. Later we will weed out representations that lie outside the safe | 
|---|
| 319 | // interval and thus _might_ lie outside the correct interval. | 
|---|
| 320 | uint64_t unit = 1; | 
|---|
| 321 | DiyFp too_low = DiyFp(low.f() - unit, low.e()); | 
|---|
| 322 | DiyFp too_high = DiyFp(high.f() + unit, high.e()); | 
|---|
| 323 | // too_low and too_high are guaranteed to lie outside the interval we want the | 
|---|
| 324 | // generated number in. | 
|---|
| 325 | DiyFp unsafe_interval = DiyFp::Minus(too_high, too_low); | 
|---|
| 326 | // We now cut the input number into two parts: the integral digits and the | 
|---|
| 327 | // fractionals. We will not write any decimal separator though, but adapt | 
|---|
| 328 | // kappa instead. | 
|---|
| 329 | // Reminder: we are currently computing the digits (stored inside the buffer) | 
|---|
| 330 | // such that:   too_low < buffer * 10^kappa < too_high | 
|---|
| 331 | // We use too_high for the digit_generation and stop as soon as possible. | 
|---|
| 332 | // If we stop early we effectively round down. | 
|---|
| 333 | DiyFp one = DiyFp(static_cast<uint64_t>(1) << -w.e(), w.e()); | 
|---|
| 334 | // Division by one is a shift. | 
|---|
| 335 | uint32_t integrals = static_cast<uint32_t>(too_high.f() >> -one.e()); | 
|---|
| 336 | // Modulo by one is an and. | 
|---|
| 337 | uint64_t fractionals = too_high.f() & (one.f() - 1); | 
|---|
| 338 | uint32_t divisor; | 
|---|
| 339 | int divisor_exponent_plus_one; | 
|---|
| 340 | BiggestPowerTen(integrals, DiyFp::kSignificandSize - (-one.e()), | 
|---|
| 341 | &divisor, &divisor_exponent_plus_one); | 
|---|
| 342 | *kappa = divisor_exponent_plus_one; | 
|---|
| 343 | *length = 0; | 
|---|
| 344 | // Loop invariant: buffer = too_high / 10^kappa  (integer division) | 
|---|
| 345 | // The invariant holds for the first iteration: kappa has been initialized | 
|---|
| 346 | // with the divisor exponent + 1. And the divisor is the biggest power of ten | 
|---|
| 347 | // that is smaller than integrals. | 
|---|
| 348 | while (*kappa > 0) { | 
|---|
| 349 | int digit = integrals / divisor; | 
|---|
| 350 | ASSERT(digit <= 9); | 
|---|
| 351 | buffer[*length] = static_cast<char>('0' + digit); | 
|---|
| 352 | (*length)++; | 
|---|
| 353 | integrals %= divisor; | 
|---|
| 354 | (*kappa)--; | 
|---|
| 355 | // Note that kappa now equals the exponent of the divisor and that the | 
|---|
| 356 | // invariant thus holds again. | 
|---|
| 357 | uint64_t rest = | 
|---|
| 358 | (static_cast<uint64_t>(integrals) << -one.e()) + fractionals; | 
|---|
| 359 | // Invariant: too_high = buffer * 10^kappa + DiyFp(rest, one.e()) | 
|---|
| 360 | // Reminder: unsafe_interval.e() == one.e() | 
|---|
| 361 | if (rest < unsafe_interval.f()) { | 
|---|
| 362 | // Rounding down (by not emitting the remaining digits) yields a number | 
|---|
| 363 | // that lies within the unsafe interval. | 
|---|
| 364 | return RoundWeed(buffer, *length, DiyFp::Minus(too_high, w).f(), | 
|---|
| 365 | unsafe_interval.f(), rest, | 
|---|
| 366 | static_cast<uint64_t>(divisor) << -one.e(), unit); | 
|---|
| 367 | } | 
|---|
| 368 | divisor /= 10; | 
|---|
| 369 | } | 
|---|
| 370 |  | 
|---|
| 371 | // The integrals have been generated. We are at the point of the decimal | 
|---|
| 372 | // separator. In the following loop we simply multiply the remaining digits by | 
|---|
| 373 | // 10 and divide by one. We just need to pay attention to multiply associated | 
|---|
| 374 | // data (like the interval or 'unit'), too. | 
|---|
| 375 | // Note that the multiplication by 10 does not overflow, because w.e >= -60 | 
|---|
| 376 | // and thus one.e >= -60. | 
|---|
| 377 | ASSERT(one.e() >= -60); | 
|---|
| 378 | ASSERT(fractionals < one.f()); | 
|---|
| 379 | ASSERT(UINT64_2PART_C(0xFFFFFFFF, FFFFFFFF) / 10 >= one.f()); | 
|---|
| 380 | for (;;) { | 
|---|
| 381 | fractionals *= 10; | 
|---|
| 382 | unit *= 10; | 
|---|
| 383 | unsafe_interval.set_f(unsafe_interval.f() * 10); | 
|---|
| 384 | // Integer division by one. | 
|---|
| 385 | int digit = static_cast<int>(fractionals >> -one.e()); | 
|---|
| 386 | ASSERT(digit <= 9); | 
|---|
| 387 | buffer[*length] = static_cast<char>('0' + digit); | 
|---|
| 388 | (*length)++; | 
|---|
| 389 | fractionals &= one.f() - 1;  // Modulo by one. | 
|---|
| 390 | (*kappa)--; | 
|---|
| 391 | if (fractionals < unsafe_interval.f()) { | 
|---|
| 392 | return RoundWeed(buffer, *length, DiyFp::Minus(too_high, w).f() * unit, | 
|---|
| 393 | unsafe_interval.f(), fractionals, one.f(), unit); | 
|---|
| 394 | } | 
|---|
| 395 | } | 
|---|
| 396 | } | 
|---|
| 397 |  | 
|---|
| 398 |  | 
|---|
| 399 |  | 
|---|
| 400 | // Generates (at most) requested_digits digits of input number w. | 
|---|
| 401 | // w is a floating-point number (DiyFp), consisting of a significand and an | 
|---|
| 402 | // exponent. Its exponent is bounded by kMinimalTargetExponent and | 
|---|
| 403 | // kMaximalTargetExponent. | 
|---|
| 404 | //       Hence -60 <= w.e() <= -32. | 
|---|
| 405 | // | 
|---|
| 406 | // Returns false if it fails, in which case the generated digits in the buffer | 
|---|
| 407 | // should not be used. | 
|---|
| 408 | // Preconditions: | 
|---|
| 409 | //  * w is correct up to 1 ulp (unit in the last place). That | 
|---|
| 410 | //    is, its error must be strictly less than a unit of its last digit. | 
|---|
| 411 | //  * kMinimalTargetExponent <= w.e() <= kMaximalTargetExponent | 
|---|
| 412 | // | 
|---|
| 413 | // Postconditions: returns false if procedure fails. | 
|---|
| 414 | //   otherwise: | 
|---|
| 415 | //     * buffer is not null-terminated, but length contains the number of | 
|---|
| 416 | //       digits. | 
|---|
| 417 | //     * the representation in buffer is the most precise representation of | 
|---|
| 418 | //       requested_digits digits. | 
|---|
| 419 | //     * buffer contains at most requested_digits digits of w. If there are less | 
|---|
| 420 | //       than requested_digits digits then some trailing '0's have been removed. | 
|---|
| 421 | //     * kappa is such that | 
|---|
| 422 | //            w = buffer * 10^kappa + eps with |eps| < 10^kappa / 2. | 
|---|
| 423 | // | 
|---|
| 424 | // Remark: This procedure takes into account the imprecision of its input | 
|---|
| 425 | //   numbers. If the precision is not enough to guarantee all the postconditions | 
|---|
| 426 | //   then false is returned. This usually happens rarely, but the failure-rate | 
|---|
| 427 | //   increases with higher requested_digits. | 
|---|
| 428 | static bool DigitGenCounted(DiyFp w, | 
|---|
| 429 | int requested_digits, | 
|---|
| 430 | Vector<char> buffer, | 
|---|
| 431 | int* length, | 
|---|
| 432 | int* kappa) { | 
|---|
| 433 | ASSERT(kMinimalTargetExponent <= w.e() && w.e() <= kMaximalTargetExponent); | 
|---|
| 434 | ASSERT(kMinimalTargetExponent >= -60); | 
|---|
| 435 | ASSERT(kMaximalTargetExponent <= -32); | 
|---|
| 436 | // w is assumed to have an error less than 1 unit. Whenever w is scaled we | 
|---|
| 437 | // also scale its error. | 
|---|
| 438 | uint64_t w_error = 1; | 
|---|
| 439 | // We cut the input number into two parts: the integral digits and the | 
|---|
| 440 | // fractional digits. We don't emit any decimal separator, but adapt kappa | 
|---|
| 441 | // instead. Example: instead of writing "1.2" we put "12" into the buffer and | 
|---|
| 442 | // increase kappa by 1. | 
|---|
| 443 | DiyFp one = DiyFp(static_cast<uint64_t>(1) << -w.e(), w.e()); | 
|---|
| 444 | // Division by one is a shift. | 
|---|
| 445 | uint32_t integrals = static_cast<uint32_t>(w.f() >> -one.e()); | 
|---|
| 446 | // Modulo by one is an and. | 
|---|
| 447 | uint64_t fractionals = w.f() & (one.f() - 1); | 
|---|
| 448 | uint32_t divisor; | 
|---|
| 449 | int divisor_exponent_plus_one; | 
|---|
| 450 | BiggestPowerTen(integrals, DiyFp::kSignificandSize - (-one.e()), | 
|---|
| 451 | &divisor, &divisor_exponent_plus_one); | 
|---|
| 452 | *kappa = divisor_exponent_plus_one; | 
|---|
| 453 | *length = 0; | 
|---|
| 454 |  | 
|---|
| 455 | // Loop invariant: buffer = w / 10^kappa  (integer division) | 
|---|
| 456 | // The invariant holds for the first iteration: kappa has been initialized | 
|---|
| 457 | // with the divisor exponent + 1. And the divisor is the biggest power of ten | 
|---|
| 458 | // that is smaller than 'integrals'. | 
|---|
| 459 | while (*kappa > 0) { | 
|---|
| 460 | int digit = integrals / divisor; | 
|---|
| 461 | ASSERT(digit <= 9); | 
|---|
| 462 | buffer[*length] = static_cast<char>('0' + digit); | 
|---|
| 463 | (*length)++; | 
|---|
| 464 | requested_digits--; | 
|---|
| 465 | integrals %= divisor; | 
|---|
| 466 | (*kappa)--; | 
|---|
| 467 | // Note that kappa now equals the exponent of the divisor and that the | 
|---|
| 468 | // invariant thus holds again. | 
|---|
| 469 | if (requested_digits == 0) break; | 
|---|
| 470 | divisor /= 10; | 
|---|
| 471 | } | 
|---|
| 472 |  | 
|---|
| 473 | if (requested_digits == 0) { | 
|---|
| 474 | uint64_t rest = | 
|---|
| 475 | (static_cast<uint64_t>(integrals) << -one.e()) + fractionals; | 
|---|
| 476 | return RoundWeedCounted(buffer, *length, rest, | 
|---|
| 477 | static_cast<uint64_t>(divisor) << -one.e(), w_error, | 
|---|
| 478 | kappa); | 
|---|
| 479 | } | 
|---|
| 480 |  | 
|---|
| 481 | // The integrals have been generated. We are at the point of the decimal | 
|---|
| 482 | // separator. In the following loop we simply multiply the remaining digits by | 
|---|
| 483 | // 10 and divide by one. We just need to pay attention to multiply associated | 
|---|
| 484 | // data (the 'unit'), too. | 
|---|
| 485 | // Note that the multiplication by 10 does not overflow, because w.e >= -60 | 
|---|
| 486 | // and thus one.e >= -60. | 
|---|
| 487 | ASSERT(one.e() >= -60); | 
|---|
| 488 | ASSERT(fractionals < one.f()); | 
|---|
| 489 | ASSERT(UINT64_2PART_C(0xFFFFFFFF, FFFFFFFF) / 10 >= one.f()); | 
|---|
| 490 | while (requested_digits > 0 && fractionals > w_error) { | 
|---|
| 491 | fractionals *= 10; | 
|---|
| 492 | w_error *= 10; | 
|---|
| 493 | // Integer division by one. | 
|---|
| 494 | int digit = static_cast<int>(fractionals >> -one.e()); | 
|---|
| 495 | ASSERT(digit <= 9); | 
|---|
| 496 | buffer[*length] = static_cast<char>('0' + digit); | 
|---|
| 497 | (*length)++; | 
|---|
| 498 | requested_digits--; | 
|---|
| 499 | fractionals &= one.f() - 1;  // Modulo by one. | 
|---|
| 500 | (*kappa)--; | 
|---|
| 501 | } | 
|---|
| 502 | if (requested_digits != 0) return false; | 
|---|
| 503 | return RoundWeedCounted(buffer, *length, fractionals, one.f(), w_error, | 
|---|
| 504 | kappa); | 
|---|
| 505 | } | 
|---|
| 506 |  | 
|---|
| 507 |  | 
|---|
| 508 | // Provides a decimal representation of v. | 
|---|
| 509 | // Returns true if it succeeds, otherwise the result cannot be trusted. | 
|---|
| 510 | // There will be *length digits inside the buffer (not null-terminated). | 
|---|
| 511 | // If the function returns true then | 
|---|
| 512 | //        v == (double) (buffer * 10^decimal_exponent). | 
|---|
| 513 | // The digits in the buffer are the shortest representation possible: no | 
|---|
| 514 | // 0.09999999999999999 instead of 0.1. The shorter representation will even be | 
|---|
| 515 | // chosen even if the longer one would be closer to v. | 
|---|
| 516 | // The last digit will be closest to the actual v. That is, even if several | 
|---|
| 517 | // digits might correctly yield 'v' when read again, the closest will be | 
|---|
| 518 | // computed. | 
|---|
| 519 | static bool Grisu3(double v, | 
|---|
| 520 | FastDtoaMode mode, | 
|---|
| 521 | Vector<char> buffer, | 
|---|
| 522 | int* length, | 
|---|
| 523 | int* decimal_exponent) { | 
|---|
| 524 | DiyFp w = Double(v).AsNormalizedDiyFp(); | 
|---|
| 525 | // boundary_minus and boundary_plus are the boundaries between v and its | 
|---|
| 526 | // closest floating-point neighbors. Any number strictly between | 
|---|
| 527 | // boundary_minus and boundary_plus will round to v when convert to a double. | 
|---|
| 528 | // Grisu3 will never output representations that lie exactly on a boundary. | 
|---|
| 529 | DiyFp boundary_minus, boundary_plus; | 
|---|
| 530 | if (mode == FAST_DTOA_SHORTEST) { | 
|---|
| 531 | Double(v).NormalizedBoundaries(&boundary_minus, &boundary_plus); | 
|---|
| 532 | } else { | 
|---|
| 533 | ASSERT(mode == FAST_DTOA_SHORTEST_SINGLE); | 
|---|
| 534 | float single_v = static_cast<float>(v); | 
|---|
| 535 | Single(single_v).NormalizedBoundaries(&boundary_minus, &boundary_plus); | 
|---|
| 536 | } | 
|---|
| 537 | ASSERT(boundary_plus.e() == w.e()); | 
|---|
| 538 | DiyFp ten_mk;  // Cached power of ten: 10^-k | 
|---|
| 539 | int mk;        // -k | 
|---|
| 540 | int ten_mk_minimal_binary_exponent = | 
|---|
| 541 | kMinimalTargetExponent - (w.e() + DiyFp::kSignificandSize); | 
|---|
| 542 | int ten_mk_maximal_binary_exponent = | 
|---|
| 543 | kMaximalTargetExponent - (w.e() + DiyFp::kSignificandSize); | 
|---|
| 544 | PowersOfTenCache::GetCachedPowerForBinaryExponentRange( | 
|---|
| 545 | ten_mk_minimal_binary_exponent, | 
|---|
| 546 | ten_mk_maximal_binary_exponent, | 
|---|
| 547 | &ten_mk, &mk); | 
|---|
| 548 | ASSERT((kMinimalTargetExponent <= w.e() + ten_mk.e() + | 
|---|
| 549 | DiyFp::kSignificandSize) && | 
|---|
| 550 | (kMaximalTargetExponent >= w.e() + ten_mk.e() + | 
|---|
| 551 | DiyFp::kSignificandSize)); | 
|---|
| 552 | // Note that ten_mk is only an approximation of 10^-k. A DiyFp only contains a | 
|---|
| 553 | // 64 bit significand and ten_mk is thus only precise up to 64 bits. | 
|---|
| 554 |  | 
|---|
| 555 | // The DiyFp::Times procedure rounds its result, and ten_mk is approximated | 
|---|
| 556 | // too. The variable scaled_w (as well as scaled_boundary_minus/plus) are now | 
|---|
| 557 | // off by a small amount. | 
|---|
| 558 | // In fact: scaled_w - w*10^k < 1ulp (unit in the last place) of scaled_w. | 
|---|
| 559 | // In other words: let f = scaled_w.f() and e = scaled_w.e(), then | 
|---|
| 560 | //           (f-1) * 2^e < w*10^k < (f+1) * 2^e | 
|---|
| 561 | DiyFp scaled_w = DiyFp::Times(w, ten_mk); | 
|---|
| 562 | ASSERT(scaled_w.e() == | 
|---|
| 563 | boundary_plus.e() + ten_mk.e() + DiyFp::kSignificandSize); | 
|---|
| 564 | // In theory it would be possible to avoid some recomputations by computing | 
|---|
| 565 | // the difference between w and boundary_minus/plus (a power of 2) and to | 
|---|
| 566 | // compute scaled_boundary_minus/plus by subtracting/adding from | 
|---|
| 567 | // scaled_w. However the code becomes much less readable and the speed | 
|---|
| 568 | // enhancements are not terriffic. | 
|---|
| 569 | DiyFp scaled_boundary_minus = DiyFp::Times(boundary_minus, ten_mk); | 
|---|
| 570 | DiyFp scaled_boundary_plus  = DiyFp::Times(boundary_plus,  ten_mk); | 
|---|
| 571 |  | 
|---|
| 572 | // DigitGen will generate the digits of scaled_w. Therefore we have | 
|---|
| 573 | // v == (double) (scaled_w * 10^-mk). | 
|---|
| 574 | // Set decimal_exponent == -mk and pass it to DigitGen. If scaled_w is not an | 
|---|
| 575 | // integer than it will be updated. For instance if scaled_w == 1.23 then | 
|---|
| 576 | // the buffer will be filled with "123" und the decimal_exponent will be | 
|---|
| 577 | // decreased by 2. | 
|---|
| 578 | int kappa; | 
|---|
| 579 | bool result = DigitGen(scaled_boundary_minus, scaled_w, scaled_boundary_plus, | 
|---|
| 580 | buffer, length, &kappa); | 
|---|
| 581 | *decimal_exponent = -mk + kappa; | 
|---|
| 582 | return result; | 
|---|
| 583 | } | 
|---|
| 584 |  | 
|---|
| 585 |  | 
|---|
| 586 | // The "counted" version of grisu3 (see above) only generates requested_digits | 
|---|
| 587 | // number of digits. This version does not generate the shortest representation, | 
|---|
| 588 | // and with enough requested digits 0.1 will at some point print as 0.9999999... | 
|---|
| 589 | // Grisu3 is too imprecise for real halfway cases (1.5 will not work) and | 
|---|
| 590 | // therefore the rounding strategy for halfway cases is irrelevant. | 
|---|
| 591 | static bool Grisu3Counted(double v, | 
|---|
| 592 | int requested_digits, | 
|---|
| 593 | Vector<char> buffer, | 
|---|
| 594 | int* length, | 
|---|
| 595 | int* decimal_exponent) { | 
|---|
| 596 | DiyFp w = Double(v).AsNormalizedDiyFp(); | 
|---|
| 597 | DiyFp ten_mk;  // Cached power of ten: 10^-k | 
|---|
| 598 | int mk;        // -k | 
|---|
| 599 | int ten_mk_minimal_binary_exponent = | 
|---|
| 600 | kMinimalTargetExponent - (w.e() + DiyFp::kSignificandSize); | 
|---|
| 601 | int ten_mk_maximal_binary_exponent = | 
|---|
| 602 | kMaximalTargetExponent - (w.e() + DiyFp::kSignificandSize); | 
|---|
| 603 | PowersOfTenCache::GetCachedPowerForBinaryExponentRange( | 
|---|
| 604 | ten_mk_minimal_binary_exponent, | 
|---|
| 605 | ten_mk_maximal_binary_exponent, | 
|---|
| 606 | &ten_mk, &mk); | 
|---|
| 607 | ASSERT((kMinimalTargetExponent <= w.e() + ten_mk.e() + | 
|---|
| 608 | DiyFp::kSignificandSize) && | 
|---|
| 609 | (kMaximalTargetExponent >= w.e() + ten_mk.e() + | 
|---|
| 610 | DiyFp::kSignificandSize)); | 
|---|
| 611 | // Note that ten_mk is only an approximation of 10^-k. A DiyFp only contains a | 
|---|
| 612 | // 64 bit significand and ten_mk is thus only precise up to 64 bits. | 
|---|
| 613 |  | 
|---|
| 614 | // The DiyFp::Times procedure rounds its result, and ten_mk is approximated | 
|---|
| 615 | // too. The variable scaled_w (as well as scaled_boundary_minus/plus) are now | 
|---|
| 616 | // off by a small amount. | 
|---|
| 617 | // In fact: scaled_w - w*10^k < 1ulp (unit in the last place) of scaled_w. | 
|---|
| 618 | // In other words: let f = scaled_w.f() and e = scaled_w.e(), then | 
|---|
| 619 | //           (f-1) * 2^e < w*10^k < (f+1) * 2^e | 
|---|
| 620 | DiyFp scaled_w = DiyFp::Times(w, ten_mk); | 
|---|
| 621 |  | 
|---|
| 622 | // We now have (double) (scaled_w * 10^-mk). | 
|---|
| 623 | // DigitGen will generate the first requested_digits digits of scaled_w and | 
|---|
| 624 | // return together with a kappa such that scaled_w ~= buffer * 10^kappa. (It | 
|---|
| 625 | // will not always be exactly the same since DigitGenCounted only produces a | 
|---|
| 626 | // limited number of digits.) | 
|---|
| 627 | int kappa; | 
|---|
| 628 | bool result = DigitGenCounted(scaled_w, requested_digits, | 
|---|
| 629 | buffer, length, &kappa); | 
|---|
| 630 | *decimal_exponent = -mk + kappa; | 
|---|
| 631 | return result; | 
|---|
| 632 | } | 
|---|
| 633 |  | 
|---|
| 634 |  | 
|---|
| 635 | bool FastDtoa(double v, | 
|---|
| 636 | FastDtoaMode mode, | 
|---|
| 637 | int requested_digits, | 
|---|
| 638 | Vector<char> buffer, | 
|---|
| 639 | int* length, | 
|---|
| 640 | int* decimal_point) { | 
|---|
| 641 | ASSERT(v > 0); | 
|---|
| 642 | ASSERT(!Double(v).IsSpecial()); | 
|---|
| 643 |  | 
|---|
| 644 | bool result = false; | 
|---|
| 645 | int decimal_exponent = 0; | 
|---|
| 646 | switch (mode) { | 
|---|
| 647 | case FAST_DTOA_SHORTEST: | 
|---|
| 648 | case FAST_DTOA_SHORTEST_SINGLE: | 
|---|
| 649 | result = Grisu3(v, mode, buffer, length, &decimal_exponent); | 
|---|
| 650 | break; | 
|---|
| 651 | case FAST_DTOA_PRECISION: | 
|---|
| 652 | result = Grisu3Counted(v, requested_digits, | 
|---|
| 653 | buffer, length, &decimal_exponent); | 
|---|
| 654 | break; | 
|---|
| 655 | default: | 
|---|
| 656 | UNREACHABLE(); | 
|---|
| 657 | } | 
|---|
| 658 | if (result) { | 
|---|
| 659 | *decimal_point = *length + decimal_exponent; | 
|---|
| 660 | buffer[*length] = '\0'; | 
|---|
| 661 | } | 
|---|
| 662 | return result; | 
|---|
| 663 | } | 
|---|
| 664 |  | 
|---|
| 665 | }  // namespace double_conversion | 
|---|
| 666 |  | 
|---|