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