OLD | NEW |
1 // Copyright 2010 the V8 project authors. All rights reserved. | 1 // Copyright 2010 the V8 project authors. All rights reserved. |
2 // Redistribution and use in source and binary forms, with or without | 2 // Redistribution and use in source and binary forms, with or without |
3 // modification, are permitted provided that the following conditions are | 3 // modification, are permitted provided that the following conditions are |
4 // met: | 4 // met: |
5 // | 5 // |
6 // * Redistributions of source code must retain the above copyright | 6 // * Redistributions of source code must retain the above copyright |
7 // notice, this list of conditions and the following disclaimer. | 7 // notice, this list of conditions and the following disclaimer. |
8 // * Redistributions in binary form must reproduce the above | 8 // * Redistributions in binary form must reproduce the above |
9 // copyright notice, this list of conditions and the following | 9 // copyright notice, this list of conditions and the following |
10 // disclaimer in the documentation and/or other materials provided | 10 // disclaimer in the documentation and/or other materials provided |
11 // with the distribution. | 11 // with the distribution. |
12 // * Neither the name of Google Inc. nor the names of its | 12 // * Neither the name of Google Inc. nor the names of its |
13 // contributors may be used to endorse or promote products derived | 13 // contributors may be used to endorse or promote products derived |
14 // from this software without specific prior written permission. | 14 // from this software without specific prior written permission. |
15 // | 15 // |
16 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS | 16 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
17 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT | 17 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
18 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR | 18 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
19 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT | 19 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT |
20 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, | 20 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, |
21 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT | 21 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT |
22 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, | 22 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, |
23 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY | 23 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY |
24 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT | 24 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |
25 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE | 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. | 26 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
27 | 27 |
28 #include "bignum-dtoa.h" | 28 #include "bignum-dtoa.h" |
29 | 29 |
| 30 #include <math.h> |
30 #include "bignum.h" | 31 #include "bignum.h" |
31 #include "double.h" | 32 #include "double.h" |
32 #include <math.h> | |
33 | 33 |
34 namespace WTF { | 34 namespace WTF { |
35 | 35 |
36 namespace double_conversion { | 36 namespace double_conversion { |
37 | 37 |
38 static int NormalizedExponent(uint64_t significand, int exponent) { | 38 static int NormalizedExponent(uint64_t significand, int exponent) { |
39 ASSERT(significand != 0); | 39 ASSERT(significand != 0); |
40 while ((significand & Double::kHiddenBit) == 0) { | 40 while ((significand & Double::kHiddenBit) == 0) { |
41 significand = significand << 1; | 41 significand = significand << 1; |
42 exponent = exponent - 1; | 42 exponent = exponent - 1; |
| 43 } |
| 44 return exponent; |
| 45 } |
| 46 |
| 47 // Forward declarations: |
| 48 // Returns an estimation of k such that 10^(k-1) <= v < 10^k. |
| 49 static int EstimatePower(int exponent); |
| 50 // Computes v / 10^estimated_power exactly, as a ratio of two bignums, numerator |
| 51 // and denominator. |
| 52 static void InitialScaledStartValues(double v, |
| 53 int estimated_power, |
| 54 bool need_boundary_deltas, |
| 55 Bignum* numerator, |
| 56 Bignum* denominator, |
| 57 Bignum* delta_minus, |
| 58 Bignum* delta_plus); |
| 59 // Multiplies numerator/denominator so that its values lies in the range 1-10. |
| 60 // Returns decimal_point s.t. |
| 61 // v = numerator'/denominator' * 10^(decimal_point-1) |
| 62 // where numerator' and denominator' are the values of numerator and |
| 63 // denominator after the call to this function. |
| 64 static void FixupMultiply10(int estimated_power, |
| 65 bool is_even, |
| 66 int* decimal_point, |
| 67 Bignum* numerator, |
| 68 Bignum* denominator, |
| 69 Bignum* delta_minus, |
| 70 Bignum* delta_plus); |
| 71 // Generates digits from the left to the right and stops when the generated |
| 72 // digits yield the shortest decimal representation of v. |
| 73 static void GenerateShortestDigits(Bignum* numerator, |
| 74 Bignum* denominator, |
| 75 Bignum* delta_minus, |
| 76 Bignum* delta_plus, |
| 77 bool is_even, |
| 78 Vector<char> buffer, |
| 79 int* length); |
| 80 // Generates 'requested_digits' after the decimal point. |
| 81 static void BignumToFixed(int requested_digits, |
| 82 int* decimal_point, |
| 83 Bignum* numerator, |
| 84 Bignum* denominator, |
| 85 Vector<char>(buffer), |
| 86 int* length); |
| 87 // Generates 'count' digits of numerator/denominator. |
| 88 // Once 'count' digits have been produced rounds the result depending on the |
| 89 // remainder (remainders of exactly .5 round upwards). Might update the |
| 90 // decimal_point when rounding up (for example for 0.9999). |
| 91 static void GenerateCountedDigits(int count, |
| 92 int* decimal_point, |
| 93 Bignum* numerator, |
| 94 Bignum* denominator, |
| 95 Vector<char>(buffer), |
| 96 int* length); |
| 97 |
| 98 void BignumDtoa(double v, |
| 99 BignumDtoaMode mode, |
| 100 int requested_digits, |
| 101 Vector<char> buffer, |
| 102 int* length, |
| 103 int* decimal_point) { |
| 104 ASSERT(v > 0); |
| 105 ASSERT(!Double(v).IsSpecial()); |
| 106 uint64_t significand = Double(v).Significand(); |
| 107 bool is_even = (significand & 1) == 0; |
| 108 int exponent = Double(v).Exponent(); |
| 109 int normalized_exponent = NormalizedExponent(significand, exponent); |
| 110 // estimated_power might be too low by 1. |
| 111 int estimated_power = EstimatePower(normalized_exponent); |
| 112 |
| 113 // Shortcut for Fixed. |
| 114 // The requested digits correspond to the digits after the point. If the |
| 115 // number is much too small, then there is no need in trying to get any |
| 116 // digits. |
| 117 if (mode == BIGNUM_DTOA_FIXED && -estimated_power - 1 > requested_digits) { |
| 118 buffer[0] = '\0'; |
| 119 *length = 0; |
| 120 // Set decimal-point to -requested_digits. This is what Gay does. |
| 121 // Note that it should not have any effect anyways since the string is |
| 122 // empty. |
| 123 *decimal_point = -requested_digits; |
| 124 return; |
| 125 } |
| 126 |
| 127 Bignum numerator; |
| 128 Bignum denominator; |
| 129 Bignum delta_minus; |
| 130 Bignum delta_plus; |
| 131 // Make sure the bignum can grow large enough. The smallest double equals |
| 132 // 4e-324. In this case the denominator needs fewer than 324*4 binary digits. |
| 133 // The maximum double is 1.7976931348623157e308 which needs fewer than |
| 134 // 308*4 binary digits. |
| 135 ASSERT(Bignum::kMaxSignificantBits >= 324 * 4); |
| 136 bool need_boundary_deltas = (mode == BIGNUM_DTOA_SHORTEST); |
| 137 InitialScaledStartValues(v, estimated_power, need_boundary_deltas, &numerator, |
| 138 &denominator, &delta_minus, &delta_plus); |
| 139 // We now have v = (numerator / denominator) * 10^estimated_power. |
| 140 FixupMultiply10(estimated_power, is_even, decimal_point, &numerator, |
| 141 &denominator, &delta_minus, &delta_plus); |
| 142 // We now have v = (numerator / denominator) * 10^(decimal_point-1), and |
| 143 // 1 <= (numerator + delta_plus) / denominator < 10 |
| 144 switch (mode) { |
| 145 case BIGNUM_DTOA_SHORTEST: |
| 146 GenerateShortestDigits(&numerator, &denominator, &delta_minus, |
| 147 &delta_plus, is_even, buffer, length); |
| 148 break; |
| 149 case BIGNUM_DTOA_FIXED: |
| 150 BignumToFixed(requested_digits, decimal_point, &numerator, &denominator, |
| 151 buffer, length); |
| 152 break; |
| 153 case BIGNUM_DTOA_PRECISION: |
| 154 GenerateCountedDigits(requested_digits, decimal_point, &numerator, |
| 155 &denominator, buffer, length); |
| 156 break; |
| 157 default: |
| 158 UNREACHABLE(); |
| 159 } |
| 160 buffer[*length] = '\0'; |
| 161 } |
| 162 |
| 163 // The procedure starts generating digits from the left to the right and stops |
| 164 // when the generated digits yield the shortest decimal representation of v. A |
| 165 // decimal representation of v is a number lying closer to v than to any other |
| 166 // double, so it converts to v when read. |
| 167 // |
| 168 // This is true if d, the decimal representation, is between m- and m+, the |
| 169 // upper and lower boundaries. d must be strictly between them if !is_even. |
| 170 // m- := (numerator - delta_minus) / denominator |
| 171 // m+ := (numerator + delta_plus) / denominator |
| 172 // |
| 173 // Precondition: 0 <= (numerator+delta_plus) / denominator < 10. |
| 174 // If 1 <= (numerator+delta_plus) / denominator < 10 then no leading 0 digit |
| 175 // will be produced. This should be the standard precondition. |
| 176 static void GenerateShortestDigits(Bignum* numerator, |
| 177 Bignum* denominator, |
| 178 Bignum* delta_minus, |
| 179 Bignum* delta_plus, |
| 180 bool is_even, |
| 181 Vector<char> buffer, |
| 182 int* length) { |
| 183 // Small optimization: if delta_minus and delta_plus are the same just reuse |
| 184 // one of the two bignums. |
| 185 if (Bignum::Equal(*delta_minus, *delta_plus)) { |
| 186 delta_plus = delta_minus; |
| 187 } |
| 188 *length = 0; |
| 189 while (true) { |
| 190 uint16_t digit; |
| 191 digit = numerator->DivideModuloIntBignum(*denominator); |
| 192 ASSERT(digit <= 9); // digit is a uint16_t and therefore always positive. |
| 193 // digit = numerator / denominator (integer division). |
| 194 // numerator = numerator % denominator. |
| 195 buffer[(*length)++] = static_cast<char>(digit + '0'); |
| 196 |
| 197 // Can we stop already? |
| 198 // If the remainder of the division is less than the distance to the lower |
| 199 // boundary we can stop. In this case we simply round down (discarding the |
| 200 // remainder). |
| 201 // Similarly we test if we can round up (using the upper boundary). |
| 202 bool in_delta_room_minus; |
| 203 bool in_delta_room_plus; |
| 204 if (is_even) { |
| 205 in_delta_room_minus = Bignum::LessEqual(*numerator, *delta_minus); |
| 206 } else { |
| 207 in_delta_room_minus = Bignum::Less(*numerator, *delta_minus); |
| 208 } |
| 209 if (is_even) { |
| 210 in_delta_room_plus = |
| 211 Bignum::PlusCompare(*numerator, *delta_plus, *denominator) >= 0; |
| 212 } else { |
| 213 in_delta_room_plus = |
| 214 Bignum::PlusCompare(*numerator, *delta_plus, *denominator) > 0; |
| 215 } |
| 216 if (!in_delta_room_minus && !in_delta_room_plus) { |
| 217 // Prepare for next iteration. |
| 218 numerator->Times10(); |
| 219 delta_minus->Times10(); |
| 220 // We optimized delta_plus to be equal to delta_minus (if they share the |
| 221 // same value). So don't multiply delta_plus if they point to the same |
| 222 // object. |
| 223 if (delta_minus != delta_plus) { |
| 224 delta_plus->Times10(); |
| 225 } |
| 226 } else if (in_delta_room_minus && in_delta_room_plus) { |
| 227 // Let's see if 2*numerator < denominator. |
| 228 // If yes, then the next digit would be < 5 and we can round down. |
| 229 int compare = Bignum::PlusCompare(*numerator, *numerator, *denominator); |
| 230 if (compare < 0) { |
| 231 // Remaining digits are less than .5. -> Round down (== do nothing). |
| 232 } else if (compare > 0) { |
| 233 // Remaining digits are more than .5 of denominator. -> Round up. |
| 234 // Note that the last digit could not be a '9' as otherwise the whole |
| 235 // loop would have stopped earlier. |
| 236 // We still have an assert here in case the preconditions were not |
| 237 // satisfied. |
| 238 ASSERT(buffer[(*length) - 1] != '9'); |
| 239 buffer[(*length) - 1]++; |
| 240 } else { |
| 241 // Halfway case. |
| 242 // TODO(floitsch): need a way to solve half-way cases. |
| 243 // For now let's round towards even (since this is what Gay seems to |
| 244 // do). |
| 245 |
| 246 if ((buffer[(*length) - 1] - '0') % 2 == 0) { |
| 247 // Round down => Do nothing. |
| 248 } else { |
| 249 ASSERT(buffer[(*length) - 1] != '9'); |
| 250 buffer[(*length) - 1]++; |
43 } | 251 } |
44 return exponent; | 252 } |
45 } | 253 return; |
46 | 254 } else if (in_delta_room_minus) { |
47 | 255 // Round down (== do nothing). |
48 // Forward declarations: | 256 return; |
49 // Returns an estimation of k such that 10^(k-1) <= v < 10^k. | 257 } else { // in_delta_room_plus |
50 static int EstimatePower(int exponent); | 258 // Round up. |
51 // Computes v / 10^estimated_power exactly, as a ratio of two bignums, numer
ator | 259 // Note again that the last digit could not be '9' since this would have |
52 // and denominator. | 260 // stopped the loop earlier. |
53 static void InitialScaledStartValues(double v, | 261 // We still have an ASSERT here, in case the preconditions were not |
54 int estimated_power, | 262 // satisfied. |
55 bool need_boundary_deltas, | 263 ASSERT(buffer[(*length) - 1] != '9'); |
56 Bignum* numerator, | 264 buffer[(*length) - 1]++; |
57 Bignum* denominator, | 265 return; |
58 Bignum* delta_minus, | 266 } |
59 Bignum* delta_plus); | 267 } |
60 // Multiplies numerator/denominator so that its values lies in the range 1-1
0. | 268 } |
61 // Returns decimal_point s.t. | 269 |
62 // v = numerator'/denominator' * 10^(decimal_point-1) | 270 // Let v = numerator / denominator < 10. |
63 // where numerator' and denominator' are the values of numerator and | 271 // Then we generate 'count' digits of d = x.xxxxx... (without the decimal point) |
64 // denominator after the call to this function. | 272 // from left to right. Once 'count' digits have been produced we decide wether |
65 static void FixupMultiply10(int estimated_power, bool is_even, | 273 // to round up or down. Remainders of exactly .5 round upwards. Numbers such |
66 int* decimal_point, | 274 // as 9.999999 propagate a carry all the way, and change the |
67 Bignum* numerator, Bignum* denominator, | 275 // exponent (decimal_point), when rounding upwards. |
68 Bignum* delta_minus, Bignum* delta_plus); | 276 static void GenerateCountedDigits(int count, |
69 // Generates digits from the left to the right and stops when the generated | 277 int* decimal_point, |
70 // digits yield the shortest decimal representation of v. | 278 Bignum* numerator, |
71 static void GenerateShortestDigits(Bignum* numerator, Bignum* denominator, | 279 Bignum* denominator, |
72 Bignum* delta_minus, Bignum* delta_plus, | 280 Vector<char>(buffer), |
73 bool is_even, | 281 int* length) { |
74 Vector<char> buffer, int* length); | 282 ASSERT(count >= 0); |
75 // Generates 'requested_digits' after the decimal point. | 283 for (int i = 0; i < count - 1; ++i) { |
76 static void BignumToFixed(int requested_digits, int* decimal_point, | 284 uint16_t digit; |
77 Bignum* numerator, Bignum* denominator, | 285 digit = numerator->DivideModuloIntBignum(*denominator); |
78 Vector<char>(buffer), int* length); | 286 ASSERT(digit <= 9); // digit is a uint16_t and therefore always positive. |
79 // Generates 'count' digits of numerator/denominator. | 287 // digit = numerator / denominator (integer division). |
80 // Once 'count' digits have been produced rounds the result depending on the | 288 // numerator = numerator % denominator. |
81 // remainder (remainders of exactly .5 round upwards). Might update the | 289 buffer[i] = static_cast<char>(digit + '0'); |
82 // decimal_point when rounding up (for example for 0.9999). | 290 // Prepare for next iteration. |
83 static void GenerateCountedDigits(int count, int* decimal_point, | 291 numerator->Times10(); |
84 Bignum* numerator, Bignum* denominator, | 292 } |
85 Vector<char>(buffer), int* length); | 293 // Generate the last digit. |
86 | 294 uint16_t digit; |
87 | 295 digit = numerator->DivideModuloIntBignum(*denominator); |
88 void BignumDtoa(double v, BignumDtoaMode mode, int requested_digits, | 296 if (Bignum::PlusCompare(*numerator, *numerator, *denominator) >= 0) { |
89 Vector<char> buffer, int* length, int* decimal_point) { | 297 digit++; |
90 ASSERT(v > 0); | 298 } |
91 ASSERT(!Double(v).IsSpecial()); | 299 buffer[count - 1] = static_cast<char>(digit + '0'); |
92 uint64_t significand = Double(v).Significand(); | 300 // Correct bad digits (in case we had a sequence of '9's). Propagate the |
93 bool is_even = (significand & 1) == 0; | 301 // carry until we hat a non-'9' or til we reach the first digit. |
94 int exponent = Double(v).Exponent(); | 302 for (int i = count - 1; i > 0; --i) { |
95 int normalized_exponent = NormalizedExponent(significand, exponent); | 303 if (buffer[i] != '0' + 10) |
96 // estimated_power might be too low by 1. | 304 break; |
97 int estimated_power = EstimatePower(normalized_exponent); | 305 buffer[i] = '0'; |
98 | 306 buffer[i - 1]++; |
99 // Shortcut for Fixed. | 307 } |
100 // The requested digits correspond to the digits after the point. If the | 308 if (buffer[0] == '0' + 10) { |
101 // number is much too small, then there is no need in trying to get any | 309 // Propagate a carry past the top place. |
102 // digits. | 310 buffer[0] = '1'; |
103 if (mode == BIGNUM_DTOA_FIXED && -estimated_power - 1 > requested_digits
) { | 311 (*decimal_point)++; |
104 buffer[0] = '\0'; | 312 } |
105 *length = 0; | 313 *length = count; |
106 // Set decimal-point to -requested_digits. This is what Gay does. | 314 } |
107 // Note that it should not have any effect anyways since the string
is | 315 |
108 // empty. | 316 // Generates 'requested_digits' after the decimal point. It might omit |
109 *decimal_point = -requested_digits; | 317 // trailing '0's. If the input number is too small then no digits at all are |
110 return; | 318 // generated (ex.: 2 fixed digits for 0.00001). |
111 } | 319 // |
112 | 320 // Input verifies: 1 <= (numerator + delta) / denominator < 10. |
113 Bignum numerator; | 321 static void BignumToFixed(int requested_digits, |
114 Bignum denominator; | 322 int* decimal_point, |
115 Bignum delta_minus; | 323 Bignum* numerator, |
116 Bignum delta_plus; | 324 Bignum* denominator, |
117 // Make sure the bignum can grow large enough. The smallest double equal
s | 325 Vector<char>(buffer), |
118 // 4e-324. In this case the denominator needs fewer than 324*4 binary di
gits. | 326 int* length) { |
119 // The maximum double is 1.7976931348623157e308 which needs fewer than | 327 // Note that we have to look at more than just the requested_digits, since |
120 // 308*4 binary digits. | 328 // a number could be rounded up. Example: v=0.5 with requested_digits=0. |
121 ASSERT(Bignum::kMaxSignificantBits >= 324*4); | 329 // Even though the power of v equals 0 we can't just stop here. |
122 bool need_boundary_deltas = (mode == BIGNUM_DTOA_SHORTEST); | 330 if (-(*decimal_point) > requested_digits) { |
123 InitialScaledStartValues(v, estimated_power, need_boundary_deltas, | 331 // The number is definitively too small. |
124 &numerator, &denominator, | 332 // Ex: 0.001 with requested_digits == 1. |
125 &delta_minus, &delta_plus); | 333 // Set decimal-point to -requested_digits. This is what Gay does. |
126 // We now have v = (numerator / denominator) * 10^estimated_power. | 334 // Note that it should not have any effect anyways since the string is |
127 FixupMultiply10(estimated_power, is_even, decimal_point, | 335 // empty. |
128 &numerator, &denominator, | 336 *decimal_point = -requested_digits; |
129 &delta_minus, &delta_plus); | 337 *length = 0; |
130 // We now have v = (numerator / denominator) * 10^(decimal_point-1), and | 338 return; |
131 // 1 <= (numerator + delta_plus) / denominator < 10 | 339 } else if (-(*decimal_point) == requested_digits) { |
132 switch (mode) { | 340 // We only need to verify if the number rounds down or up. |
133 case BIGNUM_DTOA_SHORTEST: | 341 // Ex: 0.04 and 0.06 with requested_digits == 1. |
134 GenerateShortestDigits(&numerator, &denominator, | 342 ASSERT(*decimal_point == -requested_digits); |
135 &delta_minus, &delta_plus, | 343 // Initially the fraction lies in range (1, 10]. Multiply the denominator |
136 is_even, buffer, length); | 344 // by 10 so that we can compare more easily. |
137 break; | 345 denominator->Times10(); |
138 case BIGNUM_DTOA_FIXED: | 346 if (Bignum::PlusCompare(*numerator, *numerator, *denominator) >= 0) { |
139 BignumToFixed(requested_digits, decimal_point, | 347 // If the fraction is >= 0.5 then we have to include the rounded |
140 &numerator, &denominator, | 348 // digit. |
141 buffer, length); | 349 buffer[0] = '1'; |
142 break; | 350 *length = 1; |
143 case BIGNUM_DTOA_PRECISION: | 351 (*decimal_point)++; |
144 GenerateCountedDigits(requested_digits, decimal_point, | 352 } else { |
145 &numerator, &denominator, | 353 // Note that we caught most of similar cases earlier. |
146 buffer, length); | 354 *length = 0; |
147 break; | 355 } |
148 default: | 356 return; |
149 UNREACHABLE(); | 357 } else { |
150 } | 358 // The requested digits correspond to the digits after the point. |
151 buffer[*length] = '\0'; | 359 // The variable 'needed_digits' includes the digits before the point. |
152 } | 360 int needed_digits = (*decimal_point) + requested_digits; |
153 | 361 GenerateCountedDigits(needed_digits, decimal_point, numerator, denominator, |
154 | 362 buffer, length); |
155 // The procedure starts generating digits from the left to the right and sto
ps | 363 } |
156 // when the generated digits yield the shortest decimal representation of v.
A | 364 } |
157 // decimal representation of v is a number lying closer to v than to any oth
er | 365 |
158 // double, so it converts to v when read. | 366 // Returns an estimation of k such that 10^(k-1) <= v < 10^k where |
159 // | 367 // v = f * 2^exponent and 2^52 <= f < 2^53. |
160 // This is true if d, the decimal representation, is between m- and m+, the | 368 // v is hence a normalized double with the given exponent. The output is an |
161 // upper and lower boundaries. d must be strictly between them if !is_even. | 369 // approximation for the exponent of the decimal approimation .digits * 10^k. |
162 // m- := (numerator - delta_minus) / denominator | 370 // |
163 // m+ := (numerator + delta_plus) / denominator | 371 // The result might undershoot by 1 in which case 10^k <= v < 10^k+1. |
164 // | 372 // Note: this property holds for v's upper boundary m+ too. |
165 // Precondition: 0 <= (numerator+delta_plus) / denominator < 10. | 373 // 10^k <= m+ < 10^k+1. |
166 // If 1 <= (numerator+delta_plus) / denominator < 10 then no leading 0 dig
it | 374 // (see explanation below). |
167 // will be produced. This should be the standard precondition. | 375 // |
168 static void GenerateShortestDigits(Bignum* numerator, Bignum* denominator, | 376 // Examples: |
169 Bignum* delta_minus, Bignum* delta_plus, | 377 // EstimatePower(0) => 16 |
170 bool is_even, | 378 // EstimatePower(-52) => 0 |
171 Vector<char> buffer, int* length) { | 379 // |
172 // Small optimization: if delta_minus and delta_plus are the same just r
euse | 380 // Note: e >= 0 => EstimatedPower(e) > 0. No similar claim can be made for e<0. |
173 // one of the two bignums. | 381 static int EstimatePower(int exponent) { |
174 if (Bignum::Equal(*delta_minus, *delta_plus)) { | 382 // This function estimates log10 of v where v = f*2^e (with e == exponent). |
175 delta_plus = delta_minus; | 383 // Note that 10^floor(log10(v)) <= v, but v <= 10^ceil(log10(v)). |
176 } | 384 // Note that f is bounded by its container size. Let p = 53 (the double's |
177 *length = 0; | 385 // significand size). Then 2^(p-1) <= f < 2^p. |
178 while (true) { | 386 // |
179 uint16_t digit; | 387 // Given that log10(v) == log2(v)/log2(10) and e+(len(f)-1) is quite close |
180 digit = numerator->DivideModuloIntBignum(*denominator); | 388 // to log2(v) the function is simplified to (e+(len(f)-1)/log2(10)). |
181 ASSERT(digit <= 9); // digit is a uint16_t and therefore always pos
itive. | 389 // The computed number undershoots by less than 0.631 (when we compute log3 |
182 // digit = numerator / denominator (integer division). | 390 // and not log10). |
183 // numerator = numerator % denominator. | 391 // |
184 buffer[(*length)++] = static_cast<char>(digit + '0'); | 392 // Optimization: since we only need an approximated result this computation |
185 | 393 // can be performed on 64 bit integers. On x86/x64 architecture the speedup is |
186 // Can we stop already? | 394 // not really measurable, though. |
187 // If the remainder of the division is less than the distance to the
lower | 395 // |
188 // boundary we can stop. In this case we simply round down (discardi
ng the | 396 // Since we want to avoid overshooting we decrement by 1e10 so that |
189 // remainder). | 397 // floating-point imprecisions don't affect us. |
190 // Similarly we test if we can round up (using the upper boundary). | 398 // |
191 bool in_delta_room_minus; | 399 // Explanation for v's boundary m+: the computation takes advantage of |
192 bool in_delta_room_plus; | 400 // the fact that 2^(p-1) <= f < 2^p. Boundaries still satisfy this requirement |
193 if (is_even) { | 401 // (even for denormals where the delta can be much more important). |
194 in_delta_room_minus = Bignum::LessEqual(*numerator, *delta_minus
); | 402 |
195 } else { | 403 const double k1Log10 = 0.30102999566398114; // 1/lg(10) |
196 in_delta_room_minus = Bignum::Less(*numerator, *delta_minus); | 404 |
197 } | 405 // For doubles len(f) == 53 (don't forget the hidden bit). |
198 if (is_even) { | 406 const int kSignificandSize = 53; |
199 in_delta_room_plus = | 407 double estimate = ceil((exponent + kSignificandSize - 1) * k1Log10 - 1e-10); |
200 Bignum::PlusCompare(*numerator, *delta_plus, *denominator) >= 0; | 408 return static_cast<int>(estimate); |
201 } else { | 409 } |
202 in_delta_room_plus = | 410 |
203 Bignum::PlusCompare(*numerator, *delta_plus, *denominator) > 0; | 411 // See comments for InitialScaledStartValues. |
204 } | 412 static void InitialScaledStartValuesPositiveExponent(double v, |
205 if (!in_delta_room_minus && !in_delta_room_plus) { | 413 int estimated_power, |
206 // Prepare for next iteration. | 414 bool need_boundary_deltas, |
207 numerator->Times10(); | 415 Bignum* numerator, |
208 delta_minus->Times10(); | 416 Bignum* denominator, |
209 // We optimized delta_plus to be equal to delta_minus (if they s
hare the | 417 Bignum* delta_minus, |
210 // same value). So don't multiply delta_plus if they point to th
e same | 418 Bignum* delta_plus) { |
211 // object. | 419 // A positive exponent implies a positive power. |
212 if (delta_minus != delta_plus) { | 420 ASSERT(estimated_power >= 0); |
213 delta_plus->Times10(); | 421 // Since the estimated_power is positive we simply multiply the denominator |
214 } | 422 // by 10^estimated_power. |
215 } else if (in_delta_room_minus && in_delta_room_plus) { | 423 |
216 // Let's see if 2*numerator < denominator. | 424 // numerator = v. |
217 // If yes, then the next digit would be < 5 and we can round dow
n. | 425 numerator->AssignUInt64(Double(v).Significand()); |
218 int compare = Bignum::PlusCompare(*numerator, *numerator, *denom
inator); | 426 numerator->ShiftLeft(Double(v).Exponent()); |
219 if (compare < 0) { | 427 // denominator = 10^estimated_power. |
220 // Remaining digits are less than .5. -> Round down (== do n
othing). | 428 denominator->AssignPowerUInt16(10, estimated_power); |
221 } else if (compare > 0) { | 429 |
222 // Remaining digits are more than .5 of denominator. -> Roun
d up. | 430 if (need_boundary_deltas) { |
223 // Note that the last digit could not be a '9' as otherwise
the whole | 431 // Introduce a common denominator so that the deltas to the boundaries are |
224 // loop would have stopped earlier. | 432 // integers. |
225 // We still have an assert here in case the preconditions we
re not | 433 denominator->ShiftLeft(1); |
226 // satisfied. | 434 numerator->ShiftLeft(1); |
227 ASSERT(buffer[(*length) - 1] != '9'); | 435 // Let v = f * 2^e, then m+ - v = 1/2 * 2^e; With the common |
228 buffer[(*length) - 1]++; | 436 // denominator (of 2) delta_plus equals 2^e. |
229 } else { | 437 delta_plus->AssignUInt16(1); |
230 // Halfway case. | 438 delta_plus->ShiftLeft(Double(v).Exponent()); |
231 // TODO(floitsch): need a way to solve half-way cases. | 439 // Same for delta_minus (with adjustments below if f == 2^p-1). |
232 // For now let's round towards even (since this is what Ga
y seems to | 440 delta_minus->AssignUInt16(1); |
233 // do). | 441 delta_minus->ShiftLeft(Double(v).Exponent()); |
234 | 442 |
235 if ((buffer[(*length) - 1] - '0') % 2 == 0) { | 443 // If the significand (without the hidden bit) is 0, then the lower |
236 // Round down => Do nothing. | 444 // boundary is closer than just half a ulp (unit in the last place). |
237 } else { | 445 // There is only one exception: if the next lower number is a denormal then |
238 ASSERT(buffer[(*length) - 1] != '9'); | 446 // the distance is 1 ulp. This cannot be the case for exponent >= 0 (but we |
239 buffer[(*length) - 1]++; | 447 // have to test it in the other function where exponent < 0). |
240 } | 448 uint64_t v_bits = Double(v).AsUint64(); |
241 } | 449 if ((v_bits & Double::kSignificandMask) == 0) { |
242 return; | 450 // The lower boundary is closer at half the distance of "normal" numbers. |
243 } else if (in_delta_room_minus) { | 451 // Increase the common denominator and adapt all but the delta_minus. |
244 // Round down (== do nothing). | 452 denominator->ShiftLeft(1); // *2 |
245 return; | 453 numerator->ShiftLeft(1); // *2 |
246 } else { // in_delta_room_plus | 454 delta_plus->ShiftLeft(1); // *2 |
247 // Round up. | 455 } |
248 // Note again that the last digit could not be '9' since this wo
uld have | 456 } |
249 // stopped the loop earlier. | 457 } |
250 // We still have an ASSERT here, in case the preconditions were
not | 458 |
251 // satisfied. | 459 // See comments for InitialScaledStartValues |
252 ASSERT(buffer[(*length) -1] != '9'); | 460 static void InitialScaledStartValuesNegativeExponentPositivePower( |
253 buffer[(*length) - 1]++; | 461 double v, |
254 return; | 462 int estimated_power, |
255 } | 463 bool need_boundary_deltas, |
256 } | 464 Bignum* numerator, |
257 } | 465 Bignum* denominator, |
258 | 466 Bignum* delta_minus, |
259 | 467 Bignum* delta_plus) { |
260 // Let v = numerator / denominator < 10. | 468 uint64_t significand = Double(v).Significand(); |
261 // Then we generate 'count' digits of d = x.xxxxx... (without the decimal po
int) | 469 int exponent = Double(v).Exponent(); |
262 // from left to right. Once 'count' digits have been produced we decide weth
er | 470 // v = f * 2^e with e < 0, and with estimated_power >= 0. |
263 // to round up or down. Remainders of exactly .5 round upwards. Numbers such | 471 // This means that e is close to 0 (have a look at how estimated_power is |
264 // as 9.999999 propagate a carry all the way, and change the | 472 // computed). |
265 // exponent (decimal_point), when rounding upwards. | 473 |
266 static void GenerateCountedDigits(int count, int* decimal_point, | 474 // numerator = significand |
267 Bignum* numerator, Bignum* denominator, | 475 // since v = significand * 2^exponent this is equivalent to |
268 Vector<char>(buffer), int* length) { | 476 // numerator = v * / 2^-exponent |
269 ASSERT(count >= 0); | 477 numerator->AssignUInt64(significand); |
270 for (int i = 0; i < count - 1; ++i) { | 478 // denominator = 10^estimated_power * 2^-exponent (with exponent < 0) |
271 uint16_t digit; | 479 denominator->AssignPowerUInt16(10, estimated_power); |
272 digit = numerator->DivideModuloIntBignum(*denominator); | 480 denominator->ShiftLeft(-exponent); |
273 ASSERT(digit <= 9); // digit is a uint16_t and therefore always pos
itive. | 481 |
274 // digit = numerator / denominator (integer division). | 482 if (need_boundary_deltas) { |
275 // numerator = numerator % denominator. | 483 // Introduce a common denominator so that the deltas to the boundaries are |
276 buffer[i] = static_cast<char>(digit + '0'); | 484 // integers. |
277 // Prepare for next iteration. | 485 denominator->ShiftLeft(1); |
278 numerator->Times10(); | 486 numerator->ShiftLeft(1); |
279 } | 487 // Let v = f * 2^e, then m+ - v = 1/2 * 2^e; With the common |
280 // Generate the last digit. | 488 // denominator (of 2) delta_plus equals 2^e. |
281 uint16_t digit; | 489 // Given that the denominator already includes v's exponent the distance |
282 digit = numerator->DivideModuloIntBignum(*denominator); | 490 // to the boundaries is simply 1. |
283 if (Bignum::PlusCompare(*numerator, *numerator, *denominator) >= 0) { | 491 delta_plus->AssignUInt16(1); |
284 digit++; | 492 // Same for delta_minus (with adjustments below if f == 2^p-1). |
285 } | 493 delta_minus->AssignUInt16(1); |
286 buffer[count - 1] = static_cast<char>(digit + '0'); | 494 |
287 // Correct bad digits (in case we had a sequence of '9's). Propagate the | 495 // If the significand (without the hidden bit) is 0, then the lower |
288 // carry until we hat a non-'9' or til we reach the first digit. | 496 // boundary is closer than just one ulp (unit in the last place). |
289 for (int i = count - 1; i > 0; --i) { | 497 // There is only one exception: if the next lower number is a denormal |
290 if (buffer[i] != '0' + 10) break; | 498 // then the distance is 1 ulp. Since the exponent is close to zero |
291 buffer[i] = '0'; | 499 // (otherwise estimated_power would have been negative) this cannot happen |
292 buffer[i - 1]++; | 500 // here either. |
293 } | 501 uint64_t v_bits = Double(v).AsUint64(); |
294 if (buffer[0] == '0' + 10) { | 502 if ((v_bits & Double::kSignificandMask) == 0) { |
295 // Propagate a carry past the top place. | 503 // The lower boundary is closer at half the distance of "normal" numbers. |
296 buffer[0] = '1'; | 504 // Increase the denominator and adapt all but the delta_minus. |
297 (*decimal_point)++; | 505 denominator->ShiftLeft(1); // *2 |
298 } | 506 numerator->ShiftLeft(1); // *2 |
299 *length = count; | 507 delta_plus->ShiftLeft(1); // *2 |
300 } | 508 } |
301 | 509 } |
302 | 510 } |
303 // Generates 'requested_digits' after the decimal point. It might omit | 511 |
304 // trailing '0's. If the input number is too small then no digits at all are | 512 // See comments for InitialScaledStartValues |
305 // generated (ex.: 2 fixed digits for 0.00001). | 513 static void InitialScaledStartValuesNegativeExponentNegativePower( |
306 // | 514 double v, |
307 // Input verifies: 1 <= (numerator + delta) / denominator < 10. | 515 int estimated_power, |
308 static void BignumToFixed(int requested_digits, int* decimal_point, | 516 bool need_boundary_deltas, |
309 Bignum* numerator, Bignum* denominator, | 517 Bignum* numerator, |
310 Vector<char>(buffer), int* length) { | 518 Bignum* denominator, |
311 // Note that we have to look at more than just the requested_digits, sin
ce | 519 Bignum* delta_minus, |
312 // a number could be rounded up. Example: v=0.5 with requested_digits=0. | 520 Bignum* delta_plus) { |
313 // Even though the power of v equals 0 we can't just stop here. | 521 const uint64_t kMinimalNormalizedExponent = |
314 if (-(*decimal_point) > requested_digits) { | 522 UINT64_2PART_C(0x00100000, 00000000); |
315 // The number is definitively too small. | 523 uint64_t significand = Double(v).Significand(); |
316 // Ex: 0.001 with requested_digits == 1. | 524 int exponent = Double(v).Exponent(); |
317 // Set decimal-point to -requested_digits. This is what Gay does. | 525 // Instead of multiplying the denominator with 10^estimated_power we |
318 // Note that it should not have any effect anyways since the string
is | 526 // multiply all values (numerator and deltas) by 10^-estimated_power. |
319 // empty. | 527 |
320 *decimal_point = -requested_digits; | 528 // Use numerator as temporary container for power_ten. |
321 *length = 0; | 529 Bignum* power_ten = numerator; |
322 return; | 530 power_ten->AssignPowerUInt16(10, -estimated_power); |
323 } else if (-(*decimal_point) == requested_digits) { | 531 |
324 // We only need to verify if the number rounds down or up. | 532 if (need_boundary_deltas) { |
325 // Ex: 0.04 and 0.06 with requested_digits == 1. | 533 // Since power_ten == numerator we must make a copy of 10^estimated_power |
326 ASSERT(*decimal_point == -requested_digits); | 534 // before we complete the computation of the numerator. |
327 // Initially the fraction lies in range (1, 10]. Multiply the denomi
nator | 535 // delta_plus = delta_minus = 10^estimated_power |
328 // by 10 so that we can compare more easily. | 536 delta_plus->AssignBignum(*power_ten); |
329 denominator->Times10(); | 537 delta_minus->AssignBignum(*power_ten); |
330 if (Bignum::PlusCompare(*numerator, *numerator, *denominator) >= 0)
{ | 538 } |
331 // If the fraction is >= 0.5 then we have to include the rounded | 539 |
332 // digit. | 540 // numerator = significand * 2 * 10^-estimated_power |
333 buffer[0] = '1'; | 541 // since v = significand * 2^exponent this is equivalent to |
334 *length = 1; | 542 // numerator = v * 10^-estimated_power * 2 * 2^-exponent. |
335 (*decimal_point)++; | 543 // Remember: numerator has been abused as power_ten. So no need to assign it |
336 } else { | 544 // to itself. |
337 // Note that we caught most of similar cases earlier. | 545 ASSERT(numerator == power_ten); |
338 *length = 0; | 546 numerator->MultiplyByUInt64(significand); |
339 } | 547 |
340 return; | 548 // denominator = 2 * 2^-exponent with exponent < 0. |
341 } else { | 549 denominator->AssignUInt16(1); |
342 // The requested digits correspond to the digits after the point. | 550 denominator->ShiftLeft(-exponent); |
343 // The variable 'needed_digits' includes the digits before the point
. | 551 |
344 int needed_digits = (*decimal_point) + requested_digits; | 552 if (need_boundary_deltas) { |
345 GenerateCountedDigits(needed_digits, decimal_point, | 553 // Introduce a common denominator so that the deltas to the boundaries are |
346 numerator, denominator, | 554 // integers. |
347 buffer, length); | 555 numerator->ShiftLeft(1); |
348 } | 556 denominator->ShiftLeft(1); |
349 } | 557 // With this shift the boundaries have their correct value, since |
350 | 558 // delta_plus = 10^-estimated_power, and |
351 | 559 // delta_minus = 10^-estimated_power. |
352 // Returns an estimation of k such that 10^(k-1) <= v < 10^k where | 560 // These assignments have been done earlier. |
353 // v = f * 2^exponent and 2^52 <= f < 2^53. | 561 |
354 // v is hence a normalized double with the given exponent. The output is an | 562 // The special case where the lower boundary is twice as close. |
355 // approximation for the exponent of the decimal approimation .digits * 10^k
. | 563 // This time we have to look out for the exception too. |
356 // | 564 uint64_t v_bits = Double(v).AsUint64(); |
357 // The result might undershoot by 1 in which case 10^k <= v < 10^k+1. | 565 if ((v_bits & Double::kSignificandMask) == 0 && |
358 // Note: this property holds for v's upper boundary m+ too. | 566 // The only exception where a significand == 0 has its boundaries at |
359 // 10^k <= m+ < 10^k+1. | 567 // "normal" distances: |
360 // (see explanation below). | 568 (v_bits & Double::kExponentMask) != kMinimalNormalizedExponent) { |
361 // | 569 numerator->ShiftLeft(1); // *2 |
362 // Examples: | 570 denominator->ShiftLeft(1); // *2 |
363 // EstimatePower(0) => 16 | 571 delta_plus->ShiftLeft(1); // *2 |
364 // EstimatePower(-52) => 0 | 572 } |
365 // | 573 } |
366 // Note: e >= 0 => EstimatedPower(e) > 0. No similar claim can be made for e
<0. | 574 } |
367 static int EstimatePower(int exponent) { | 575 |
368 // This function estimates log10 of v where v = f*2^e (with e == exponen
t). | 576 // Let v = significand * 2^exponent. |
369 // Note that 10^floor(log10(v)) <= v, but v <= 10^ceil(log10(v)). | 577 // Computes v / 10^estimated_power exactly, as a ratio of two bignums, numerator |
370 // Note that f is bounded by its container size. Let p = 53 (the double'
s | 578 // and denominator. The functions GenerateShortestDigits and |
371 // significand size). Then 2^(p-1) <= f < 2^p. | 579 // GenerateCountedDigits will then convert this ratio to its decimal |
372 // | 580 // representation d, with the required accuracy. |
373 // Given that log10(v) == log2(v)/log2(10) and e+(len(f)-1) is quite clo
se | 581 // Then d * 10^estimated_power is the representation of v. |
374 // to log2(v) the function is simplified to (e+(len(f)-1)/log2(10)). | 582 // (Note: the fraction and the estimated_power might get adjusted before |
375 // The computed number undershoots by less than 0.631 (when we compute l
og3 | 583 // generating the decimal representation.) |
376 // and not log10). | 584 // |
377 // | 585 // The initial start values consist of: |
378 // Optimization: since we only need an approximated result this computat
ion | 586 // - a scaled numerator: s.t. numerator/denominator == v / 10^estimated_power. |
379 // can be performed on 64 bit integers. On x86/x64 architecture the spee
dup is | 587 // - a scaled (common) denominator. |
380 // not really measurable, though. | 588 // optionally (used by GenerateShortestDigits to decide if it has the shortest |
381 // | 589 // decimal converting back to v): |
382 // Since we want to avoid overshooting we decrement by 1e10 so that | 590 // - v - m-: the distance to the lower boundary. |
383 // floating-point imprecisions don't affect us. | 591 // - m+ - v: the distance to the upper boundary. |
384 // | 592 // |
385 // Explanation for v's boundary m+: the computation takes advantage of | 593 // v, m+, m-, and therefore v - m- and m+ - v all share the same denominator. |
386 // the fact that 2^(p-1) <= f < 2^p. Boundaries still satisfy this requi
rement | 594 // |
387 // (even for denormals where the delta can be much more important). | 595 // Let ep == estimated_power, then the returned values will satisfy: |
388 | 596 // v / 10^ep = numerator / denominator. |
389 const double k1Log10 = 0.30102999566398114; // 1/lg(10) | 597 // v's boundarys m- and m+: |
390 | 598 // m- / 10^ep == v / 10^ep - delta_minus / denominator |
391 // For doubles len(f) == 53 (don't forget the hidden bit). | 599 // m+ / 10^ep == v / 10^ep + delta_plus / denominator |
392 const int kSignificandSize = 53; | 600 // Or in other words: |
393 double estimate = ceil((exponent + kSignificandSize - 1) * k1Log10 - 1e-
10); | 601 // m- == v - delta_minus * 10^ep / denominator; |
394 return static_cast<int>(estimate); | 602 // m+ == v + delta_plus * 10^ep / denominator; |
395 } | 603 // |
396 | 604 // Since 10^(k-1) <= v < 10^k (with k == estimated_power) |
397 | 605 // or 10^k <= v < 10^(k+1) |
398 // See comments for InitialScaledStartValues. | 606 // we then have 0.1 <= numerator/denominator < 1 |
399 static void InitialScaledStartValuesPositiveExponent( | 607 // or 1 <= numerator/denominator < 10 |
400 double v, int estimated
_power, bool need_boundary_deltas, | 608 // |
401 Bignum* numerator, Bign
um* denominator, | 609 // It is then easy to kickstart the digit-generation routine. |
402 Bignum* delta_minus, Bi
gnum* delta_plus) { | 610 // |
403 // A positive exponent implies a positive power. | 611 // The boundary-deltas are only filled if need_boundary_deltas is set. |
404 ASSERT(estimated_power >= 0); | 612 static void InitialScaledStartValues(double v, |
405 // Since the estimated_power is positive we simply multiply the denomina
tor | 613 int estimated_power, |
406 // by 10^estimated_power. | 614 bool need_boundary_deltas, |
407 | 615 Bignum* numerator, |
408 // numerator = v. | 616 Bignum* denominator, |
409 numerator->AssignUInt64(Double(v).Significand()); | 617 Bignum* delta_minus, |
410 numerator->ShiftLeft(Double(v).Exponent()); | 618 Bignum* delta_plus) { |
411 // denominator = 10^estimated_power. | 619 if (Double(v).Exponent() >= 0) { |
412 denominator->AssignPowerUInt16(10, estimated_power); | 620 InitialScaledStartValuesPositiveExponent( |
413 | 621 v, estimated_power, need_boundary_deltas, numerator, denominator, |
414 if (need_boundary_deltas) { | 622 delta_minus, delta_plus); |
415 // Introduce a common denominator so that the deltas to the boundari
es are | 623 } else if (estimated_power >= 0) { |
416 // integers. | 624 InitialScaledStartValuesNegativeExponentPositivePower( |
417 denominator->ShiftLeft(1); | 625 v, estimated_power, need_boundary_deltas, numerator, denominator, |
418 numerator->ShiftLeft(1); | 626 delta_minus, delta_plus); |
419 // Let v = f * 2^e, then m+ - v = 1/2 * 2^e; With the common | 627 } else { |
420 // denominator (of 2) delta_plus equals 2^e. | 628 InitialScaledStartValuesNegativeExponentNegativePower( |
421 delta_plus->AssignUInt16(1); | 629 v, estimated_power, need_boundary_deltas, numerator, denominator, |
422 delta_plus->ShiftLeft(Double(v).Exponent()); | 630 delta_minus, delta_plus); |
423 // Same for delta_minus (with adjustments below if f == 2^p-1). | 631 } |
424 delta_minus->AssignUInt16(1); | 632 } |
425 delta_minus->ShiftLeft(Double(v).Exponent()); | 633 |
426 | 634 // This routine multiplies numerator/denominator so that its values lies in the |
427 // If the significand (without the hidden bit) is 0, then the lower | 635 // range 1-10. That is after a call to this function we have: |
428 // boundary is closer than just half a ulp (unit in the last place). | 636 // 1 <= (numerator + delta_plus) /denominator < 10. |
429 // There is only one exception: if the next lower number is a denorm
al then | 637 // Let numerator the input before modification and numerator' the argument |
430 // the distance is 1 ulp. This cannot be the case for exponent >= 0
(but we | 638 // after modification, then the output-parameter decimal_point is such that |
431 // have to test it in the other function where exponent < 0). | 639 // numerator / denominator * 10^estimated_power == |
432 uint64_t v_bits = Double(v).AsUint64(); | 640 // numerator' / denominator' * 10^(decimal_point - 1) |
433 if ((v_bits & Double::kSignificandMask) == 0) { | 641 // In some cases estimated_power was too low, and this is already the case. We |
434 // The lower boundary is closer at half the distance of "normal"
numbers. | 642 // then simply adjust the power so that 10^(k-1) <= v < 10^k (with k == |
435 // Increase the common denominator and adapt all but the delta_m
inus. | 643 // estimated_power) but do not touch the numerator or denominator. |
436 denominator->ShiftLeft(1); // *2 | 644 // Otherwise the routine multiplies the numerator and the deltas by 10. |
437 numerator->ShiftLeft(1); // *2 | 645 static void FixupMultiply10(int estimated_power, |
438 delta_plus->ShiftLeft(1); // *2 | 646 bool is_even, |
439 } | 647 int* decimal_point, |
440 } | 648 Bignum* numerator, |
441 } | 649 Bignum* denominator, |
442 | 650 Bignum* delta_minus, |
443 | 651 Bignum* delta_plus) { |
444 // See comments for InitialScaledStartValues | 652 bool in_range; |
445 static void InitialScaledStartValuesNegativeExponentPositivePower( | 653 if (is_even) { |
446 double v,
int estimated_power, bool need_boundary_deltas, | 654 // For IEEE doubles half-way cases (in decimal system numbers ending with 5) |
447 Bignum* nu
merator, Bignum* denominator, | 655 // are rounded to the closest floating-point number with even significand. |
448 Bignum* de
lta_minus, Bignum* delta_plus) { | 656 in_range = Bignum::PlusCompare(*numerator, *delta_plus, *denominator) >= 0; |
449 uint64_t significand = Double(v).Significand(); | 657 } else { |
450 int exponent = Double(v).Exponent(); | 658 in_range = Bignum::PlusCompare(*numerator, *delta_plus, *denominator) > 0; |
451 // v = f * 2^e with e < 0, and with estimated_power >= 0. | 659 } |
452 // This means that e is close to 0 (have a look at how estimated_power i
s | 660 if (in_range) { |
453 // computed). | 661 // Since numerator + delta_plus >= denominator we already have |
454 | 662 // 1 <= numerator/denominator < 10. Simply update the estimated_power. |
455 // numerator = significand | 663 *decimal_point = estimated_power + 1; |
456 // since v = significand * 2^exponent this is equivalent to | 664 } else { |
457 // numerator = v * / 2^-exponent | 665 *decimal_point = estimated_power; |
458 numerator->AssignUInt64(significand); | 666 numerator->Times10(); |
459 // denominator = 10^estimated_power * 2^-exponent (with exponent < 0) | 667 if (Bignum::Equal(*delta_minus, *delta_plus)) { |
460 denominator->AssignPowerUInt16(10, estimated_power); | 668 delta_minus->Times10(); |
461 denominator->ShiftLeft(-exponent); | 669 delta_plus->AssignBignum(*delta_minus); |
462 | 670 } else { |
463 if (need_boundary_deltas) { | 671 delta_minus->Times10(); |
464 // Introduce a common denominator so that the deltas to the boundari
es are | 672 delta_plus->Times10(); |
465 // integers. | 673 } |
466 denominator->ShiftLeft(1); | 674 } |
467 numerator->ShiftLeft(1); | 675 } |
468 // Let v = f * 2^e, then m+ - v = 1/2 * 2^e; With the common | |
469 // denominator (of 2) delta_plus equals 2^e. | |
470 // Given that the denominator already includes v's exponent the dist
ance | |
471 // to the boundaries is simply 1. | |
472 delta_plus->AssignUInt16(1); | |
473 // Same for delta_minus (with adjustments below if f == 2^p-1). | |
474 delta_minus->AssignUInt16(1); | |
475 | |
476 // If the significand (without the hidden bit) is 0, then the lower | |
477 // boundary is closer than just one ulp (unit in the last place). | |
478 // There is only one exception: if the next lower number is a denorm
al | |
479 // then the distance is 1 ulp. Since the exponent is close to zero | |
480 // (otherwise estimated_power would have been negative) this cannot
happen | |
481 // here either. | |
482 uint64_t v_bits = Double(v).AsUint64(); | |
483 if ((v_bits & Double::kSignificandMask) == 0) { | |
484 // The lower boundary is closer at half the distance of "normal"
numbers. | |
485 // Increase the denominator and adapt all but the delta_minus. | |
486 denominator->ShiftLeft(1); // *2 | |
487 numerator->ShiftLeft(1); // *2 | |
488 delta_plus->ShiftLeft(1); // *2 | |
489 } | |
490 } | |
491 } | |
492 | |
493 | |
494 // See comments for InitialScaledStartValues | |
495 static void InitialScaledStartValuesNegativeExponentNegativePower( | |
496 double v,
int estimated_power, bool need_boundary_deltas, | |
497 Bignum* nu
merator, Bignum* denominator, | |
498 Bignum* de
lta_minus, Bignum* delta_plus) { | |
499 const uint64_t kMinimalNormalizedExponent = | |
500 UINT64_2PART_C(0x00100000, 00000000); | |
501 uint64_t significand = Double(v).Significand(); | |
502 int exponent = Double(v).Exponent(); | |
503 // Instead of multiplying the denominator with 10^estimated_power we | |
504 // multiply all values (numerator and deltas) by 10^-estimated_power. | |
505 | |
506 // Use numerator as temporary container for power_ten. | |
507 Bignum* power_ten = numerator; | |
508 power_ten->AssignPowerUInt16(10, -estimated_power); | |
509 | |
510 if (need_boundary_deltas) { | |
511 // Since power_ten == numerator we must make a copy of 10^estimated_
power | |
512 // before we complete the computation of the numerator. | |
513 // delta_plus = delta_minus = 10^estimated_power | |
514 delta_plus->AssignBignum(*power_ten); | |
515 delta_minus->AssignBignum(*power_ten); | |
516 } | |
517 | |
518 // numerator = significand * 2 * 10^-estimated_power | |
519 // since v = significand * 2^exponent this is equivalent to | |
520 // numerator = v * 10^-estimated_power * 2 * 2^-exponent. | |
521 // Remember: numerator has been abused as power_ten. So no need to assig
n it | |
522 // to itself. | |
523 ASSERT(numerator == power_ten); | |
524 numerator->MultiplyByUInt64(significand); | |
525 | |
526 // denominator = 2 * 2^-exponent with exponent < 0. | |
527 denominator->AssignUInt16(1); | |
528 denominator->ShiftLeft(-exponent); | |
529 | |
530 if (need_boundary_deltas) { | |
531 // Introduce a common denominator so that the deltas to the boundari
es are | |
532 // integers. | |
533 numerator->ShiftLeft(1); | |
534 denominator->ShiftLeft(1); | |
535 // With this shift the boundaries have their correct value, since | |
536 // delta_plus = 10^-estimated_power, and | |
537 // delta_minus = 10^-estimated_power. | |
538 // These assignments have been done earlier. | |
539 | |
540 // The special case where the lower boundary is twice as close. | |
541 // This time we have to look out for the exception too. | |
542 uint64_t v_bits = Double(v).AsUint64(); | |
543 if ((v_bits & Double::kSignificandMask) == 0 && | |
544 // The only exception where a significand == 0 has its boundarie
s at | |
545 // "normal" distances: | |
546 (v_bits & Double::kExponentMask) != kMinimalNormalizedExponent)
{ | |
547 numerator->ShiftLeft(1); // *2 | |
548 denominator->ShiftLeft(1); // *2 | |
549 delta_plus->ShiftLeft(1); // *2 | |
550 } | |
551 } | |
552 } | |
553 | |
554 | |
555 // Let v = significand * 2^exponent. | |
556 // Computes v / 10^estimated_power exactly, as a ratio of two bignums, numer
ator | |
557 // and denominator. The functions GenerateShortestDigits and | |
558 // GenerateCountedDigits will then convert this ratio to its decimal | |
559 // representation d, with the required accuracy. | |
560 // Then d * 10^estimated_power is the representation of v. | |
561 // (Note: the fraction and the estimated_power might get adjusted before | |
562 // generating the decimal representation.) | |
563 // | |
564 // The initial start values consist of: | |
565 // - a scaled numerator: s.t. numerator/denominator == v / 10^estimated_pow
er. | |
566 // - a scaled (common) denominator. | |
567 // optionally (used by GenerateShortestDigits to decide if it has the short
est | |
568 // decimal converting back to v): | |
569 // - v - m-: the distance to the lower boundary. | |
570 // - m+ - v: the distance to the upper boundary. | |
571 // | |
572 // v, m+, m-, and therefore v - m- and m+ - v all share the same denominator
. | |
573 // | |
574 // Let ep == estimated_power, then the returned values will satisfy: | |
575 // v / 10^ep = numerator / denominator. | |
576 // v's boundarys m- and m+: | |
577 // m- / 10^ep == v / 10^ep - delta_minus / denominator | |
578 // m+ / 10^ep == v / 10^ep + delta_plus / denominator | |
579 // Or in other words: | |
580 // m- == v - delta_minus * 10^ep / denominator; | |
581 // m+ == v + delta_plus * 10^ep / denominator; | |
582 // | |
583 // Since 10^(k-1) <= v < 10^k (with k == estimated_power) | |
584 // or 10^k <= v < 10^(k+1) | |
585 // we then have 0.1 <= numerator/denominator < 1 | |
586 // or 1 <= numerator/denominator < 10 | |
587 // | |
588 // It is then easy to kickstart the digit-generation routine. | |
589 // | |
590 // The boundary-deltas are only filled if need_boundary_deltas is set. | |
591 static void InitialScaledStartValues(double v, | |
592 int estimated_power, | |
593 bool need_boundary_deltas, | |
594 Bignum* numerator, | |
595 Bignum* denominator, | |
596 Bignum* delta_minus, | |
597 Bignum* delta_plus) { | |
598 if (Double(v).Exponent() >= 0) { | |
599 InitialScaledStartValuesPositiveExponent( | |
600 v, estimated_power, need_bo
undary_deltas, | |
601 numerator, denominator, del
ta_minus, delta_plus); | |
602 } else if (estimated_power >= 0) { | |
603 InitialScaledStartValuesNegativeExponentPositivePower( | |
604 v, estimated_p
ower, need_boundary_deltas, | |
605 numerator, den
ominator, delta_minus, delta_plus); | |
606 } else { | |
607 InitialScaledStartValuesNegativeExponentNegativePower( | |
608 v, estimated_p
ower, need_boundary_deltas, | |
609 numerator, den
ominator, delta_minus, delta_plus); | |
610 } | |
611 } | |
612 | |
613 | |
614 // This routine multiplies numerator/denominator so that its values lies in
the | |
615 // range 1-10. That is after a call to this function we have: | |
616 // 1 <= (numerator + delta_plus) /denominator < 10. | |
617 // Let numerator the input before modification and numerator' the argument | |
618 // after modification, then the output-parameter decimal_point is such that | |
619 // numerator / denominator * 10^estimated_power == | |
620 // numerator' / denominator' * 10^(decimal_point - 1) | |
621 // In some cases estimated_power was too low, and this is already the case.
We | |
622 // then simply adjust the power so that 10^(k-1) <= v < 10^k (with k == | |
623 // estimated_power) but do not touch the numerator or denominator. | |
624 // Otherwise the routine multiplies the numerator and the deltas by 10. | |
625 static void FixupMultiply10(int estimated_power, bool is_even, | |
626 int* decimal_point, | |
627 Bignum* numerator, Bignum* denominator, | |
628 Bignum* delta_minus, Bignum* delta_plus) { | |
629 bool in_range; | |
630 if (is_even) { | |
631 // For IEEE doubles half-way cases (in decimal system numbers ending
with 5) | |
632 // are rounded to the closest floating-point number with even signif
icand. | |
633 in_range = Bignum::PlusCompare(*numerator, *delta_plus, *denominator
) >= 0; | |
634 } else { | |
635 in_range = Bignum::PlusCompare(*numerator, *delta_plus, *denominator
) > 0; | |
636 } | |
637 if (in_range) { | |
638 // Since numerator + delta_plus >= denominator we already have | |
639 // 1 <= numerator/denominator < 10. Simply update the estimated_powe
r. | |
640 *decimal_point = estimated_power + 1; | |
641 } else { | |
642 *decimal_point = estimated_power; | |
643 numerator->Times10(); | |
644 if (Bignum::Equal(*delta_minus, *delta_plus)) { | |
645 delta_minus->Times10(); | |
646 delta_plus->AssignBignum(*delta_minus); | |
647 } else { | |
648 delta_minus->Times10(); | |
649 delta_plus->Times10(); | |
650 } | |
651 } | |
652 } | |
653 | 676 |
654 } // namespace double_conversion | 677 } // namespace double_conversion |
655 | 678 |
656 } // namespace WTF | 679 } // namespace WTF |
OLD | NEW |