blob: 33fc5462377a134507b8aa838a4cffda2a5b0043 [file] [log] [blame]
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001// Copyright 2020 The Wuffs Authors.
2//
3// Licensed under the Apache License, Version 2.0 (the "License");
4// you may not use this file except in compliance with the License.
5// You may obtain a copy of the License at
6//
7// https://www.apache.org/licenses/LICENSE-2.0
8//
9// Unless required by applicable law or agreed to in writing, software
10// distributed under the License is distributed on an "AS IS" BASIS,
11// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12// See the License for the specific language governing permissions and
13// limitations under the License.
14
15// ---------------- IEEE 754 Floating Point
16
Nigel Tao7bf7cf22020-07-12 16:23:15 +100017WUFFS_BASE__MAYBE_STATIC wuffs_base__lossy_value_u16 //
Nigel Taoa3931d52020-07-12 21:06:44 +100018wuffs_base__ieee_754_bit_representation__from_f64_to_u16_truncate(double f) {
Nigel Tao7bf7cf22020-07-12 16:23:15 +100019 uint64_t u = 0;
20 if (sizeof(uint64_t) == sizeof(double)) {
21 memcpy(&u, &f, sizeof(uint64_t));
22 }
Nigel Tao56d90962020-07-12 21:11:49 +100023 uint16_t neg = ((uint16_t)((u >> 63) << 15));
Nigel Tao7bf7cf22020-07-12 16:23:15 +100024 u &= 0x7FFFFFFFFFFFFFFF;
25 uint64_t exp = u >> 52;
26 uint64_t man = u & 0x000FFFFFFFFFFFFF;
27
28 if (exp == 0x7FF) {
29 if (man == 0) { // Infinity.
30 wuffs_base__lossy_value_u16 ret;
31 ret.value = neg | 0x7C00;
32 ret.lossy = false;
33 return ret;
34 }
35 // NaN. Shift the 52 mantissa bits to 10 mantissa bits, keeping the most
36 // significant mantissa bit (quiet vs signaling NaNs). Also set the low 9
37 // bits of ret.value so that the 10-bit mantissa is non-zero.
38 wuffs_base__lossy_value_u16 ret;
39 ret.value = neg | 0x7DFF | ((uint16_t)(man >> 42));
40 ret.lossy = false;
41 return ret;
42
43 } else if (exp > 0x40E) { // Truncate to the largest finite f16.
44 wuffs_base__lossy_value_u16 ret;
45 ret.value = neg | 0x7BFF;
46 ret.lossy = true;
47 return ret;
48
49 } else if (exp <= 0x3E6) { // Truncate to zero.
50 wuffs_base__lossy_value_u16 ret;
51 ret.value = neg;
52 ret.lossy = (u != 0);
53 return ret;
54
55 } else if (exp <= 0x3F0) { // Normal f64, subnormal f16.
56 // Convert from a 53-bit mantissa (after realizing the implicit bit) to a
57 // 10-bit mantissa and then adjust for the exponent.
58 man |= 0x0010000000000000;
Nigel Tao56d90962020-07-12 21:11:49 +100059 uint32_t shift = ((uint32_t)(1051 - exp)); // 1051 = 0x3F0 + 53 - 10.
Nigel Tao7bf7cf22020-07-12 16:23:15 +100060 uint64_t shifted_man = man >> shift;
61 wuffs_base__lossy_value_u16 ret;
62 ret.value = neg | ((uint16_t)shifted_man);
63 ret.lossy = (shifted_man << shift) != man;
64 return ret;
65 }
66
67 // Normal f64, normal f16.
68
69 // Re-bias from 1023 to 15 and shift above f16's 10 mantissa bits.
70 exp = (exp - 1008) << 10; // 1008 = 1023 - 15 = 0x3FF - 0xF.
71
72 // Convert from a 52-bit mantissa (excluding the implicit bit) to a 10-bit
73 // mantissa (again excluding the implicit bit). We lose some information if
74 // any of the bottom 42 bits are non-zero.
75 wuffs_base__lossy_value_u16 ret;
76 ret.value = neg | ((uint16_t)exp) | ((uint16_t)(man >> 42));
77 ret.lossy = (man << 22) != 0;
78 return ret;
79}
80
81WUFFS_BASE__MAYBE_STATIC wuffs_base__lossy_value_u32 //
Nigel Taoa3931d52020-07-12 21:06:44 +100082wuffs_base__ieee_754_bit_representation__from_f64_to_u32_truncate(double f) {
Nigel Tao7bf7cf22020-07-12 16:23:15 +100083 uint64_t u = 0;
84 if (sizeof(uint64_t) == sizeof(double)) {
85 memcpy(&u, &f, sizeof(uint64_t));
86 }
87 uint32_t neg = ((uint32_t)(u >> 63)) << 31;
88 u &= 0x7FFFFFFFFFFFFFFF;
89 uint64_t exp = u >> 52;
90 uint64_t man = u & 0x000FFFFFFFFFFFFF;
91
92 if (exp == 0x7FF) {
93 if (man == 0) { // Infinity.
94 wuffs_base__lossy_value_u32 ret;
95 ret.value = neg | 0x7F800000;
96 ret.lossy = false;
97 return ret;
98 }
99 // NaN. Shift the 52 mantissa bits to 23 mantissa bits, keeping the most
100 // significant mantissa bit (quiet vs signaling NaNs). Also set the low 22
101 // bits of ret.value so that the 23-bit mantissa is non-zero.
102 wuffs_base__lossy_value_u32 ret;
103 ret.value = neg | 0x7FBFFFFF | ((uint32_t)(man >> 29));
104 ret.lossy = false;
105 return ret;
106
107 } else if (exp > 0x47E) { // Truncate to the largest finite f32.
108 wuffs_base__lossy_value_u32 ret;
109 ret.value = neg | 0x7F7FFFFF;
110 ret.lossy = true;
111 return ret;
112
113 } else if (exp <= 0x369) { // Truncate to zero.
114 wuffs_base__lossy_value_u32 ret;
115 ret.value = neg;
116 ret.lossy = (u != 0);
117 return ret;
118
119 } else if (exp <= 0x380) { // Normal f64, subnormal f32.
120 // Convert from a 53-bit mantissa (after realizing the implicit bit) to a
121 // 23-bit mantissa and then adjust for the exponent.
122 man |= 0x0010000000000000;
Nigel Tao56d90962020-07-12 21:11:49 +1000123 uint32_t shift = ((uint32_t)(926 - exp)); // 926 = 0x380 + 53 - 23.
Nigel Tao7bf7cf22020-07-12 16:23:15 +1000124 uint64_t shifted_man = man >> shift;
125 wuffs_base__lossy_value_u32 ret;
126 ret.value = neg | ((uint32_t)shifted_man);
127 ret.lossy = (shifted_man << shift) != man;
128 return ret;
129 }
130
131 // Normal f64, normal f32.
132
133 // Re-bias from 1023 to 127 and shift above f32's 23 mantissa bits.
134 exp = (exp - 896) << 23; // 896 = 1023 - 127 = 0x3FF - 0x7F.
135
136 // Convert from a 52-bit mantissa (excluding the implicit bit) to a 23-bit
137 // mantissa (again excluding the implicit bit). We lose some information if
138 // any of the bottom 29 bits are non-zero.
139 wuffs_base__lossy_value_u32 ret;
140 ret.value = neg | ((uint32_t)exp) | ((uint32_t)(man >> 29));
141 ret.lossy = (man << 35) != 0;
142 return ret;
143}
144
145// --------
146
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000147#define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE 2047
148#define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION 800
149
150// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL is the largest N
151// such that ((10 << N) < (1 << 64)).
152#define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL 60
153
154// wuffs_base__private_implementation__high_prec_dec (abbreviated as HPD) is a
155// fixed precision floating point decimal number, augmented with ±infinity
156// values, but it cannot represent NaN (Not a Number).
157//
158// "High precision" means that the mantissa holds 800 decimal digits. 800 is
159// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION.
160//
161// An HPD isn't for general purpose arithmetic, only for conversions to and
162// from IEEE 754 double-precision floating point, where the largest and
163// smallest positive, finite values are approximately 1.8e+308 and 4.9e-324.
164// HPD exponents above +2047 mean infinity, below -2047 mean zero. The ±2047
165// bounds are further away from zero than ±(324 + 800), where 800 and 2047 is
166// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION and
167// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE.
168//
169// digits[.. num_digits] are the number's digits in big-endian order. The
170// uint8_t values are in the range [0 ..= 9], not ['0' ..= '9'], where e.g. '7'
171// is the ASCII value 0x37.
172//
173// decimal_point is the index (within digits) of the decimal point. It may be
174// negative or be larger than num_digits, in which case the explicit digits are
175// padded with implicit zeroes.
176//
177// For example, if num_digits is 3 and digits is "\x07\x08\x09":
Nigel Taof148f5c2021-11-05 16:50:38 +1100178// - A decimal_point of -2 means ".00789"
179// - A decimal_point of -1 means ".0789"
180// - A decimal_point of +0 means ".789"
181// - A decimal_point of +1 means "7.89"
182// - A decimal_point of +2 means "78.9"
183// - A decimal_point of +3 means "789."
184// - A decimal_point of +4 means "7890."
185// - A decimal_point of +5 means "78900."
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000186//
187// As above, a decimal_point higher than +2047 means that the overall value is
188// infinity, lower than -2047 means zero.
189//
190// negative is a sign bit. An HPD can distinguish positive and negative zero.
191//
192// truncated is whether there are more than
193// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION digits, and at
194// least one of those extra digits are non-zero. The existence of long-tail
195// digits can affect rounding.
196//
197// The "all fields are zero" value is valid, and represents the number +0.
Nigel Tao4f1d24c2020-09-23 22:02:53 +1000198typedef struct wuffs_base__private_implementation__high_prec_dec__struct {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000199 uint32_t num_digits;
200 int32_t decimal_point;
201 bool negative;
202 bool truncated;
203 uint8_t digits[WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION];
204} wuffs_base__private_implementation__high_prec_dec;
205
206// wuffs_base__private_implementation__high_prec_dec__trim trims trailing
207// zeroes from the h->digits[.. h->num_digits] slice. They have no benefit,
208// since we explicitly track h->decimal_point.
209//
210// Preconditions:
211// - h is non-NULL.
212static inline void //
213wuffs_base__private_implementation__high_prec_dec__trim(
214 wuffs_base__private_implementation__high_prec_dec* h) {
215 while ((h->num_digits > 0) && (h->digits[h->num_digits - 1] == 0)) {
216 h->num_digits--;
217 }
218}
219
220// wuffs_base__private_implementation__high_prec_dec__assign sets h to
221// represent the number x.
222//
223// Preconditions:
224// - h is non-NULL.
225static void //
226wuffs_base__private_implementation__high_prec_dec__assign(
227 wuffs_base__private_implementation__high_prec_dec* h,
228 uint64_t x,
229 bool negative) {
230 uint32_t n = 0;
231
232 // Set h->digits.
233 if (x > 0) {
234 // Calculate the digits, working right-to-left. After we determine n (how
235 // many digits there are), copy from buf to h->digits.
236 //
237 // UINT64_MAX, 18446744073709551615, is 20 digits long. It can be faster to
238 // copy a constant number of bytes than a variable number (20 instead of
239 // n). Make buf large enough (and start writing to it from the middle) so
240 // that can we always copy 20 bytes: the slice buf[(20-n) .. (40-n)].
241 uint8_t buf[40] = {0};
242 uint8_t* ptr = &buf[20];
243 do {
244 uint64_t remaining = x / 10;
245 x -= remaining * 10;
246 ptr--;
247 *ptr = (uint8_t)x;
248 n++;
249 x = remaining;
250 } while (x > 0);
251 memcpy(h->digits, ptr, 20);
252 }
253
254 // Set h's other fields.
255 h->num_digits = n;
256 h->decimal_point = (int32_t)n;
257 h->negative = negative;
258 h->truncated = false;
259 wuffs_base__private_implementation__high_prec_dec__trim(h);
260}
261
262static wuffs_base__status //
263wuffs_base__private_implementation__high_prec_dec__parse(
264 wuffs_base__private_implementation__high_prec_dec* h,
Nigel Taoe0c5de92020-07-11 11:48:17 +1000265 wuffs_base__slice_u8 s,
266 uint32_t options) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000267 if (!h) {
268 return wuffs_base__make_status(wuffs_base__error__bad_receiver);
269 }
270 h->num_digits = 0;
271 h->decimal_point = 0;
272 h->negative = false;
273 h->truncated = false;
274
275 uint8_t* p = s.ptr;
276 uint8_t* q = s.ptr + s.len;
277
Nigel Taoc5c98852020-07-11 13:10:14 +1000278 if (options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES) {
279 for (;; p++) {
280 if (p >= q) {
281 return wuffs_base__make_status(wuffs_base__error__bad_argument);
282 } else if (*p != '_') {
283 break;
284 }
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000285 }
286 }
287
288 // Parse sign.
289 do {
290 if (*p == '+') {
291 p++;
292 } else if (*p == '-') {
293 h->negative = true;
294 p++;
295 } else {
296 break;
297 }
Nigel Taoc5c98852020-07-11 13:10:14 +1000298 if (options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES) {
299 for (;; p++) {
300 if (p >= q) {
301 return wuffs_base__make_status(wuffs_base__error__bad_argument);
302 } else if (*p != '_') {
303 break;
304 }
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000305 }
306 }
307 } while (0);
308
309 // Parse digits, up to (and including) a '.', 'E' or 'e'. Examples for each
310 // limb in this if-else chain:
311 // - "0.789"
312 // - "1002.789"
313 // - ".789"
314 // - Other (invalid input).
315 uint32_t nd = 0;
316 int32_t dp = 0;
317 bool no_digits_before_separator = false;
Nigel Taoe82bc8e2020-07-11 12:49:15 +1000318 if (('0' == *p) &&
319 !(options &
320 WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_MULTIPLE_LEADING_ZEROES)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000321 p++;
322 for (;; p++) {
323 if (p >= q) {
324 goto after_all;
Nigel Taoe0c5de92020-07-11 11:48:17 +1000325 } else if (*p ==
326 ((options &
327 WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
328 ? ','
329 : '.')) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000330 p++;
331 goto after_sep;
332 } else if ((*p == 'E') || (*p == 'e')) {
333 p++;
334 goto after_exp;
Nigel Taoc5c98852020-07-11 13:10:14 +1000335 } else if ((*p != '_') ||
336 !(options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000337 return wuffs_base__make_status(wuffs_base__error__bad_argument);
338 }
339 }
340
Nigel Taoe82bc8e2020-07-11 12:49:15 +1000341 } else if (('0' <= *p) && (*p <= '9')) {
342 if (*p == '0') {
343 for (; (p < q) && (*p == '0'); p++) {
344 }
345 } else {
346 h->digits[nd++] = (uint8_t)(*p - '0');
347 dp = (int32_t)nd;
348 p++;
349 }
350
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000351 for (;; p++) {
352 if (p >= q) {
353 goto after_all;
354 } else if (('0' <= *p) && (*p <= '9')) {
355 if (nd < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
356 h->digits[nd++] = (uint8_t)(*p - '0');
357 dp = (int32_t)nd;
358 } else if ('0' != *p) {
359 // Long-tail non-zeroes set the truncated bit.
360 h->truncated = true;
361 }
Nigel Taoe0c5de92020-07-11 11:48:17 +1000362 } else if (*p ==
363 ((options &
364 WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
365 ? ','
366 : '.')) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000367 p++;
368 goto after_sep;
369 } else if ((*p == 'E') || (*p == 'e')) {
370 p++;
371 goto after_exp;
Nigel Taoc5c98852020-07-11 13:10:14 +1000372 } else if ((*p != '_') ||
373 !(options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000374 return wuffs_base__make_status(wuffs_base__error__bad_argument);
375 }
376 }
377
Nigel Taoe0c5de92020-07-11 11:48:17 +1000378 } else if (*p == ((options &
379 WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
380 ? ','
381 : '.')) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000382 p++;
383 no_digits_before_separator = true;
384
385 } else {
386 return wuffs_base__make_status(wuffs_base__error__bad_argument);
387 }
388
389after_sep:
390 for (;; p++) {
391 if (p >= q) {
392 goto after_all;
393 } else if ('0' == *p) {
394 if (nd == 0) {
395 // Track leading zeroes implicitly.
396 dp--;
397 } else if (nd <
398 WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
399 h->digits[nd++] = (uint8_t)(*p - '0');
400 }
401 } else if (('0' < *p) && (*p <= '9')) {
402 if (nd < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
403 h->digits[nd++] = (uint8_t)(*p - '0');
404 } else {
405 // Long-tail non-zeroes set the truncated bit.
406 h->truncated = true;
407 }
408 } else if ((*p == 'E') || (*p == 'e')) {
409 p++;
410 goto after_exp;
Nigel Taoc5c98852020-07-11 13:10:14 +1000411 } else if ((*p != '_') ||
412 !(options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000413 return wuffs_base__make_status(wuffs_base__error__bad_argument);
414 }
415 }
416
417after_exp:
418 do {
Nigel Taoc5c98852020-07-11 13:10:14 +1000419 if (options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES) {
420 for (;; p++) {
421 if (p >= q) {
422 return wuffs_base__make_status(wuffs_base__error__bad_argument);
423 } else if (*p != '_') {
424 break;
425 }
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000426 }
427 }
428
429 int32_t exp_sign = +1;
430 if (*p == '+') {
431 p++;
432 } else if (*p == '-') {
433 exp_sign = -1;
434 p++;
435 }
436
437 int32_t exp = 0;
438 const int32_t exp_large =
439 WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE +
440 WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION;
441 bool saw_exp_digits = false;
442 for (; p < q; p++) {
Nigel Taoc5c98852020-07-11 13:10:14 +1000443 if ((*p == '_') &&
444 (options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000445 // No-op.
446 } else if (('0' <= *p) && (*p <= '9')) {
447 saw_exp_digits = true;
448 if (exp < exp_large) {
449 exp = (10 * exp) + ((int32_t)(*p - '0'));
450 }
451 } else {
452 break;
453 }
454 }
455 if (!saw_exp_digits) {
456 return wuffs_base__make_status(wuffs_base__error__bad_argument);
457 }
458 dp += exp_sign * exp;
459 } while (0);
460
461after_all:
462 if (p != q) {
463 return wuffs_base__make_status(wuffs_base__error__bad_argument);
464 }
465 h->num_digits = nd;
466 if (nd == 0) {
467 if (no_digits_before_separator) {
468 return wuffs_base__make_status(wuffs_base__error__bad_argument);
469 }
470 h->decimal_point = 0;
471 } else if (dp <
472 -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
473 h->decimal_point =
474 -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE - 1;
475 } else if (dp >
476 +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
477 h->decimal_point =
478 +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE + 1;
479 } else {
480 h->decimal_point = dp;
481 }
482 wuffs_base__private_implementation__high_prec_dec__trim(h);
483 return wuffs_base__make_status(NULL);
484}
485
486// --------
487
488// wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits
489// returns the number of additional decimal digits when left-shifting by shift.
490//
491// See below for preconditions.
492static uint32_t //
493wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits(
494 wuffs_base__private_implementation__high_prec_dec* h,
495 uint32_t shift) {
496 // Masking with 0x3F should be unnecessary (assuming the preconditions) but
497 // it's cheap and ensures that we don't overflow the
498 // wuffs_base__private_implementation__hpd_left_shift array.
499 shift &= 63;
500
501 uint32_t x_a = wuffs_base__private_implementation__hpd_left_shift[shift];
502 uint32_t x_b = wuffs_base__private_implementation__hpd_left_shift[shift + 1];
503 uint32_t num_new_digits = x_a >> 11;
504 uint32_t pow5_a = 0x7FF & x_a;
505 uint32_t pow5_b = 0x7FF & x_b;
506
507 const uint8_t* pow5 =
508 &wuffs_base__private_implementation__powers_of_5[pow5_a];
509 uint32_t i = 0;
510 uint32_t n = pow5_b - pow5_a;
511 for (; i < n; i++) {
512 if (i >= h->num_digits) {
513 return num_new_digits - 1;
514 } else if (h->digits[i] == pow5[i]) {
515 continue;
516 } else if (h->digits[i] < pow5[i]) {
517 return num_new_digits - 1;
518 } else {
519 return num_new_digits;
520 }
521 }
522 return num_new_digits;
523}
524
525// --------
526
527// wuffs_base__private_implementation__high_prec_dec__rounded_integer returns
528// the integral (non-fractional) part of h, provided that it is 18 or fewer
529// decimal digits. For 19 or more digits, it returns UINT64_MAX. Note that:
Nigel Taof148f5c2021-11-05 16:50:38 +1100530// - (1 << 53) is 9007199254740992, which has 16 decimal digits.
531// - (1 << 56) is 72057594037927936, which has 17 decimal digits.
532// - (1 << 59) is 576460752303423488, which has 18 decimal digits.
533// - (1 << 63) is 9223372036854775808, which has 19 decimal digits.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000534// and that IEEE 754 double precision has 52 mantissa bits.
535//
536// That integral part is rounded-to-even: rounding 7.5 or 8.5 both give 8.
537//
538// h's negative bit is ignored: rounding -8.6 returns 9.
539//
540// See below for preconditions.
541static uint64_t //
542wuffs_base__private_implementation__high_prec_dec__rounded_integer(
543 wuffs_base__private_implementation__high_prec_dec* h) {
544 if ((h->num_digits == 0) || (h->decimal_point < 0)) {
545 return 0;
546 } else if (h->decimal_point > 18) {
547 return UINT64_MAX;
548 }
549
550 uint32_t dp = (uint32_t)(h->decimal_point);
551 uint64_t n = 0;
552 uint32_t i = 0;
553 for (; i < dp; i++) {
554 n = (10 * n) + ((i < h->num_digits) ? h->digits[i] : 0);
555 }
556
557 bool round_up = false;
558 if (dp < h->num_digits) {
559 round_up = h->digits[dp] >= 5;
560 if ((h->digits[dp] == 5) && (dp + 1 == h->num_digits)) {
561 // We are exactly halfway. If we're truncated, round up, otherwise round
562 // to even.
563 round_up = h->truncated || //
564 ((dp > 0) && (1 & h->digits[dp - 1]));
565 }
566 }
567 if (round_up) {
568 n++;
569 }
570
571 return n;
572}
573
574// wuffs_base__private_implementation__high_prec_dec__small_xshift shifts h's
575// number (where 'x' is 'l' or 'r' for left or right) by a small shift value.
576//
577// Preconditions:
578// - h is non-NULL.
579// - h->decimal_point is "not extreme".
580// - shift is non-zero.
581// - shift is "a small shift".
582//
583// "Not extreme" means within
584// ±WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE.
585//
586// "A small shift" means not more than
587// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL.
588//
589// wuffs_base__private_implementation__high_prec_dec__rounded_integer and
590// wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits
591// have the same preconditions.
592//
593// wuffs_base__private_implementation__high_prec_dec__lshift keeps the first
594// two preconditions but not the last two. Its shift argument is signed and
595// does not need to be "small": zero is a no-op, positive means left shift and
596// negative means right shift.
597
598static void //
599wuffs_base__private_implementation__high_prec_dec__small_lshift(
600 wuffs_base__private_implementation__high_prec_dec* h,
601 uint32_t shift) {
602 if (h->num_digits == 0) {
603 return;
604 }
605 uint32_t num_new_digits =
606 wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits(
607 h, shift);
608 uint32_t rx = h->num_digits - 1; // Read index.
609 uint32_t wx = h->num_digits - 1 + num_new_digits; // Write index.
610 uint64_t n = 0;
611
612 // Repeat: pick up a digit, put down a digit, right to left.
613 while (((int32_t)rx) >= 0) {
614 n += ((uint64_t)(h->digits[rx])) << shift;
615 uint64_t quo = n / 10;
616 uint64_t rem = n - (10 * quo);
617 if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
618 h->digits[wx] = (uint8_t)rem;
619 } else if (rem > 0) {
620 h->truncated = true;
621 }
622 n = quo;
623 wx--;
624 rx--;
625 }
626
627 // Put down leading digits, right to left.
628 while (n > 0) {
629 uint64_t quo = n / 10;
630 uint64_t rem = n - (10 * quo);
631 if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
632 h->digits[wx] = (uint8_t)rem;
633 } else if (rem > 0) {
634 h->truncated = true;
635 }
636 n = quo;
637 wx--;
638 }
639
640 // Finish.
641 h->num_digits += num_new_digits;
642 if (h->num_digits >
643 WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
644 h->num_digits = WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION;
645 }
646 h->decimal_point += (int32_t)num_new_digits;
647 wuffs_base__private_implementation__high_prec_dec__trim(h);
648}
649
650static void //
651wuffs_base__private_implementation__high_prec_dec__small_rshift(
652 wuffs_base__private_implementation__high_prec_dec* h,
653 uint32_t shift) {
654 uint32_t rx = 0; // Read index.
655 uint32_t wx = 0; // Write index.
656 uint64_t n = 0;
657
658 // Pick up enough leading digits to cover the first shift.
659 while ((n >> shift) == 0) {
660 if (rx < h->num_digits) {
661 // Read a digit.
662 n = (10 * n) + h->digits[rx++];
663 } else if (n == 0) {
664 // h's number used to be zero and remains zero.
665 return;
666 } else {
667 // Read sufficient implicit trailing zeroes.
668 while ((n >> shift) == 0) {
669 n = 10 * n;
670 rx++;
671 }
672 break;
673 }
674 }
675 h->decimal_point -= ((int32_t)(rx - 1));
676 if (h->decimal_point <
677 -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
678 // After the shift, h's number is effectively zero.
679 h->num_digits = 0;
680 h->decimal_point = 0;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000681 h->truncated = false;
682 return;
683 }
684
685 // Repeat: pick up a digit, put down a digit, left to right.
686 uint64_t mask = (((uint64_t)(1)) << shift) - 1;
687 while (rx < h->num_digits) {
688 uint8_t new_digit = ((uint8_t)(n >> shift));
689 n = (10 * (n & mask)) + h->digits[rx++];
690 h->digits[wx++] = new_digit;
691 }
692
693 // Put down trailing digits, left to right.
694 while (n > 0) {
695 uint8_t new_digit = ((uint8_t)(n >> shift));
696 n = 10 * (n & mask);
697 if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
698 h->digits[wx++] = new_digit;
699 } else if (new_digit > 0) {
700 h->truncated = true;
701 }
702 }
703
704 // Finish.
705 h->num_digits = wx;
706 wuffs_base__private_implementation__high_prec_dec__trim(h);
707}
708
709static void //
710wuffs_base__private_implementation__high_prec_dec__lshift(
711 wuffs_base__private_implementation__high_prec_dec* h,
712 int32_t shift) {
713 if (shift > 0) {
714 while (shift > +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) {
715 wuffs_base__private_implementation__high_prec_dec__small_lshift(
716 h, WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL);
717 shift -= WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
718 }
719 wuffs_base__private_implementation__high_prec_dec__small_lshift(
720 h, ((uint32_t)(+shift)));
721 } else if (shift < 0) {
722 while (shift < -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) {
723 wuffs_base__private_implementation__high_prec_dec__small_rshift(
724 h, WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL);
725 shift += WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
726 }
727 wuffs_base__private_implementation__high_prec_dec__small_rshift(
728 h, ((uint32_t)(-shift)));
729 }
730}
731
732// --------
733
734// wuffs_base__private_implementation__high_prec_dec__round_etc rounds h's
735// number. For those functions that take an n argument, rounding produces at
736// most n digits (which is not necessarily at most n decimal places). Negative
737// n values are ignored, as well as any n greater than or equal to h's number
738// of digits. The etc__round_just_enough function implicitly chooses an n to
739// implement WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION.
740//
741// Preconditions:
742// - h is non-NULL.
743// - h->decimal_point is "not extreme".
744//
745// "Not extreme" means within
746// ±WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE.
747
748static void //
749wuffs_base__private_implementation__high_prec_dec__round_down(
750 wuffs_base__private_implementation__high_prec_dec* h,
751 int32_t n) {
752 if ((n < 0) || (h->num_digits <= (uint32_t)n)) {
753 return;
754 }
755 h->num_digits = (uint32_t)(n);
756 wuffs_base__private_implementation__high_prec_dec__trim(h);
757}
758
759static void //
760wuffs_base__private_implementation__high_prec_dec__round_up(
761 wuffs_base__private_implementation__high_prec_dec* h,
762 int32_t n) {
763 if ((n < 0) || (h->num_digits <= (uint32_t)n)) {
764 return;
765 }
766
767 for (n--; n >= 0; n--) {
768 if (h->digits[n] < 9) {
769 h->digits[n]++;
770 h->num_digits = (uint32_t)(n + 1);
771 return;
772 }
773 }
774
775 // The number is all 9s. Change to a single 1 and adjust the decimal point.
776 h->digits[0] = 1;
777 h->num_digits = 1;
778 h->decimal_point++;
779}
780
781static void //
782wuffs_base__private_implementation__high_prec_dec__round_nearest(
783 wuffs_base__private_implementation__high_prec_dec* h,
784 int32_t n) {
785 if ((n < 0) || (h->num_digits <= (uint32_t)n)) {
786 return;
787 }
788 bool up = h->digits[n] >= 5;
789 if ((h->digits[n] == 5) && ((n + 1) == ((int32_t)(h->num_digits)))) {
790 up = h->truncated || //
791 ((n > 0) && ((h->digits[n - 1] & 1) != 0));
792 }
793
794 if (up) {
795 wuffs_base__private_implementation__high_prec_dec__round_up(h, n);
796 } else {
797 wuffs_base__private_implementation__high_prec_dec__round_down(h, n);
798 }
799}
800
801static void //
802wuffs_base__private_implementation__high_prec_dec__round_just_enough(
803 wuffs_base__private_implementation__high_prec_dec* h,
804 int32_t exp2,
805 uint64_t mantissa) {
806 // The magic numbers 52 and 53 in this function are because IEEE 754 double
807 // precision has 52 mantissa bits.
808 //
809 // Let f be the floating point number represented by exp2 and mantissa (and
810 // also the number in h): the number (mantissa * (2 ** (exp2 - 52))).
811 //
812 // If f is zero or a small integer, we can return early.
813 if ((mantissa == 0) ||
814 ((exp2 < 53) && (h->decimal_point >= ((int32_t)(h->num_digits))))) {
815 return;
816 }
817
818 // The smallest normal f has an exp2 of -1022 and a mantissa of (1 << 52).
819 // Subnormal numbers have the same exp2 but a smaller mantissa.
820 static const int32_t min_incl_normal_exp2 = -1022;
821 static const uint64_t min_incl_normal_mantissa = 0x0010000000000000ul;
822
823 // Compute lower and upper bounds such that any number between them (possibly
824 // inclusive) will round to f. First, the lower bound. Our number f is:
825 // ((mantissa + 0) * (2 ** ( exp2 - 52)))
826 //
827 // The next lowest floating point number is:
828 // ((mantissa - 1) * (2 ** ( exp2 - 52)))
829 // unless (mantissa - 1) drops the (1 << 52) bit and exp2 is not the
830 // min_incl_normal_exp2. Either way, call it:
831 // ((l_mantissa) * (2 ** (l_exp2 - 52)))
832 //
833 // The lower bound is halfway between them (noting that 52 became 53):
834 // (((2 * l_mantissa) + 1) * (2 ** (l_exp2 - 53)))
835 int32_t l_exp2 = exp2;
836 uint64_t l_mantissa = mantissa - 1;
837 if ((exp2 > min_incl_normal_exp2) && (mantissa <= min_incl_normal_mantissa)) {
838 l_exp2 = exp2 - 1;
839 l_mantissa = (2 * mantissa) - 1;
840 }
841 wuffs_base__private_implementation__high_prec_dec lower;
842 wuffs_base__private_implementation__high_prec_dec__assign(
843 &lower, (2 * l_mantissa) + 1, false);
844 wuffs_base__private_implementation__high_prec_dec__lshift(&lower,
845 l_exp2 - 53);
846
847 // Next, the upper bound. Our number f is:
848 // ((mantissa + 0) * (2 ** (exp2 - 52)))
849 //
850 // The next highest floating point number is:
851 // ((mantissa + 1) * (2 ** (exp2 - 52)))
852 //
853 // The upper bound is halfway between them (noting that 52 became 53):
854 // (((2 * mantissa) + 1) * (2 ** (exp2 - 53)))
855 wuffs_base__private_implementation__high_prec_dec upper;
856 wuffs_base__private_implementation__high_prec_dec__assign(
857 &upper, (2 * mantissa) + 1, false);
858 wuffs_base__private_implementation__high_prec_dec__lshift(&upper, exp2 - 53);
859
860 // The lower and upper bounds are possible outputs only if the original
861 // mantissa is even, so that IEEE round-to-even would round to the original
862 // mantissa and not its neighbors.
863 bool inclusive = (mantissa & 1) == 0;
864
865 // As we walk the digits, we want to know whether rounding up would fall
866 // within the upper bound. This is tracked by upper_delta:
867 // - When -1, the digits of h and upper are the same so far.
868 // - When +0, we saw a difference of 1 between h and upper on a previous
869 // digit and subsequently only 9s for h and 0s for upper. Thus, rounding
870 // up may fall outside of the bound if !inclusive.
871 // - When +1, the difference is greater than 1 and we know that rounding up
872 // falls within the bound.
873 //
874 // This is a state machine with three states. The numerical value for each
875 // state (-1, +0 or +1) isn't important, other than their order.
876 int upper_delta = -1;
877
878 // We can now figure out the shortest number of digits required. Walk the
879 // digits until h has distinguished itself from lower or upper.
880 //
881 // The zi and zd variables are indexes and digits, for z in l (lower), h (the
882 // number) and u (upper).
883 //
884 // The lower, h and upper numbers may have their decimal points at different
885 // places. In this case, upper is the longest, so we iterate ui starting from
886 // 0 and iterate li and hi starting from either 0 or -1.
887 int32_t ui = 0;
888 for (;; ui++) {
889 // Calculate hd, the middle number's digit.
890 int32_t hi = ui - upper.decimal_point + h->decimal_point;
891 if (hi >= ((int32_t)(h->num_digits))) {
892 break;
893 }
894 uint8_t hd = (((uint32_t)hi) < h->num_digits) ? h->digits[hi] : 0;
895
896 // Calculate ld, the lower bound's digit.
897 int32_t li = ui - upper.decimal_point + lower.decimal_point;
898 uint8_t ld = (((uint32_t)li) < lower.num_digits) ? lower.digits[li] : 0;
899
900 // We can round down (truncate) if lower has a different digit than h or if
901 // lower is inclusive and is exactly the result of rounding down (i.e. we
902 // have reached the final digit of lower).
903 bool can_round_down =
904 (ld != hd) || //
905 (inclusive && ((li + 1) == ((int32_t)(lower.num_digits))));
906
907 // Calculate ud, the upper bound's digit, and update upper_delta.
908 uint8_t ud = (((uint32_t)ui) < upper.num_digits) ? upper.digits[ui] : 0;
909 if (upper_delta < 0) {
910 if ((hd + 1) < ud) {
911 // For example:
912 // h = 12345???
913 // upper = 12347???
914 upper_delta = +1;
915 } else if (hd != ud) {
916 // For example:
917 // h = 12345???
918 // upper = 12346???
919 upper_delta = +0;
920 }
921 } else if (upper_delta == 0) {
922 if ((hd != 9) || (ud != 0)) {
923 // For example:
924 // h = 1234598?
925 // upper = 1234600?
926 upper_delta = +1;
927 }
928 }
929
930 // We can round up if upper has a different digit than h and either upper
931 // is inclusive or upper is bigger than the result of rounding up.
932 bool can_round_up =
933 (upper_delta > 0) || //
934 ((upper_delta == 0) && //
935 (inclusive || ((ui + 1) < ((int32_t)(upper.num_digits)))));
936
937 // If we can round either way, round to nearest. If we can round only one
938 // way, do it. If we can't round, continue the loop.
939 if (can_round_down) {
940 if (can_round_up) {
941 wuffs_base__private_implementation__high_prec_dec__round_nearest(
942 h, hi + 1);
943 return;
944 } else {
945 wuffs_base__private_implementation__high_prec_dec__round_down(h,
946 hi + 1);
947 return;
948 }
949 } else {
950 if (can_round_up) {
951 wuffs_base__private_implementation__high_prec_dec__round_up(h, hi + 1);
952 return;
953 }
954 }
955 }
956}
957
958// --------
959
Nigel Taoc4fa8e22020-07-18 17:35:13 +1000960// wuffs_base__private_implementation__parse_number_f64_eisel_lemire produces
961// the IEEE 754 double-precision value for an exact mantissa and base-10
962// exponent. For example:
Nigel Taob15a0fc2020-07-08 10:50:14 +1000963// - when parsing "12345.678e+02", man is 12345678 and exp10 is -1.
964// - when parsing "-12", man is 12 and exp10 is 0. Processing the leading
965// minus sign is the responsibility of the caller, not this function.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000966//
967// On success, it returns a non-negative int64_t such that the low 63 bits hold
968// the 11-bit exponent and 52-bit mantissa.
969//
970// On failure, it returns a negative value.
971//
Nigel Taoc4fa8e22020-07-18 17:35:13 +1000972// The algorithm is based on an original idea by Michael Eisel that was refined
973// by Daniel Lemire. See
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000974// https://lemire.me/blog/2020/03/10/fast-float-parsing-in-practice/
Nigel Tao1d8d18f2020-10-07 22:13:51 +1100975// and
976// https://nigeltao.github.io/blog/2020/eisel-lemire.html
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000977//
978// Preconditions:
979// - man is non-zero.
Nigel Tao09872962020-09-15 22:22:51 +1000980// - exp10 is in the range [-307 ..= 288], the same range of the
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000981// wuffs_base__private_implementation__powers_of_10 array.
Nigel Tao8b45db02020-09-15 21:50:32 +1000982//
983// The exp10 range (and the fact that man is in the range [1 ..= UINT64_MAX],
984// approximately [1 ..= 1.85e+19]) means that (man * (10 ** exp10)) is in the
985// range [1e-307 ..= 1.85e+307]. This is entirely within the range of normal
986// (neither subnormal nor non-finite) f64 values: DBL_MIN and DBL_MAX are
987// approximately 2.23e–308 and 1.80e+308.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000988static int64_t //
Nigel Taoc4fa8e22020-07-18 17:35:13 +1000989wuffs_base__private_implementation__parse_number_f64_eisel_lemire(
990 uint64_t man,
991 int32_t exp10) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000992 // Look up the (possibly truncated) base-2 representation of (10 ** exp10).
993 // The look-up table was constructed so that it is already normalized: the
994 // table entry's mantissa's MSB (most significant bit) is on.
Nigel Taoafe7f272020-09-23 15:52:13 +1000995 const uint64_t* po10 =
996 &wuffs_base__private_implementation__powers_of_10[exp10 + 307][0];
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000997
998 // Normalize the man argument. The (man != 0) precondition means that a
999 // non-zero bit exists.
1000 uint32_t clz = wuffs_base__count_leading_zeroes_u64(man);
1001 man <<= clz;
1002
1003 // Calculate the return value's base-2 exponent. We might tweak it by ±1
Nigel Taob6d85522020-09-23 15:21:47 +10001004 // later, but its initial value comes from a linear scaling of exp10,
1005 // converting from power-of-10 to power-of-2, and adjusting by clz.
1006 //
1007 // The magic constants are:
1008 // - 1087 = 1023 + 64. The 1023 is the f64 exponent bias. The 64 is because
1009 // the look-up table uses 64-bit mantissas.
1010 // - 217706 is such that the ratio 217706 / 65536 ≈ 3.321930 is close enough
1011 // (over the practical range of exp10) to log(10) / log(2) ≈ 3.321928.
1012 // - 65536 = 1<<16 is arbitrary but a power of 2, so division is a shift.
1013 //
1014 // Equality of the linearly-scaled value and the actual power-of-2, over the
1015 // range of exp10 arguments that this function accepts, is confirmed by
1016 // script/print-mpb-powers-of-10.go
1017 uint64_t ret_exp2 =
1018 ((uint64_t)(((217706 * exp10) >> 16) + 1087)) - ((uint64_t)clz);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001019
1020 // Multiply the two mantissas. Normalization means that both mantissas are at
1021 // least (1<<63), so the 128-bit product must be at least (1<<126). The high
Nigel Tao74d4af62020-07-10 11:27:17 +10001022 // 64 bits of the product, x_hi, must therefore be at least (1<<62).
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001023 //
Nigel Tao74d4af62020-07-10 11:27:17 +10001024 // As a consequence, x_hi has either 0 or 1 leading zeroes. Shifting x_hi
1025 // right by either 9 or 10 bits (depending on x_hi's MSB) will therefore
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001026 // leave the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on.
Nigel Taoafe7f272020-09-23 15:52:13 +10001027 wuffs_base__multiply_u64__output x = wuffs_base__multiply_u64(man, po10[1]);
Nigel Tao74d4af62020-07-10 11:27:17 +10001028 uint64_t x_hi = x.hi;
1029 uint64_t x_lo = x.lo;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001030
1031 // Before we shift right by at least 9 bits, recall that the look-up table
1032 // entry was possibly truncated. We have so far only calculated a lower bound
1033 // for the product (man * e), where e is (10 ** exp10). The upper bound would
1034 // add a further (man * 1) to the 128-bit product, which overflows the lower
Nigel Tao74d4af62020-07-10 11:27:17 +10001035 // 64-bit limb if ((x_lo + man) < man).
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001036 //
Nigel Tao74d4af62020-07-10 11:27:17 +10001037 // If overflow occurs, that adds 1 to x_hi. Since we're about to shift right
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001038 // by at least 9 bits, that carried 1 can be ignored unless the higher 64-bit
1039 // limb's low 9 bits are all on.
Nigel Taoba3818c2020-09-28 12:51:45 +10001040 //
1041 // For example, parsing "9999999999999999999" will take the if-true branch
1042 // here, since:
1043 // - x_hi = 0x4563918244F3FFFF
1044 // - x_lo = 0x8000000000000000
1045 // - man = 0x8AC7230489E7FFFF
Nigel Tao74d4af62020-07-10 11:27:17 +10001046 if (((x_hi & 0x1FF) == 0x1FF) && ((x_lo + man) < man)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001047 // Refine our calculation of (man * e). Before, our approximation of e used
1048 // a "low resolution" 64-bit mantissa. Now use a "high resolution" 128-bit
1049 // mantissa. We've already calculated x = (man * bits_0_to_63_incl_of_e).
1050 // Now calculate y = (man * bits_64_to_127_incl_of_e).
Nigel Taoafe7f272020-09-23 15:52:13 +10001051 wuffs_base__multiply_u64__output y = wuffs_base__multiply_u64(man, po10[0]);
Nigel Tao74d4af62020-07-10 11:27:17 +10001052 uint64_t y_hi = y.hi;
1053 uint64_t y_lo = y.lo;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001054
1055 // Merge the 128-bit x and 128-bit y, which overlap by 64 bits, to
1056 // calculate the 192-bit product of the 64-bit man by the 128-bit e.
1057 // As we exit this if-block, we only care about the high 128 bits
1058 // (merged_hi and merged_lo) of that 192-bit product.
Nigel Taoba3818c2020-09-28 12:51:45 +10001059 //
1060 // For example, parsing "1.234e-45" will take the if-true branch here,
1061 // since:
1062 // - x_hi = 0x70B7E3696DB29FFF
1063 // - x_lo = 0xE040000000000000
1064 // - y_hi = 0x33718BBEAB0E0D7A
1065 // - y_lo = 0xA880000000000000
Nigel Tao74d4af62020-07-10 11:27:17 +10001066 uint64_t merged_hi = x_hi;
1067 uint64_t merged_lo = x_lo + y_hi;
1068 if (merged_lo < x_lo) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001069 merged_hi++; // Carry the overflow bit.
1070 }
1071
1072 // The "high resolution" approximation of e is still a lower bound. Once
1073 // again, see if the upper bound is large enough to produce a different
1074 // result. This time, if it does, give up instead of reaching for an even
1075 // more precise approximation to e.
1076 //
1077 // This three-part check is similar to the two-part check that guarded the
1078 // if block that we're now in, but it has an extra term for the middle 64
1079 // bits (checking that adding 1 to merged_lo would overflow).
Nigel Taoba3818c2020-09-28 12:51:45 +10001080 //
1081 // For example, parsing "5.9604644775390625e-8" will take the if-true
1082 // branch here, since:
1083 // - merged_hi = 0x7FFFFFFFFFFFFFFF
1084 // - merged_lo = 0xFFFFFFFFFFFFFFFF
1085 // - y_lo = 0x4DB3FFC120988200
1086 // - man = 0xD3C21BCECCEDA100
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001087 if (((merged_hi & 0x1FF) == 0x1FF) && ((merged_lo + 1) == 0) &&
Nigel Tao74d4af62020-07-10 11:27:17 +10001088 (y_lo + man < man)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001089 return -1;
1090 }
1091
1092 // Replace the 128-bit x with merged.
Nigel Tao74d4af62020-07-10 11:27:17 +10001093 x_hi = merged_hi;
1094 x_lo = merged_lo;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001095 }
1096
Nigel Tao74d4af62020-07-10 11:27:17 +10001097 // As mentioned above, shifting x_hi right by either 9 or 10 bits will leave
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001098 // the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on. If the
1099 // MSB (before shifting) was on, adjust ret_exp2 for the larger shift.
1100 //
1101 // Having bit 53 on (and higher bits off) means that ret_mantissa is a 54-bit
1102 // number.
Nigel Tao74d4af62020-07-10 11:27:17 +10001103 uint64_t msb = x_hi >> 63;
1104 uint64_t ret_mantissa = x_hi >> (msb + 9);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001105 ret_exp2 -= 1 ^ msb;
1106
1107 // IEEE 754 rounds to-nearest with ties rounded to-even. Rounding to-even can
1108 // be tricky. If we're half-way between two exactly representable numbers
1109 // (x's low 73 bits are zero and the next 2 bits that matter are "01"), give
1110 // up instead of trying to pick the winner.
1111 //
1112 // Technically, we could tighten the condition by changing "73" to "73 or 74,
1113 // depending on msb", but a flat "73" is simpler.
Nigel Taoba3818c2020-09-28 12:51:45 +10001114 //
1115 // For example, parsing "1e+23" will take the if-true branch here, since:
1116 // - x_hi = 0x54B40B1F852BDA00
1117 // - ret_mantissa = 0x002A5A058FC295ED
Nigel Tao74d4af62020-07-10 11:27:17 +10001118 if ((x_lo == 0) && ((x_hi & 0x1FF) == 0) && ((ret_mantissa & 3) == 1)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001119 return -1;
1120 }
1121
1122 // If we're not halfway then it's rounding to-nearest. Starting with a 54-bit
1123 // number, carry the lowest bit (bit 0) up if it's on. Regardless of whether
1124 // it was on or off, shifting right by one then produces a 53-bit number. If
1125 // carrying up overflowed, shift again.
1126 ret_mantissa += ret_mantissa & 1;
1127 ret_mantissa >>= 1;
Nigel Tao8b45db02020-09-15 21:50:32 +10001128 // This if block is equivalent to (but benchmarks slightly faster than) the
1129 // following branchless form:
1130 // uint64_t overflow_adjustment = ret_mantissa >> 53;
1131 // ret_mantissa >>= overflow_adjustment;
1132 // ret_exp2 += overflow_adjustment;
Nigel Taoba3818c2020-09-28 12:51:45 +10001133 //
1134 // For example, parsing "7.2057594037927933e+16" will take the if-true
1135 // branch here, since:
1136 // - x_hi = 0x7FFFFFFFFFFFFE80
1137 // - ret_mantissa = 0x0020000000000000
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001138 if ((ret_mantissa >> 53) > 0) {
1139 ret_mantissa >>= 1;
1140 ret_exp2++;
1141 }
1142
1143 // Starting with a 53-bit number, IEEE 754 double-precision normal numbers
1144 // have an implicit mantissa bit. Mask that away and keep the low 52 bits.
1145 ret_mantissa &= 0x000FFFFFFFFFFFFF;
1146
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001147 // Pack the bits and return.
1148 return ((int64_t)(ret_mantissa | (ret_exp2 << 52)));
1149}
1150
1151// --------
1152
1153static wuffs_base__result_f64 //
Nigel Taoe0c5de92020-07-11 11:48:17 +10001154wuffs_base__private_implementation__parse_number_f64_special(
1155 wuffs_base__slice_u8 s,
Nigel Tao4d61a052020-07-11 12:34:40 +10001156 uint32_t options) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001157 do {
Nigel Tao4d61a052020-07-11 12:34:40 +10001158 if (options & WUFFS_BASE__PARSE_NUMBER_FXX__REJECT_INF_AND_NAN) {
1159 goto fail;
1160 }
1161
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001162 uint8_t* p = s.ptr;
1163 uint8_t* q = s.ptr + s.len;
1164
1165 for (; (p < q) && (*p == '_'); p++) {
1166 }
1167 if (p >= q) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001168 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001169 }
1170
1171 // Parse sign.
1172 bool negative = false;
1173 do {
1174 if (*p == '+') {
1175 p++;
1176 } else if (*p == '-') {
1177 negative = true;
1178 p++;
1179 } else {
1180 break;
1181 }
1182 for (; (p < q) && (*p == '_'); p++) {
1183 }
1184 } while (0);
1185 if (p >= q) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001186 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001187 }
1188
1189 bool nan = false;
1190 switch (p[0]) {
1191 case 'I':
1192 case 'i':
1193 if (((q - p) < 3) || //
1194 ((p[1] != 'N') && (p[1] != 'n')) || //
1195 ((p[2] != 'F') && (p[2] != 'f'))) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001196 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001197 }
1198 p += 3;
1199
1200 if ((p >= q) || (*p == '_')) {
1201 break;
1202 } else if (((q - p) < 5) || //
1203 ((p[0] != 'I') && (p[0] != 'i')) || //
1204 ((p[1] != 'N') && (p[1] != 'n')) || //
1205 ((p[2] != 'I') && (p[2] != 'i')) || //
1206 ((p[3] != 'T') && (p[3] != 't')) || //
1207 ((p[4] != 'Y') && (p[4] != 'y'))) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001208 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001209 }
1210 p += 5;
1211
1212 if ((p >= q) || (*p == '_')) {
1213 break;
1214 }
Nigel Tao4d61a052020-07-11 12:34:40 +10001215 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001216
1217 case 'N':
1218 case 'n':
1219 if (((q - p) < 3) || //
1220 ((p[1] != 'A') && (p[1] != 'a')) || //
1221 ((p[2] != 'N') && (p[2] != 'n'))) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001222 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001223 }
1224 p += 3;
1225
1226 if ((p >= q) || (*p == '_')) {
1227 nan = true;
1228 break;
1229 }
Nigel Tao4d61a052020-07-11 12:34:40 +10001230 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001231
1232 default:
Nigel Tao4d61a052020-07-11 12:34:40 +10001233 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001234 }
1235
1236 // Finish.
1237 for (; (p < q) && (*p == '_'); p++) {
1238 }
1239 if (p != q) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001240 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001241 }
1242 wuffs_base__result_f64 ret;
1243 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001244 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001245 (nan ? 0x7FFFFFFFFFFFFFFF : 0x7FF0000000000000) |
1246 (negative ? 0x8000000000000000 : 0));
1247 return ret;
1248 } while (0);
1249
Nigel Tao4d61a052020-07-11 12:34:40 +10001250fail:
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001251 do {
1252 wuffs_base__result_f64 ret;
Nigel Tao4d61a052020-07-11 12:34:40 +10001253 ret.status.repr = wuffs_base__error__bad_argument;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001254 ret.value = 0;
1255 return ret;
1256 } while (0);
1257}
1258
1259WUFFS_BASE__MAYBE_STATIC wuffs_base__result_f64 //
Nigel Taoe0c5de92020-07-11 11:48:17 +10001260wuffs_base__private_implementation__high_prec_dec__to_f64(
Nigel Tao4d61a052020-07-11 12:34:40 +10001261 wuffs_base__private_implementation__high_prec_dec* h,
1262 uint32_t options) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001263 do {
1264 // powers converts decimal powers of 10 to binary powers of 2. For example,
1265 // (10000 >> 13) is 1. It stops before the elements exceed 60, also known
1266 // as WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL.
1267 static const uint32_t num_powers = 19;
1268 static const uint8_t powers[19] = {
1269 0, 3, 6, 9, 13, 16, 19, 23, 26, 29, //
1270 33, 36, 39, 43, 46, 49, 53, 56, 59, //
1271 };
1272
1273 // Handle zero and obvious extremes. The largest and smallest positive
1274 // finite f64 values are approximately 1.8e+308 and 4.9e-324.
1275 if ((h->num_digits == 0) || (h->decimal_point < -326)) {
1276 goto zero;
1277 } else if (h->decimal_point > 310) {
1278 goto infinity;
1279 }
1280
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001281 // Try the fast Eisel-Lemire algorithm again. Calculating the (man, exp10)
1282 // pair from the high_prec_dec h is more correct but slower than the
1283 // approach taken in wuffs_base__parse_number_f64. The latter is optimized
1284 // for the common cases (e.g. assuming no underscores or a leading '+'
1285 // sign) rather than the full set of cases allowed by the Wuffs API.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001286 if (h->num_digits <= 19) {
1287 uint64_t man = 0;
1288 uint32_t i;
1289 for (i = 0; i < h->num_digits; i++) {
1290 man = (10 * man) + h->digits[i];
1291 }
1292 int32_t exp10 = h->decimal_point - ((int32_t)(h->num_digits));
Nigel Tao8b45db02020-09-15 21:50:32 +10001293 if ((man != 0) && (-307 <= exp10) && (exp10 <= 288)) {
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001294 int64_t r =
1295 wuffs_base__private_implementation__parse_number_f64_eisel_lemire(
1296 man, exp10);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001297 if (r >= 0) {
1298 wuffs_base__result_f64 ret;
1299 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001300 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001301 ((uint64_t)r) | (((uint64_t)(h->negative)) << 63));
1302 return ret;
1303 }
1304 }
1305 }
1306
Nigel Taoce685a62020-11-03 15:24:02 +11001307 // When Eisel-Lemire fails, fall back to Simple Decimal Conversion. See
1308 // https://nigeltao.github.io/blog/2020/parse-number-f64-simple.html
1309 //
sarastro-nld46220c2022-11-21 21:29:03 +01001310 // Scale by powers of 2 until we're in the range [1 .. 2]. First we shift
1311 // right, possibly a little too far, ending with a value certainly below
1312 // 10 ...
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001313 const int32_t f64_bias = -1023;
1314 int32_t exp2 = 0;
sarastro-nld46220c2022-11-21 21:29:03 +01001315 while (h->decimal_point > 1) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001316 uint32_t n = (uint32_t)(+h->decimal_point);
1317 uint32_t shift =
1318 (n < num_powers)
1319 ? powers[n]
1320 : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1321
1322 wuffs_base__private_implementation__high_prec_dec__small_rshift(h, shift);
1323 if (h->decimal_point <
1324 -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
1325 goto zero;
1326 }
1327 exp2 += (int32_t)shift;
1328 }
sarastro-nld46220c2022-11-21 21:29:03 +01001329 // ... then we shift left, putting us in [0.1 .. 10].
1330 while (h->decimal_point < 0) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001331 uint32_t shift;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001332 uint32_t n = (uint32_t)(-h->decimal_point);
1333 shift = (n < num_powers)
sarastro-nld46220c2022-11-21 21:29:03 +01001334 ? powers[n] + 1 // shift 1 extra since we're now multiplying
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001335 : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001336
1337 wuffs_base__private_implementation__high_prec_dec__small_lshift(h, shift);
1338 if (h->decimal_point >
1339 +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
1340 goto infinity;
1341 }
1342 exp2 -= (int32_t)shift;
1343 }
1344
sarastro-nld46220c2022-11-21 21:29:03 +01001345 // To get into the range [1 .. 2] we need to look at the first 3 digits of
1346 // the mantisse to determine the final left shift. If we are already between 1
1347 // and 2 then just shift left by 52.
1348 uint64_t man = 0;
1349 uint32_t i;
1350 for (i = 0; i < 3; i++) {
1351 man = (10 * man) + ((i < h->num_digits) ? h->digits[i] : 0);
1352 }
1353 int32_t final_lshift = 52;
1354 int32_t additional_shift = 0;
1355 if (h->decimal_point == 0) { // value is in [0.1 .. 1]
1356 if (man < 125) {
1357 additional_shift = -4;
1358 } else if (man < 250) {
1359 additional_shift = -3;
1360 } else if (man < 500) {
1361 additional_shift = -2;
1362 } else {
1363 additional_shift = -1;
1364 }
1365 } else { // value is in [1 .. 10]
1366 if (man < 200) {
1367 additional_shift = 0;
1368 } else if (man < 400) {
1369 additional_shift = 1;
1370 } else if (man < 800) {
1371 additional_shift = 2;
1372 } else {
1373 additional_shift = 3;
1374 }
1375 }
1376 exp2 += additional_shift;
1377 final_lshift -= additional_shift;
1378
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001379 // The minimum normal exponent is (f64_bias + 1).
1380 while ((f64_bias + 1) > exp2) {
1381 uint32_t n = (uint32_t)((f64_bias + 1) - exp2);
1382 if (n > WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) {
1383 n = WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1384 }
1385 wuffs_base__private_implementation__high_prec_dec__small_rshift(h, n);
1386 exp2 += (int32_t)n;
1387 }
1388
1389 // Check for overflow.
1390 if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1.
1391 goto infinity;
1392 }
1393
1394 // Extract 53 bits for the mantissa (in base-2).
sarastro-nld46220c2022-11-21 21:29:03 +01001395 wuffs_base__private_implementation__high_prec_dec__small_lshift(h, final_lshift);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001396 uint64_t man2 =
1397 wuffs_base__private_implementation__high_prec_dec__rounded_integer(h);
1398
1399 // Rounding might have added one bit. If so, shift and re-check overflow.
1400 if ((man2 >> 53) != 0) {
1401 man2 >>= 1;
1402 exp2++;
1403 if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1.
1404 goto infinity;
1405 }
1406 }
1407
1408 // Handle subnormal numbers.
1409 if ((man2 >> 52) == 0) {
1410 exp2 = f64_bias;
1411 }
1412
1413 // Pack the bits and return.
1414 uint64_t exp2_bits =
1415 (uint64_t)((exp2 - f64_bias) & 0x07FF); // (1 << 11) - 1.
1416 uint64_t bits = (man2 & 0x000FFFFFFFFFFFFF) | // (1 << 52) - 1.
1417 (exp2_bits << 52) | //
1418 (h->negative ? 0x8000000000000000 : 0); // (1 << 63).
1419
1420 wuffs_base__result_f64 ret;
1421 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001422 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001423 return ret;
1424 } while (0);
1425
1426zero:
1427 do {
1428 uint64_t bits = h->negative ? 0x8000000000000000 : 0;
1429
1430 wuffs_base__result_f64 ret;
1431 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001432 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001433 return ret;
1434 } while (0);
1435
1436infinity:
1437 do {
Nigel Tao4d61a052020-07-11 12:34:40 +10001438 if (options & WUFFS_BASE__PARSE_NUMBER_FXX__REJECT_INF_AND_NAN) {
1439 wuffs_base__result_f64 ret;
1440 ret.status.repr = wuffs_base__error__bad_argument;
1441 ret.value = 0;
1442 return ret;
1443 }
1444
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001445 uint64_t bits = h->negative ? 0xFFF0000000000000 : 0x7FF0000000000000;
1446
1447 wuffs_base__result_f64 ret;
1448 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001449 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001450 return ret;
1451 } while (0);
1452}
1453
1454static inline bool //
1455wuffs_base__private_implementation__is_decimal_digit(uint8_t c) {
1456 return ('0' <= c) && (c <= '9');
1457}
1458
1459WUFFS_BASE__MAYBE_STATIC wuffs_base__result_f64 //
1460wuffs_base__parse_number_f64(wuffs_base__slice_u8 s, uint32_t options) {
1461 // In practice, almost all "dd.ddddE±xxx" numbers can be represented
1462 // losslessly by a uint64_t mantissa "dddddd" and an int32_t base-10
1463 // exponent, adjusting "xxx" for the position (if present) of the decimal
1464 // separator '.' or ','.
1465 //
1466 // This (u64 man, i32 exp10) data structure is superficially similar to the
1467 // "Do It Yourself Floating Point" type from Loitsch (†), but the exponent
1468 // here is base-10, not base-2.
1469 //
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001470 // If s's number fits in a (man, exp10), parse that pair with the
1471 // Eisel-Lemire algorithm. If not, or if Eisel-Lemire fails, parsing s with
1472 // the fallback algorithm is slower but comprehensive.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001473 //
1474 // † "Printing Floating-Point Numbers Quickly and Accurately with Integers"
1475 // (https://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf).
1476 // Florian Loitsch is also the primary contributor to
1477 // https://github.com/google/double-conversion
1478 do {
1479 // Calculating that (man, exp10) pair needs to stay within s's bounds.
1480 // Provided that s isn't extremely long, work on a NUL-terminated copy of
1481 // s's contents. The NUL byte isn't a valid part of "±dd.ddddE±xxx".
1482 //
1483 // As the pointer p walks the contents, it's faster to repeatedly check "is
1484 // *p a valid digit" than "is p within bounds and *p a valid digit".
1485 if (s.len >= 256) {
1486 goto fallback;
1487 }
1488 uint8_t z[256];
1489 memcpy(&z[0], s.ptr, s.len);
1490 z[s.len] = 0;
1491 const uint8_t* p = &z[0];
1492
1493 // Look for a leading minus sign. Technically, we could also look for an
1494 // optional plus sign, but the "script/process-json-numbers.c with -p"
1495 // benchmark is noticably slower if we do. It's optional and, in practice,
1496 // usually absent. Let the fallback catch it.
1497 bool negative = (*p == '-');
1498 if (negative) {
1499 p++;
1500 }
1501
1502 // After walking "dd.dddd", comparing p later with p now will produce the
1503 // number of "d"s and "."s.
1504 const uint8_t* const start_of_digits_ptr = p;
1505
1506 // Walk the "d"s before a '.', 'E', NUL byte, etc. If it starts with '0',
1507 // it must be a single '0'. If it starts with a non-zero decimal digit, it
1508 // can be a sequence of decimal digits.
1509 //
1510 // Update the man variable during the walk. It's OK if man overflows now.
1511 // We'll detect that later.
1512 uint64_t man;
1513 if (*p == '0') {
1514 man = 0;
1515 p++;
1516 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1517 goto fallback;
1518 }
1519 } else if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1520 man = ((uint8_t)(*p - '0'));
1521 p++;
1522 for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) {
1523 man = (10 * man) + ((uint8_t)(*p - '0'));
1524 }
1525 } else {
1526 goto fallback;
1527 }
1528
1529 // Walk the "d"s after the optional decimal separator ('.' or ','),
1530 // updating the man and exp10 variables.
1531 int32_t exp10 = 0;
Nigel Taoe0c5de92020-07-11 11:48:17 +10001532 if (*p ==
1533 ((options & WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1534 ? ','
1535 : '.')) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001536 p++;
1537 const uint8_t* first_after_separator_ptr = p;
1538 if (!wuffs_base__private_implementation__is_decimal_digit(*p)) {
1539 goto fallback;
1540 }
1541 man = (10 * man) + ((uint8_t)(*p - '0'));
1542 p++;
1543 for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) {
1544 man = (10 * man) + ((uint8_t)(*p - '0'));
1545 }
1546 exp10 = ((int32_t)(first_after_separator_ptr - p));
1547 }
1548
1549 // Count the number of digits:
1550 // - for an input of "314159", digit_count is 6.
1551 // - for an input of "3.14159", digit_count is 7.
1552 //
1553 // This is off-by-one if there is a decimal separator. That's OK for now.
1554 // We'll correct for that later. The "script/process-json-numbers.c with
1555 // -p" benchmark is noticably slower if we try to correct for that now.
1556 uint32_t digit_count = (uint32_t)(p - start_of_digits_ptr);
1557
1558 // Update exp10 for the optional exponent, starting with 'E' or 'e'.
1559 if ((*p | 0x20) == 'e') {
1560 p++;
1561 int32_t exp_sign = +1;
1562 if (*p == '-') {
1563 p++;
1564 exp_sign = -1;
1565 } else if (*p == '+') {
1566 p++;
1567 }
1568 if (!wuffs_base__private_implementation__is_decimal_digit(*p)) {
1569 goto fallback;
1570 }
1571 int32_t exp_num = ((uint8_t)(*p - '0'));
1572 p++;
1573 // The rest of the exp_num walking has a peculiar control flow but, once
1574 // again, the "script/process-json-numbers.c with -p" benchmark is
1575 // sensitive to alternative formulations.
1576 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1577 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1578 p++;
1579 }
1580 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1581 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1582 p++;
1583 }
1584 while (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1585 if (exp_num > 0x1000000) {
1586 goto fallback;
1587 }
1588 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1589 p++;
1590 }
1591 exp10 += exp_sign * exp_num;
1592 }
1593
1594 // The Wuffs API is that the original slice has no trailing data. It also
1595 // allows underscores, which we don't catch here but the fallback should.
1596 if (p != &z[s.len]) {
1597 goto fallback;
1598 }
1599
1600 // Check that the uint64_t typed man variable has not overflowed, based on
1601 // digit_count.
1602 //
1603 // For reference:
1604 // - (1 << 63) is 9223372036854775808, which has 19 decimal digits.
1605 // - (1 << 64) is 18446744073709551616, which has 20 decimal digits.
1606 // - 19 nines, 9999999999999999999, is 0x8AC7230489E7FFFF, which has 64
1607 // bits and 16 hexadecimal digits.
1608 // - 20 nines, 99999999999999999999, is 0x56BC75E2D630FFFFF, which has 67
1609 // bits and 17 hexadecimal digits.
1610 if (digit_count > 19) {
1611 // Even if we have more than 19 pseudo-digits, it's not yet definitely an
1612 // overflow. Recall that digit_count might be off-by-one (too large) if
1613 // there's a decimal separator. It will also over-report the number of
1614 // meaningful digits if the input looks something like "0.000dddExxx".
1615 //
1616 // We adjust by the number of leading '0's and '.'s and re-compare to 19.
1617 // Once again, technically, we could skip ','s too, but that perturbs the
1618 // "script/process-json-numbers.c with -p" benchmark.
1619 const uint8_t* q = start_of_digits_ptr;
1620 for (; (*q == '0') || (*q == '.'); q++) {
1621 }
1622 digit_count -= (uint32_t)(q - start_of_digits_ptr);
1623 if (digit_count > 19) {
1624 goto fallback;
1625 }
1626 }
1627
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001628 // The wuffs_base__private_implementation__parse_number_f64_eisel_lemire
Nigel Tao8b45db02020-09-15 21:50:32 +10001629 // preconditions include that exp10 is in the range [-307 ..= 288].
1630 if ((exp10 < -307) || (288 < exp10)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001631 goto fallback;
1632 }
1633
Nigel Tao9f22b5e2020-09-11 09:10:08 +10001634 // If both man and (10 ** exp10) are exactly representable by a double, we
1635 // don't need to run the Eisel-Lemire algorithm.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001636 if ((-22 <= exp10) && (exp10 <= 22) && ((man >> 53) == 0)) {
1637 double d = (double)man;
1638 if (exp10 >= 0) {
1639 d *= wuffs_base__private_implementation__f64_powers_of_10[+exp10];
1640 } else {
1641 d /= wuffs_base__private_implementation__f64_powers_of_10[-exp10];
1642 }
1643 wuffs_base__result_f64 ret;
1644 ret.status.repr = NULL;
1645 ret.value = negative ? -d : +d;
1646 return ret;
1647 }
1648
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001649 // The wuffs_base__private_implementation__parse_number_f64_eisel_lemire
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001650 // preconditions include that man is non-zero. Parsing "0" should be caught
Nigel Tao9f22b5e2020-09-11 09:10:08 +10001651 // by the "If both man and (10 ** exp10)" above, but "0e99" might not.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001652 if (man == 0) {
1653 goto fallback;
1654 }
1655
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001656 // Our man and exp10 are in range. Run the Eisel-Lemire algorithm.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001657 int64_t r =
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001658 wuffs_base__private_implementation__parse_number_f64_eisel_lemire(
1659 man, exp10);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001660 if (r < 0) {
1661 goto fallback;
1662 }
1663 wuffs_base__result_f64 ret;
1664 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001665 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001666 ((uint64_t)r) | (((uint64_t)negative) << 63));
1667 return ret;
1668 } while (0);
1669
1670fallback:
1671 do {
1672 wuffs_base__private_implementation__high_prec_dec h;
1673 wuffs_base__status status =
Nigel Taoe0c5de92020-07-11 11:48:17 +10001674 wuffs_base__private_implementation__high_prec_dec__parse(&h, s,
1675 options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001676 if (status.repr) {
Nigel Taoe0c5de92020-07-11 11:48:17 +10001677 return wuffs_base__private_implementation__parse_number_f64_special(
Nigel Tao4d61a052020-07-11 12:34:40 +10001678 s, options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001679 }
Nigel Tao4d61a052020-07-11 12:34:40 +10001680 return wuffs_base__private_implementation__high_prec_dec__to_f64(&h,
1681 options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001682 } while (0);
1683}
1684
1685// --------
1686
1687static inline size_t //
1688wuffs_base__private_implementation__render_inf(wuffs_base__slice_u8 dst,
1689 bool neg,
1690 uint32_t options) {
1691 if (neg) {
1692 if (dst.len < 4) {
1693 return 0;
1694 }
Nigel Taoa1c22ca2021-01-17 22:22:49 +11001695 wuffs_base__poke_u32le__no_bounds_check(dst.ptr, 0x666E492D); // '-Inf'le.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001696 return 4;
1697 }
1698
1699 if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1700 if (dst.len < 4) {
1701 return 0;
1702 }
Nigel Taoa1c22ca2021-01-17 22:22:49 +11001703 wuffs_base__poke_u32le__no_bounds_check(dst.ptr, 0x666E492B); // '+Inf'le.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001704 return 4;
1705 }
1706
1707 if (dst.len < 3) {
1708 return 0;
1709 }
Nigel Taoa1c22ca2021-01-17 22:22:49 +11001710 wuffs_base__poke_u24le__no_bounds_check(dst.ptr, 0x666E49); // 'Inf'le.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001711 return 3;
1712}
1713
1714static inline size_t //
1715wuffs_base__private_implementation__render_nan(wuffs_base__slice_u8 dst) {
1716 if (dst.len < 3) {
1717 return 0;
1718 }
Nigel Taoa1c22ca2021-01-17 22:22:49 +11001719 wuffs_base__poke_u24le__no_bounds_check(dst.ptr, 0x4E614E); // 'NaN'le.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001720 return 3;
1721}
1722
1723static size_t //
1724wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1725 wuffs_base__slice_u8 dst,
1726 wuffs_base__private_implementation__high_prec_dec* h,
1727 uint32_t precision,
1728 uint32_t options) {
1729 size_t n = (h->negative ||
1730 (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN))
1731 ? 1
1732 : 0;
1733 if (h->decimal_point <= 0) {
1734 n += 1;
1735 } else {
1736 n += (size_t)(h->decimal_point);
1737 }
1738 if (precision > 0) {
1739 n += precision + 1; // +1 for the '.'.
1740 }
1741
1742 // Don't modify dst if the formatted number won't fit.
1743 if (n > dst.len) {
1744 return 0;
1745 }
1746
1747 // Align-left or align-right.
1748 uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT)
1749 ? &dst.ptr[dst.len - n]
1750 : &dst.ptr[0];
1751
1752 // Leading "±".
1753 if (h->negative) {
1754 *ptr++ = '-';
1755 } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1756 *ptr++ = '+';
1757 }
1758
1759 // Integral digits.
1760 if (h->decimal_point <= 0) {
1761 *ptr++ = '0';
1762 } else {
1763 uint32_t m =
1764 wuffs_base__u32__min(h->num_digits, (uint32_t)(h->decimal_point));
1765 uint32_t i = 0;
1766 for (; i < m; i++) {
1767 *ptr++ = (uint8_t)('0' | h->digits[i]);
1768 }
1769 for (; i < (uint32_t)(h->decimal_point); i++) {
1770 *ptr++ = '0';
1771 }
1772 }
1773
1774 // Separator and then fractional digits.
1775 if (precision > 0) {
1776 *ptr++ =
1777 (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1778 ? ','
1779 : '.';
1780 uint32_t i = 0;
1781 for (; i < precision; i++) {
1782 uint32_t j = ((uint32_t)(h->decimal_point)) + i;
1783 *ptr++ = (uint8_t)('0' | ((j < h->num_digits) ? h->digits[j] : 0));
1784 }
1785 }
1786
1787 return n;
1788}
1789
1790static size_t //
1791wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1792 wuffs_base__slice_u8 dst,
1793 wuffs_base__private_implementation__high_prec_dec* h,
1794 uint32_t precision,
1795 uint32_t options) {
1796 int32_t exp = 0;
1797 if (h->num_digits > 0) {
1798 exp = h->decimal_point - 1;
1799 }
1800 bool negative_exp = exp < 0;
1801 if (negative_exp) {
1802 exp = -exp;
1803 }
1804
1805 size_t n = (h->negative ||
1806 (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN))
1807 ? 4
1808 : 3; // Mininum 3 bytes: first digit and then "e±".
1809 if (precision > 0) {
1810 n += precision + 1; // +1 for the '.'.
1811 }
1812 n += (exp < 100) ? 2 : 3;
1813
1814 // Don't modify dst if the formatted number won't fit.
1815 if (n > dst.len) {
1816 return 0;
1817 }
1818
1819 // Align-left or align-right.
1820 uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT)
1821 ? &dst.ptr[dst.len - n]
1822 : &dst.ptr[0];
1823
1824 // Leading "±".
1825 if (h->negative) {
1826 *ptr++ = '-';
1827 } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1828 *ptr++ = '+';
1829 }
1830
1831 // Integral digit.
1832 if (h->num_digits > 0) {
1833 *ptr++ = (uint8_t)('0' | h->digits[0]);
1834 } else {
1835 *ptr++ = '0';
1836 }
1837
1838 // Separator and then fractional digits.
1839 if (precision > 0) {
1840 *ptr++ =
1841 (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1842 ? ','
1843 : '.';
1844 uint32_t i = 1;
1845 uint32_t j = wuffs_base__u32__min(h->num_digits, precision + 1);
1846 for (; i < j; i++) {
1847 *ptr++ = (uint8_t)('0' | h->digits[i]);
1848 }
1849 for (; i <= precision; i++) {
1850 *ptr++ = '0';
1851 }
1852 }
1853
1854 // Exponent: "e±" and then 2 or 3 digits.
1855 *ptr++ = 'e';
1856 *ptr++ = negative_exp ? '-' : '+';
1857 if (exp < 10) {
1858 *ptr++ = '0';
1859 *ptr++ = (uint8_t)('0' | exp);
1860 } else if (exp < 100) {
1861 *ptr++ = (uint8_t)('0' | (exp / 10));
1862 *ptr++ = (uint8_t)('0' | (exp % 10));
1863 } else {
1864 int32_t e = exp / 100;
1865 exp -= e * 100;
1866 *ptr++ = (uint8_t)('0' | e);
1867 *ptr++ = (uint8_t)('0' | (exp / 10));
1868 *ptr++ = (uint8_t)('0' | (exp % 10));
1869 }
1870
1871 return n;
1872}
1873
1874WUFFS_BASE__MAYBE_STATIC size_t //
1875wuffs_base__render_number_f64(wuffs_base__slice_u8 dst,
1876 double x,
1877 uint32_t precision,
1878 uint32_t options) {
1879 // Decompose x (64 bits) into negativity (1 bit), base-2 exponent (11 bits
1880 // with a -1023 bias) and mantissa (52 bits).
Nigel Tao4d449dc2020-07-12 11:00:47 +10001881 uint64_t bits = wuffs_base__ieee_754_bit_representation__from_f64_to_u64(x);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001882 bool neg = (bits >> 63) != 0;
1883 int32_t exp2 = ((int32_t)(bits >> 52)) & 0x7FF;
1884 uint64_t man = bits & 0x000FFFFFFFFFFFFFul;
1885
1886 // Apply the exponent bias and set the implicit top bit of the mantissa,
1887 // unless x is subnormal. Also take care of Inf and NaN.
1888 if (exp2 == 0x7FF) {
1889 if (man != 0) {
1890 return wuffs_base__private_implementation__render_nan(dst);
1891 }
1892 return wuffs_base__private_implementation__render_inf(dst, neg, options);
1893 } else if (exp2 == 0) {
1894 exp2 = -1022;
1895 } else {
1896 exp2 -= 1023;
1897 man |= 0x0010000000000000ul;
1898 }
1899
1900 // Ensure that precision isn't too large.
1901 if (precision > 4095) {
1902 precision = 4095;
1903 }
1904
1905 // Convert from the (neg, exp2, man) tuple to an HPD.
1906 wuffs_base__private_implementation__high_prec_dec h;
1907 wuffs_base__private_implementation__high_prec_dec__assign(&h, man, neg);
1908 if (h.num_digits > 0) {
1909 wuffs_base__private_implementation__high_prec_dec__lshift(
1910 &h, exp2 - 52); // 52 mantissa bits.
1911 }
1912
1913 // Handle the "%e" and "%f" formats.
1914 switch (options & (WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT |
1915 WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT)) {
1916 case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT: // The "%"f" format.
1917 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1918 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1919 &h, exp2, man);
1920 int32_t p = ((int32_t)(h.num_digits)) - h.decimal_point;
1921 precision = ((uint32_t)(wuffs_base__i32__max(0, p)));
1922 } else {
1923 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1924 &h, ((int32_t)precision) + h.decimal_point);
1925 }
1926 return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1927 dst, &h, precision, options);
1928
1929 case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT: // The "%e" format.
1930 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1931 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1932 &h, exp2, man);
1933 precision = (h.num_digits > 0) ? (h.num_digits - 1) : 0;
1934 } else {
1935 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1936 &h, ((int32_t)precision) + 1);
1937 }
1938 return wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1939 dst, &h, precision, options);
1940 }
1941
1942 // We have the "%g" format and so precision means the number of significant
1943 // digits, not the number of digits after the decimal separator. Perform
1944 // rounding and determine whether to use "%e" or "%f".
1945 int32_t e_threshold = 0;
1946 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1947 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1948 &h, exp2, man);
1949 precision = h.num_digits;
1950 e_threshold = 6;
1951 } else {
1952 if (precision == 0) {
1953 precision = 1;
1954 }
1955 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1956 &h, ((int32_t)precision));
1957 e_threshold = ((int32_t)precision);
1958 int32_t nd = ((int32_t)(h.num_digits));
1959 if ((e_threshold > nd) && (nd >= h.decimal_point)) {
1960 e_threshold = nd;
1961 }
1962 }
1963
1964 // Use the "%e" format if the exponent is large.
1965 int32_t e = h.decimal_point - 1;
1966 if ((e < -4) || (e_threshold <= e)) {
1967 uint32_t p = wuffs_base__u32__min(precision, h.num_digits);
1968 return wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1969 dst, &h, (p > 0) ? (p - 1) : 0, options);
1970 }
1971
1972 // Use the "%f" format otherwise.
1973 int32_t p = ((int32_t)precision);
1974 if (p > h.decimal_point) {
1975 p = ((int32_t)(h.num_digits));
1976 }
1977 precision = ((uint32_t)(wuffs_base__i32__max(0, p - h.decimal_point)));
1978 return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1979 dst, &h, precision, options);
1980}