1// Copyright 2012 the V8 project authors. All rights reserved.
2// Use of this source code is governed by a BSD-style license that can be
3// found in the LICENSE file.
4
5#include "src/strtod.h"
6
7#include <stdarg.h>
8#include <cmath>
9
10#include "src/bignum.h"
11#include "src/cached-powers.h"
12#include "src/double.h"
13#include "src/globals.h"
14#include "src/utils.h"
15
16namespace v8 {
17namespace internal {
18
19// 2^53 = 9007199254740992.
20// Any integer with at most 15 decimal digits will hence fit into a double
21// (which has a 53bit significand) without loss of precision.
22static const int kMaxExactDoubleIntegerDecimalDigits = 15;
23// 2^64 = 18446744073709551616 > 10^19
24static const int kMaxUint64DecimalDigits = 19;
25
26// Max double: 1.7976931348623157 x 10^308
27// Min non-zero double: 4.9406564584124654 x 10^-324
28// Any x >= 10^309 is interpreted as +infinity.
29// Any x <= 10^-324 is interpreted as 0.
30// Note that 2.5e-324 (despite being smaller than the min double) will be read
31// as non-zero (equal to the min non-zero double).
32static const int kMaxDecimalPower = 309;
33static const int kMinDecimalPower = -324;
34
35// 2^64 = 18446744073709551616
36static const uint64_t kMaxUint64 = V8_2PART_UINT64_C(0xFFFFFFFF, FFFFFFFF);
37
38
39static const double exact_powers_of_ten[] = {
40  1.0,  // 10^0
41  10.0,
42  100.0,
43  1000.0,
44  10000.0,
45  100000.0,
46  1000000.0,
47  10000000.0,
48  100000000.0,
49  1000000000.0,
50  10000000000.0,  // 10^10
51  100000000000.0,
52  1000000000000.0,
53  10000000000000.0,
54  100000000000000.0,
55  1000000000000000.0,
56  10000000000000000.0,
57  100000000000000000.0,
58  1000000000000000000.0,
59  10000000000000000000.0,
60  100000000000000000000.0,  // 10^20
61  1000000000000000000000.0,
62  // 10^22 = 0x21e19e0c9bab2400000 = 0x878678326eac9 * 2^22
63  10000000000000000000000.0
64};
65static const int kExactPowersOfTenSize = arraysize(exact_powers_of_ten);
66
67// Maximum number of significant digits in the decimal representation.
68// In fact the value is 772 (see conversions.cc), but to give us some margin
69// we round up to 780.
70static const int kMaxSignificantDecimalDigits = 780;
71
72static Vector<const char> TrimLeadingZeros(Vector<const char> buffer) {
73  for (int i = 0; i < buffer.length(); i++) {
74    if (buffer[i] != '0') {
75      return buffer.SubVector(i, buffer.length());
76    }
77  }
78  return Vector<const char>(buffer.start(), 0);
79}
80
81
82static Vector<const char> TrimTrailingZeros(Vector<const char> buffer) {
83  for (int i = buffer.length() - 1; i >= 0; --i) {
84    if (buffer[i] != '0') {
85      return buffer.SubVector(0, i + 1);
86    }
87  }
88  return Vector<const char>(buffer.start(), 0);
89}
90
91
92static void TrimToMaxSignificantDigits(Vector<const char> buffer,
93                                       int exponent,
94                                       char* significant_buffer,
95                                       int* significant_exponent) {
96  for (int i = 0; i < kMaxSignificantDecimalDigits - 1; ++i) {
97    significant_buffer[i] = buffer[i];
98  }
99  // The input buffer has been trimmed. Therefore the last digit must be
100  // different from '0'.
101  DCHECK(buffer[buffer.length() - 1] != '0');
102  // Set the last digit to be non-zero. This is sufficient to guarantee
103  // correct rounding.
104  significant_buffer[kMaxSignificantDecimalDigits - 1] = '1';
105  *significant_exponent =
106      exponent + (buffer.length() - kMaxSignificantDecimalDigits);
107}
108
109
110// Reads digits from the buffer and converts them to a uint64.
111// Reads in as many digits as fit into a uint64.
112// When the string starts with "1844674407370955161" no further digit is read.
113// Since 2^64 = 18446744073709551616 it would still be possible read another
114// digit if it was less or equal than 6, but this would complicate the code.
115static uint64_t ReadUint64(Vector<const char> buffer,
116                           int* number_of_read_digits) {
117  uint64_t result = 0;
118  int i = 0;
119  while (i < buffer.length() && result <= (kMaxUint64 / 10 - 1)) {
120    int digit = buffer[i++] - '0';
121    DCHECK(0 <= digit && digit <= 9);
122    result = 10 * result + digit;
123  }
124  *number_of_read_digits = i;
125  return result;
126}
127
128
129// Reads a DiyFp from the buffer.
130// The returned DiyFp is not necessarily normalized.
131// If remaining_decimals is zero then the returned DiyFp is accurate.
132// Otherwise it has been rounded and has error of at most 1/2 ulp.
133static void ReadDiyFp(Vector<const char> buffer,
134                      DiyFp* result,
135                      int* remaining_decimals) {
136  int read_digits;
137  uint64_t significand = ReadUint64(buffer, &read_digits);
138  if (buffer.length() == read_digits) {
139    *result = DiyFp(significand, 0);
140    *remaining_decimals = 0;
141  } else {
142    // Round the significand.
143    if (buffer[read_digits] >= '5') {
144      significand++;
145    }
146    // Compute the binary exponent.
147    int exponent = 0;
148    *result = DiyFp(significand, exponent);
149    *remaining_decimals = buffer.length() - read_digits;
150  }
151}
152
153
154static bool DoubleStrtod(Vector<const char> trimmed,
155                         int exponent,
156                         double* result) {
157#if (V8_TARGET_ARCH_IA32 || V8_TARGET_ARCH_X87 || defined(USE_SIMULATOR)) && \
158    !defined(_MSC_VER)
159  // On x86 the floating-point stack can be 64 or 80 bits wide. If it is
160  // 80 bits wide (as is the case on Linux) then double-rounding occurs and the
161  // result is not accurate.
162  // We know that Windows32 with MSVC, unlike with MinGW32, uses 64 bits and is
163  // therefore accurate.
164  // Note that the ARM and MIPS simulators are compiled for 32bits. They
165  // therefore exhibit the same problem.
166  return false;
167#endif
168  if (trimmed.length() <= kMaxExactDoubleIntegerDecimalDigits) {
169    int read_digits;
170    // The trimmed input fits into a double.
171    // If the 10^exponent (resp. 10^-exponent) fits into a double too then we
172    // can compute the result-double simply by multiplying (resp. dividing) the
173    // two numbers.
174    // This is possible because IEEE guarantees that floating-point operations
175    // return the best possible approximation.
176    if (exponent < 0 && -exponent < kExactPowersOfTenSize) {
177      // 10^-exponent fits into a double.
178      *result = static_cast<double>(ReadUint64(trimmed, &read_digits));
179      DCHECK(read_digits == trimmed.length());
180      *result /= exact_powers_of_ten[-exponent];
181      return true;
182    }
183    if (0 <= exponent && exponent < kExactPowersOfTenSize) {
184      // 10^exponent fits into a double.
185      *result = static_cast<double>(ReadUint64(trimmed, &read_digits));
186      DCHECK(read_digits == trimmed.length());
187      *result *= exact_powers_of_ten[exponent];
188      return true;
189    }
190    int remaining_digits =
191        kMaxExactDoubleIntegerDecimalDigits - trimmed.length();
192    if ((0 <= exponent) &&
193        (exponent - remaining_digits < kExactPowersOfTenSize)) {
194      // The trimmed string was short and we can multiply it with
195      // 10^remaining_digits. As a result the remaining exponent now fits
196      // into a double too.
197      *result = static_cast<double>(ReadUint64(trimmed, &read_digits));
198      DCHECK(read_digits == trimmed.length());
199      *result *= exact_powers_of_ten[remaining_digits];
200      *result *= exact_powers_of_ten[exponent - remaining_digits];
201      return true;
202    }
203  }
204  return false;
205}
206
207
208// Returns 10^exponent as an exact DiyFp.
209// The given exponent must be in the range [1; kDecimalExponentDistance[.
210static DiyFp AdjustmentPowerOfTen(int exponent) {
211  DCHECK(0 < exponent);
212  DCHECK(exponent < PowersOfTenCache::kDecimalExponentDistance);
213  // Simply hardcode the remaining powers for the given decimal exponent
214  // distance.
215  DCHECK(PowersOfTenCache::kDecimalExponentDistance == 8);
216  switch (exponent) {
217    case 1: return DiyFp(V8_2PART_UINT64_C(0xa0000000, 00000000), -60);
218    case 2: return DiyFp(V8_2PART_UINT64_C(0xc8000000, 00000000), -57);
219    case 3: return DiyFp(V8_2PART_UINT64_C(0xfa000000, 00000000), -54);
220    case 4: return DiyFp(V8_2PART_UINT64_C(0x9c400000, 00000000), -50);
221    case 5: return DiyFp(V8_2PART_UINT64_C(0xc3500000, 00000000), -47);
222    case 6: return DiyFp(V8_2PART_UINT64_C(0xf4240000, 00000000), -44);
223    case 7: return DiyFp(V8_2PART_UINT64_C(0x98968000, 00000000), -40);
224    default:
225      UNREACHABLE();
226      return DiyFp(0, 0);
227  }
228}
229
230
231// If the function returns true then the result is the correct double.
232// Otherwise it is either the correct double or the double that is just below
233// the correct double.
234static bool DiyFpStrtod(Vector<const char> buffer,
235                        int exponent,
236                        double* result) {
237  DiyFp input;
238  int remaining_decimals;
239  ReadDiyFp(buffer, &input, &remaining_decimals);
240  // Since we may have dropped some digits the input is not accurate.
241  // If remaining_decimals is different than 0 than the error is at most
242  // .5 ulp (unit in the last place).
243  // We don't want to deal with fractions and therefore keep a common
244  // denominator.
245  const int kDenominatorLog = 3;
246  const int kDenominator = 1 << kDenominatorLog;
247  // Move the remaining decimals into the exponent.
248  exponent += remaining_decimals;
249  int64_t error = (remaining_decimals == 0 ? 0 : kDenominator / 2);
250
251  int old_e = input.e();
252  input.Normalize();
253  error <<= old_e - input.e();
254
255  DCHECK(exponent <= PowersOfTenCache::kMaxDecimalExponent);
256  if (exponent < PowersOfTenCache::kMinDecimalExponent) {
257    *result = 0.0;
258    return true;
259  }
260  DiyFp cached_power;
261  int cached_decimal_exponent;
262  PowersOfTenCache::GetCachedPowerForDecimalExponent(exponent,
263                                                     &cached_power,
264                                                     &cached_decimal_exponent);
265
266  if (cached_decimal_exponent != exponent) {
267    int adjustment_exponent = exponent - cached_decimal_exponent;
268    DiyFp adjustment_power = AdjustmentPowerOfTen(adjustment_exponent);
269    input.Multiply(adjustment_power);
270    if (kMaxUint64DecimalDigits - buffer.length() >= adjustment_exponent) {
271      // The product of input with the adjustment power fits into a 64 bit
272      // integer.
273      DCHECK(DiyFp::kSignificandSize == 64);
274    } else {
275      // The adjustment power is exact. There is hence only an error of 0.5.
276      error += kDenominator / 2;
277    }
278  }
279
280  input.Multiply(cached_power);
281  // The error introduced by a multiplication of a*b equals
282  //   error_a + error_b + error_a*error_b/2^64 + 0.5
283  // Substituting a with 'input' and b with 'cached_power' we have
284  //   error_b = 0.5  (all cached powers have an error of less than 0.5 ulp),
285  //   error_ab = 0 or 1 / kDenominator > error_a*error_b/ 2^64
286  int error_b = kDenominator / 2;
287  int error_ab = (error == 0 ? 0 : 1);  // We round up to 1.
288  int fixed_error = kDenominator / 2;
289  error += error_b + error_ab + fixed_error;
290
291  old_e = input.e();
292  input.Normalize();
293  error <<= old_e - input.e();
294
295  // See if the double's significand changes if we add/subtract the error.
296  int order_of_magnitude = DiyFp::kSignificandSize + input.e();
297  int effective_significand_size =
298      Double::SignificandSizeForOrderOfMagnitude(order_of_magnitude);
299  int precision_digits_count =
300      DiyFp::kSignificandSize - effective_significand_size;
301  if (precision_digits_count + kDenominatorLog >= DiyFp::kSignificandSize) {
302    // This can only happen for very small denormals. In this case the
303    // half-way multiplied by the denominator exceeds the range of an uint64.
304    // Simply shift everything to the right.
305    int shift_amount = (precision_digits_count + kDenominatorLog) -
306        DiyFp::kSignificandSize + 1;
307    input.set_f(input.f() >> shift_amount);
308    input.set_e(input.e() + shift_amount);
309    // We add 1 for the lost precision of error, and kDenominator for
310    // the lost precision of input.f().
311    error = (error >> shift_amount) + 1 + kDenominator;
312    precision_digits_count -= shift_amount;
313  }
314  // We use uint64_ts now. This only works if the DiyFp uses uint64_ts too.
315  DCHECK(DiyFp::kSignificandSize == 64);
316  DCHECK(precision_digits_count < 64);
317  uint64_t one64 = 1;
318  uint64_t precision_bits_mask = (one64 << precision_digits_count) - 1;
319  uint64_t precision_bits = input.f() & precision_bits_mask;
320  uint64_t half_way = one64 << (precision_digits_count - 1);
321  precision_bits *= kDenominator;
322  half_way *= kDenominator;
323  DiyFp rounded_input(input.f() >> precision_digits_count,
324                      input.e() + precision_digits_count);
325  if (precision_bits >= half_way + error) {
326    rounded_input.set_f(rounded_input.f() + 1);
327  }
328  // If the last_bits are too close to the half-way case than we are too
329  // inaccurate and round down. In this case we return false so that we can
330  // fall back to a more precise algorithm.
331
332  *result = Double(rounded_input).value();
333  if (half_way - error < precision_bits && precision_bits < half_way + error) {
334    // Too imprecise. The caller will have to fall back to a slower version.
335    // However the returned number is guaranteed to be either the correct
336    // double, or the next-lower double.
337    return false;
338  } else {
339    return true;
340  }
341}
342
343
344// Returns the correct double for the buffer*10^exponent.
345// The variable guess should be a close guess that is either the correct double
346// or its lower neighbor (the nearest double less than the correct one).
347// Preconditions:
348//   buffer.length() + exponent <= kMaxDecimalPower + 1
349//   buffer.length() + exponent > kMinDecimalPower
350//   buffer.length() <= kMaxDecimalSignificantDigits
351static double BignumStrtod(Vector<const char> buffer,
352                           int exponent,
353                           double guess) {
354  if (guess == V8_INFINITY) {
355    return guess;
356  }
357
358  DiyFp upper_boundary = Double(guess).UpperBoundary();
359
360  DCHECK(buffer.length() + exponent <= kMaxDecimalPower + 1);
361  DCHECK(buffer.length() + exponent > kMinDecimalPower);
362  DCHECK(buffer.length() <= kMaxSignificantDecimalDigits);
363  // Make sure that the Bignum will be able to hold all our numbers.
364  // Our Bignum implementation has a separate field for exponents. Shifts will
365  // consume at most one bigit (< 64 bits).
366  // ln(10) == 3.3219...
367  DCHECK(((kMaxDecimalPower + 1) * 333 / 100) < Bignum::kMaxSignificantBits);
368  Bignum input;
369  Bignum boundary;
370  input.AssignDecimalString(buffer);
371  boundary.AssignUInt64(upper_boundary.f());
372  if (exponent >= 0) {
373    input.MultiplyByPowerOfTen(exponent);
374  } else {
375    boundary.MultiplyByPowerOfTen(-exponent);
376  }
377  if (upper_boundary.e() > 0) {
378    boundary.ShiftLeft(upper_boundary.e());
379  } else {
380    input.ShiftLeft(-upper_boundary.e());
381  }
382  int comparison = Bignum::Compare(input, boundary);
383  if (comparison < 0) {
384    return guess;
385  } else if (comparison > 0) {
386    return Double(guess).NextDouble();
387  } else if ((Double(guess).Significand() & 1) == 0) {
388    // Round towards even.
389    return guess;
390  } else {
391    return Double(guess).NextDouble();
392  }
393}
394
395
396double Strtod(Vector<const char> buffer, int exponent) {
397  Vector<const char> left_trimmed = TrimLeadingZeros(buffer);
398  Vector<const char> trimmed = TrimTrailingZeros(left_trimmed);
399  exponent += left_trimmed.length() - trimmed.length();
400  if (trimmed.length() == 0) return 0.0;
401  if (trimmed.length() > kMaxSignificantDecimalDigits) {
402    char significant_buffer[kMaxSignificantDecimalDigits];
403    int significant_exponent;
404    TrimToMaxSignificantDigits(trimmed, exponent,
405                               significant_buffer, &significant_exponent);
406    return Strtod(Vector<const char>(significant_buffer,
407                                     kMaxSignificantDecimalDigits),
408                  significant_exponent);
409  }
410  if (exponent + trimmed.length() - 1 >= kMaxDecimalPower) return V8_INFINITY;
411  if (exponent + trimmed.length() <= kMinDecimalPower) return 0.0;
412
413  double guess;
414  if (DoubleStrtod(trimmed, exponent, &guess) ||
415      DiyFpStrtod(trimmed, exponent, &guess)) {
416    return guess;
417  }
418  return BignumStrtod(trimmed, exponent, guess);
419}
420
421}  // namespace internal
422}  // namespace v8
423