OLD | NEW |
| (Empty) |
1 // Copyright 2010 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 "v8.h" | |
29 | |
30 #include "grisu3.h" | |
31 | |
32 #include "cached_powers.h" | |
33 #include "diy_fp.h" | |
34 #include "double.h" | |
35 | |
36 namespace v8 { | |
37 namespace internal { | |
38 | |
39 template <int alpha = -60, int gamma = -32> | |
40 class Grisu3 { | |
41 public: | |
42 // Provides a decimal representation of v. | |
43 // Returns true if it succeeds, otherwise the result can not be trusted. | |
44 // There will be *length digits inside the buffer (not null-terminated). | |
45 // If the function returns true then | |
46 // v == (double) (buffer * 10^decimal_exponent). | |
47 // The digits in the buffer are the shortest representation possible: no | |
48 // 0.099999999999 instead of 0.1. | |
49 // The last digit will be closest to the actual v. That is, even if several | |
50 // digits might correctly yield 'v' when read again, the closest will be | |
51 // computed. | |
52 static bool grisu3(double v, | |
53 char* buffer, int* length, int* decimal_exponent); | |
54 | |
55 private: | |
56 // Rounds the buffer according to the rest. | |
57 // If there is too much imprecision to round then false is returned. | |
58 // Similarily false is returned when the buffer is not within Delta. | |
59 static bool RoundWeed(char* buffer, int len, uint64_t wp_W, uint64_t Delta, | |
60 uint64_t rest, uint64_t ten_kappa, uint64_t ulp); | |
61 // Dispatches to the a specialized digit-generation routine. The chosen | |
62 // routine depends on w.e (which in turn depends on alpha and gamma). | |
63 // Currently there is only one digit-generation routine, but it would be easy | |
64 // to add others. | |
65 static bool DigitGen(DiyFp low, DiyFp w, DiyFp high, | |
66 char* buffer, int* len, int* kappa); | |
67 // Generates w's digits. The result is the shortest in the interval low-high. | |
68 // All DiyFp are assumed to be imprecise and this function takes this | |
69 // imprecision into account. If the function cannot compute the best | |
70 // representation (due to the imprecision) then false is returned. | |
71 static bool DigitGen_m60_m32(DiyFp low, DiyFp w, DiyFp high, | |
72 char* buffer, int* length, int* kappa); | |
73 }; | |
74 | |
75 | |
76 template<int alpha, int gamma> | |
77 bool Grisu3<alpha, gamma>::grisu3( | |
78 double v, char* buffer, int* length, int* decimal_exponent) { | |
79 DiyFp w = Double(v).AsNormalizedDiyFp(); | |
80 // boundary_minus and boundary_plus are the boundaries between v and its | |
81 // neighbors. Any number strictly between boundary_minus and boundary_plus | |
82 // will round to v when read as double. | |
83 // Grisu3 will never output representations that lie exactly on a boundary. | |
84 DiyFp boundary_minus, boundary_plus; | |
85 Double(v).NormalizedBoundaries(&boundary_minus, &boundary_plus); | |
86 ASSERT(boundary_plus.e() == w.e()); | |
87 DiyFp ten_mk; // Cached power of ten: 10^-k | |
88 int mk; // -k | |
89 GetCachedPower(w.e() + DiyFp::kSignificandSize, alpha, gamma, &mk, &ten_mk); | |
90 ASSERT(alpha <= w.e() + ten_mk.e() + DiyFp::kSignificandSize && | |
91 gamma >= w.e() + ten_mk.e() + DiyFp::kSignificandSize); | |
92 // Note that ten_mk is only an approximation of 10^-k. A DiyFp only contains a | |
93 // 64 bit significand and ten_mk is thus only precise up to 64 bits. | |
94 | |
95 // The DiyFp::Times procedure rounds its result, and ten_mk is approximated | |
96 // too. The variable scaled_w (as well as scaled_boundary_minus/plus) are now | |
97 // off by a small amount. | |
98 // In fact: scaled_w - w*10^k < 1ulp (unit in the last place) of scaled_w. | |
99 // In other words: let f = scaled_w.f() and e = scaled_w.e(), then | |
100 // (f-1) * 2^e < w*10^k < (f+1) * 2^e | |
101 DiyFp scaled_w = DiyFp::Times(w, ten_mk); | |
102 ASSERT(scaled_w.e() == | |
103 boundary_plus.e() + ten_mk.e() + DiyFp::kSignificandSize); | |
104 // In theory it would be possible to avoid some recomputations by computing | |
105 // the difference between w and boundary_minus/plus (a power of 2) and to | |
106 // compute scaled_boundary_minus/plus by subtracting/adding from | |
107 // scaled_w. However the code becomes much less readable and the speed | |
108 // enhancements are not terriffic. | |
109 DiyFp scaled_boundary_minus = DiyFp::Times(boundary_minus, ten_mk); | |
110 DiyFp scaled_boundary_plus = DiyFp::Times(boundary_plus, ten_mk); | |
111 | |
112 // DigitGen will generate the digits of scaled_w. Therefore we have | |
113 // v == (double) (scaled_w * 10^-mk). | |
114 // Set decimal_exponent == -mk and pass it to DigitGen. If scaled_w is not an | |
115 // integer than it will be updated. For instance if scaled_w == 1.23 then | |
116 // the buffer will be filled with "123" und the decimal_exponent will be | |
117 // decreased by 2. | |
118 int kappa; | |
119 bool result = DigitGen(scaled_boundary_minus, scaled_w, scaled_boundary_plus, | |
120 buffer, length, &kappa); | |
121 *decimal_exponent = -mk + kappa; | |
122 return result; | |
123 } | |
124 | |
125 // Generates the digits of input number w. | |
126 // w is a floating-point number (DiyFp), consisting of a significand and an | |
127 // exponent. Its exponent is bounded by alpha and gamma. Typically alpha >= -63 | |
128 // and gamma <= 3. | |
129 // Returns false if it fails, in which case the generated digits in the buffer | |
130 // should not be used. | |
131 // Preconditions: | |
132 // * low, w and high are correct up to 1 ulp (unit in the last place). That | |
133 // is, their error must be less that a unit of their last digits. | |
134 // * low.e() == w.e() == high.e() | |
135 // * low < w < high, and taking into account their error: low~ <= high~ | |
136 // * alpha <= w.e() <= gamma | |
137 // Postconditions: returns false if procedure fails. | |
138 // otherwise: | |
139 // * buffer is not null-terminated, but len contains the number of digits. | |
140 // * buffer contains the shortest possible decimal digit-sequence | |
141 // such that LOW < buffer * 10^kappa < HIGH, where LOW and HIGH are the | |
142 // correct values of low and high (without their error). | |
143 // * if more than one decimal representation gives the minimal number of | |
144 // decimal digits then the one closest to W (where W is the correct value | |
145 // of w) is chosen. | |
146 // Remark: this procedure takes into account the imprecision of its input | |
147 // numbers. If the precision is not enough to guarantee all the postconditions | |
148 // then false is returned. This usually happens rarely (~0.5%). | |
149 template<int alpha, int gamma> | |
150 bool Grisu3<alpha, gamma>::DigitGen( | |
151 DiyFp low, DiyFp w, DiyFp high, char* buffer, int* len, int* kappa) { | |
152 ASSERT(low.e() == w.e() && w.e() == high.e()); | |
153 ASSERT(low.f() + 1 <= high.f() - 1); | |
154 ASSERT(alpha <= w.e() && w.e() <= gamma); | |
155 // The following tests use alpha and gamma to avoid unnecessary dynamic tests. | |
156 if ((alpha >= -60 && gamma <= -32) || // -60 <= w.e() <= -32 | |
157 (alpha <= -32 && gamma >= -60 && // Alpha/gamma overlaps -60/-32 region. | |
158 -60 <= w.e() && w.e() <= -32)) { | |
159 return DigitGen_m60_m32(low, w, high, buffer, len, kappa); | |
160 } else { | |
161 // A simple adaption of the special case -60/-32 would allow greater ranges | |
162 // of alpha/gamma and thus reduce the number of precomputed cached powers of | |
163 // ten. | |
164 UNIMPLEMENTED(); | |
165 return false; | |
166 } | |
167 } | |
168 | |
169 static const uint32_t kTen4 = 10000; | |
170 static const uint32_t kTen5 = 100000; | |
171 static const uint32_t kTen6 = 1000000; | |
172 static const uint32_t kTen7 = 10000000; | |
173 static const uint32_t kTen8 = 100000000; | |
174 static const uint32_t kTen9 = 1000000000; | |
175 | |
176 // Returns the biggest power of ten that is <= than the given number. We | |
177 // furthermore receive the maximum number of bits 'number' has. | |
178 // If number_bits == 0 then 0^-1 is returned | |
179 // The number of bits must be <= 32. | |
180 static void BiggestPowerTen(uint32_t number, int number_bits, | |
181 uint32_t* power, int* exponent) { | |
182 switch (number_bits) { | |
183 case 32: | |
184 case 31: | |
185 case 30: | |
186 if (kTen9 <= number) { | |
187 *power = kTen9; | |
188 *exponent = 9; | |
189 break; | |
190 } // else fallthrough | |
191 case 29: | |
192 case 28: | |
193 case 27: | |
194 if (kTen8 <= number) { | |
195 *power = kTen8; | |
196 *exponent = 8; | |
197 break; | |
198 } // else fallthrough | |
199 case 26: | |
200 case 25: | |
201 case 24: | |
202 if (kTen7 <= number) { | |
203 *power = kTen7; | |
204 *exponent = 7; | |
205 break; | |
206 } // else fallthrough | |
207 case 23: | |
208 case 22: | |
209 case 21: | |
210 case 20: | |
211 if (kTen6 <= number) { | |
212 *power = kTen6; | |
213 *exponent = 6; | |
214 break; | |
215 } // else fallthrough | |
216 case 19: | |
217 case 18: | |
218 case 17: | |
219 if (kTen5 <= number) { | |
220 *power = kTen5; | |
221 *exponent = 5; | |
222 break; | |
223 } // else fallthrough | |
224 case 16: | |
225 case 15: | |
226 case 14: | |
227 if (kTen4 <= number) { | |
228 *power = kTen4; | |
229 *exponent = 4; | |
230 break; | |
231 } // else fallthrough | |
232 case 13: | |
233 case 12: | |
234 case 11: | |
235 case 10: | |
236 if (1000 <= number) { | |
237 *power = 1000; | |
238 *exponent = 3; | |
239 break; | |
240 } // else fallthrough | |
241 case 9: | |
242 case 8: | |
243 case 7: | |
244 if (100 <= number) { | |
245 *power = 100; | |
246 *exponent = 2; | |
247 break; | |
248 } // else fallthrough | |
249 case 6: | |
250 case 5: | |
251 case 4: | |
252 if (10 <= number) { | |
253 *power = 10; | |
254 *exponent = 1; | |
255 break; | |
256 } // else fallthrough | |
257 case 3: | |
258 case 2: | |
259 case 1: | |
260 if (1 <= number) { | |
261 *power = 1; | |
262 *exponent = 0; | |
263 break; | |
264 } // else fallthrough | |
265 case 0: | |
266 *power = 0; | |
267 *exponent = -1; | |
268 break; | |
269 default: | |
270 // Following assignments are here to silence compiler warnings. | |
271 *power = 0; | |
272 *exponent = 0; | |
273 UNREACHABLE(); | |
274 } | |
275 } | |
276 | |
277 | |
278 // Same comments as for DigitGen but with additional precondition: | |
279 // -60 <= w.e() <= -32 | |
280 // | |
281 // Say, for the sake of example, that | |
282 // w.e() == -48, and w.f() == 0x1234567890abcdef | |
283 // w's value can be computed by w.f() * 2^w.e() | |
284 // We can obtain w's integral digits by simply shifting w.f() by -w.e(). | |
285 // -> w's integral part is 0x1234 | |
286 // w's fractional part is therefore 0x567890abcdef. | |
287 // Printing w's integral part is easy (simply print 0x1234 in decimal). | |
288 // In order to print its fraction we repeatedly multiply the fraction by 10 and | |
289 // get each digit. Example the first digit after the comma would be computed by | |
290 // (0x567890abcdef * 10) >> 48. -> 3 | |
291 // The whole thing becomes slightly more complicated because we want to stop | |
292 // once we have enough digits. That is, once the digits inside the buffer | |
293 // represent 'w' we can stop. Everything inside the interval low - high | |
294 // represents w. However we have to pay attention to low, high and w's | |
295 // imprecision. | |
296 template<int alpha, int gamma> | |
297 bool Grisu3<alpha, gamma>::DigitGen_m60_m32( | |
298 DiyFp low, DiyFp w, DiyFp high, char* buffer, int* length, int* kappa) { | |
299 // low, w and high are imprecise, but by less than one ulp (unit in the last | |
300 // place). | |
301 // If we remove (resp. add) 1 ulp from low (resp. high) we are certain that | |
302 // the new numbers are outside of the interval we want the final | |
303 // representation to lie in. | |
304 // Inversely adding (resp. removing) 1 ulp from low (resp. high) would yield | |
305 // numbers that are certain to lie in the interval. We will use this fact | |
306 // later on. | |
307 // We will now start by generating the digits within the uncertain | |
308 // interval. Later we will weed out representations that lie outside the safe | |
309 // interval and thus _might_ lie outside the correct interval. | |
310 uint64_t unit = 1; | |
311 DiyFp too_low = DiyFp(low.f() - unit, low.e()); | |
312 DiyFp too_high = DiyFp(high.f() + unit, high.e()); | |
313 // too_low and too_high are guaranteed to lie outside the interval we want the | |
314 // generated number in. | |
315 DiyFp unsafe_interval = DiyFp::Minus(too_high, too_low); | |
316 // We now cut the input number into two parts: the integral digits and the | |
317 // fractionals. We will not write any decimal separator though, but adapt | |
318 // kappa instead. | |
319 // Reminder: we are currently computing the digits (stored inside the buffer) | |
320 // such that: too_low < buffer * 10^kappa < too_high | |
321 // We use too_high for the digit_generation and stop as soon as possible. | |
322 // If we stop early we effectively round down. | |
323 DiyFp one = DiyFp(static_cast<uint64_t>(1) << -w.e(), w.e()); | |
324 uint32_t integrals = too_high.f() >> -one.e(); // Division by one. | |
325 uint64_t fractionals = too_high.f() & (one.f() - 1); // Modulo by one. | |
326 uint32_t divider; | |
327 int divider_exponent; | |
328 BiggestPowerTen(integrals, DiyFp::kSignificandSize - (-one.e()), | |
329 ÷r, ÷r_exponent); | |
330 *kappa = divider_exponent + 1; | |
331 *length = 0; | |
332 // Loop invariant: buffer = too_high / 10^kappa (integer division) | |
333 // The invariant holds for the first iteration: kappa has been initialized | |
334 // with the divider exponent + 1. And the divider is the biggest power of ten | |
335 // that fits into the bits that had been reserved for the integrals. | |
336 while (*kappa > 0) { | |
337 int digit = integrals / divider; | |
338 buffer[*length] = '0' + digit; | |
339 (*length)++; | |
340 integrals %= divider; | |
341 (*kappa)--; | |
342 // Note that kappa now equals the exponent of the divider and that the | |
343 // invariant thus holds again. | |
344 uint64_t rest = | |
345 (static_cast<uint64_t>(integrals) << -one.e()) + fractionals; | |
346 // Invariant: too_high = buffer * 10^kappa + DiyFp(rest, one.e()) | |
347 // Reminder: unsafe_interval.e() == one.e() | |
348 if (rest < unsafe_interval.f()) { | |
349 // Rounding down (by not emitting the remaining digits) yields a number | |
350 // that lies within the unsafe interval. | |
351 return RoundWeed(buffer, *length, DiyFp::Minus(too_high, w).f(), | |
352 unsafe_interval.f(), rest, | |
353 static_cast<uint64_t>(divider) << -one.e(), unit); | |
354 } | |
355 divider /= 10; | |
356 } | |
357 | |
358 // The integrals have been generated. We are at the point of the decimal | |
359 // separator. In the following loop we simply multiply the remaining digits by | |
360 // 10 and divide by one. We just need to pay attention to multiply associated | |
361 // data (like the interval or 'unit'), too. | |
362 // Instead of multiplying by 10 we multiply by 5 (cheaper operation) and | |
363 // increase its (imaginary) exponent. At the same time we decrease the | |
364 // divider's (one's) exponent and shift its significand. | |
365 // Basically, if fractionals was a DiyFp (with fractionals.e == one.e): | |
366 // fractionals.f *= 10; | |
367 // fractionals.f >>= 1; fractionals.e++; // value remains unchanged. | |
368 // one.f >>= 1; one.e++; // value remains unchanged. | |
369 // and we have again fractionals.e == one.e which allows us to divide | |
370 // fractionals.f() by one.f() | |
371 // We simply combine the *= 10 and the >>= 1. | |
372 while (true) { | |
373 fractionals *= 5; | |
374 unit *= 5; | |
375 unsafe_interval.set_f(unsafe_interval.f() * 5); | |
376 unsafe_interval.set_e(unsafe_interval.e() + 1); // Will be optimized out. | |
377 one.set_f(one.f() >> 1); | |
378 one.set_e(one.e() + 1); | |
379 int digit = fractionals >> -one.e(); // Integer division by one. | |
380 buffer[*length] = '0' + digit; | |
381 (*length)++; | |
382 fractionals &= one.f() - 1; // Modulo by one. | |
383 (*kappa)--; | |
384 if (fractionals < unsafe_interval.f()) { | |
385 return RoundWeed(buffer, *length, DiyFp::Minus(too_high, w).f() * unit, | |
386 unsafe_interval.f(), fractionals, one.f(), unit); | |
387 } | |
388 } | |
389 } | |
390 | |
391 | |
392 // Rounds the given generated digits in the buffer and weeds out generated | |
393 // digits that are not in the safe interval, or where we cannot find a rounded | |
394 // representation. | |
395 // Input: * buffer containing the digits of too_high / 10^kappa | |
396 // * the buffer's length | |
397 // * distance_too_high_w == (too_high - w).f() * unit | |
398 // * unsafe_interval == (too_high - too_low).f() * unit | |
399 // * rest = (too_high - buffer * 10^kappa).f() * unit | |
400 // * ten_kappa = 10^kappa * unit | |
401 // * unit = the common multiplier | |
402 // Output: returns true on success. | |
403 // Modifies the generated digits in the buffer to approach (round towards) w. | |
404 template<int alpha, int gamma> | |
405 bool Grisu3<alpha, gamma>::RoundWeed( | |
406 char* buffer, int length, uint64_t distance_too_high_w, | |
407 uint64_t unsafe_interval, uint64_t rest, uint64_t ten_kappa, | |
408 uint64_t unit) { | |
409 uint64_t small_distance = distance_too_high_w - unit; | |
410 uint64_t big_distance = distance_too_high_w + unit; | |
411 // Let w- = too_high - big_distance, and | |
412 // w+ = too_high - small_distance. | |
413 // Note: w- < w < w+ | |
414 // | |
415 // The real w (* unit) must lie somewhere inside the interval | |
416 // ]w-; w+[ (often written as "(w-; w+)") | |
417 | |
418 // Basically the buffer currently contains a number in the unsafe interval | |
419 // ]too_low; too_high[ with too_low < w < too_high | |
420 // | |
421 // By generating the digits of too_high we got the biggest last digit. | |
422 // In the case that w+ < buffer < too_high we try to decrement the buffer. | |
423 // This way the buffer approaches (rounds towards) w. | |
424 // There are 3 conditions that stop the decrementation process: | |
425 // 1) the buffer is already below w+ | |
426 // 2) decrementing the buffer would make it leave the unsafe interval | |
427 // 3) decrementing the buffer would yield a number below w+ and farther away | |
428 // than the current number. In other words: | |
429 // (buffer{-1} < w+) && w+ - buffer{-1} > buffer - w+ | |
430 // Instead of using the buffer directly we use its distance to too_high. | |
431 // Conceptually rest ~= too_high - buffer | |
432 while (rest < small_distance && // Negated condition 1 | |
433 unsafe_interval - rest >= ten_kappa && // Negated condition 2 | |
434 (rest + ten_kappa < small_distance || // buffer{-1} > w+ | |
435 small_distance - rest >= rest + ten_kappa - small_distance)) { | |
436 buffer[length - 1]--; | |
437 rest += ten_kappa; | |
438 } | |
439 | |
440 // We have approached w+ as much as possible. We now test if approaching w- | |
441 // would require changing the buffer. If yes, then we have two possible | |
442 // representations close to w, but we cannot decide which one is closer. | |
443 if (rest < big_distance && | |
444 unsafe_interval - rest >= ten_kappa && | |
445 (rest + ten_kappa < big_distance || | |
446 big_distance - rest > rest + ten_kappa - big_distance)) { | |
447 return false; | |
448 } | |
449 | |
450 // Weeding test. | |
451 // The safe interval is [too_low + 2 ulp; too_high - 2 ulp] | |
452 // Since too_low = too_high - unsafe_interval this is equivalent too | |
453 // [too_high - unsafe_interval + 4 ulp; too_high - 2 ulp] | |
454 // Conceptually we have: rest ~= too_high - buffer | |
455 return (2 * unit <= rest) && (rest <= unsafe_interval - 4 * unit); | |
456 } | |
457 | |
458 | |
459 bool grisu3(double v, | |
460 char* buffer, int* sign, int* length, int* decimal_point) { | |
461 ASSERT(v != 0); | |
462 ASSERT(!Double(v).IsSpecial()); | |
463 | |
464 if (v < 0) { | |
465 v = -v; | |
466 *sign = 1; | |
467 } else { | |
468 *sign = 0; | |
469 } | |
470 int decimal_exponent; | |
471 bool result = Grisu3<-60, -32>::grisu3(v, buffer, length, &decimal_exponent); | |
472 *decimal_point = *length + decimal_exponent; | |
473 buffer[*length] = '\0'; | |
474 return result; | |
475 } | |
476 | |
477 } } // namespace v8::internal | |
OLD | NEW |