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