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