blob: a07e23281ce91f89365c62f20756a30aff847d56 [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.
200typedef struct {
201 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/
978//
979// Preconditions:
980// - man is non-zero.
981// - exp10 is in the range -326 ..= 310, the same range of the
982// wuffs_base__private_implementation__powers_of_10 array.
983static int64_t //
Nigel Taoc4fa8e22020-07-18 17:35:13 +1000984wuffs_base__private_implementation__parse_number_f64_eisel_lemire(
985 uint64_t man,
986 int32_t exp10) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000987 // Look up the (possibly truncated) base-2 representation of (10 ** exp10).
988 // The look-up table was constructed so that it is already normalized: the
989 // table entry's mantissa's MSB (most significant bit) is on.
990 const uint32_t* po10 =
991 &wuffs_base__private_implementation__powers_of_10[5 * (exp10 + 326)];
992
993 // Normalize the man argument. The (man != 0) precondition means that a
994 // non-zero bit exists.
995 uint32_t clz = wuffs_base__count_leading_zeroes_u64(man);
996 man <<= clz;
997
998 // Calculate the return value's base-2 exponent. We might tweak it by ±1
999 // later, but its initial value comes from the look-up table and clz.
1000 uint64_t ret_exp2 = ((uint64_t)po10[4]) - ((uint64_t)clz);
1001
1002 // Multiply the two mantissas. Normalization means that both mantissas are at
1003 // least (1<<63), so the 128-bit product must be at least (1<<126). The high
Nigel Tao74d4af62020-07-10 11:27:17 +10001004 // 64 bits of the product, x_hi, must therefore be at least (1<<62).
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001005 //
Nigel Tao74d4af62020-07-10 11:27:17 +10001006 // As a consequence, x_hi has either 0 or 1 leading zeroes. Shifting x_hi
1007 // right by either 9 or 10 bits (depending on x_hi's MSB) will therefore
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001008 // leave the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on.
Nigel Taoe722bb02020-07-10 20:33:09 +10001009#if defined(__SIZEOF_INT128__)
Nigel Taobff01d72020-07-10 12:15:09 +10001010 // See commit 18449ad75d582dd015c236abc85a16f333b796f3 "Optimize 128-bit muls
1011 // in parse_number_f64_eisel" for benchmark numbers.
Nigel Tao18449ad2020-07-10 11:48:52 +10001012 __uint128_t x =
1013 ((__uint128_t)man) * (((uint64_t)po10[2]) | (((uint64_t)po10[3]) << 32));
1014 uint64_t x_hi = ((uint64_t)(x >> 64));
1015 uint64_t x_lo = ((uint64_t)(x));
1016#else
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001017 wuffs_base__multiply_u64__output x = wuffs_base__multiply_u64(
1018 man, ((uint64_t)po10[2]) | (((uint64_t)po10[3]) << 32));
Nigel Tao74d4af62020-07-10 11:27:17 +10001019 uint64_t x_hi = x.hi;
1020 uint64_t x_lo = x.lo;
Nigel Tao18449ad2020-07-10 11:48:52 +10001021#endif
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001022
1023 // Before we shift right by at least 9 bits, recall that the look-up table
1024 // entry was possibly truncated. We have so far only calculated a lower bound
1025 // for the product (man * e), where e is (10 ** exp10). The upper bound would
1026 // add a further (man * 1) to the 128-bit product, which overflows the lower
Nigel Tao74d4af62020-07-10 11:27:17 +10001027 // 64-bit limb if ((x_lo + man) < man).
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001028 //
Nigel Tao74d4af62020-07-10 11:27:17 +10001029 // If overflow occurs, that adds 1 to x_hi. Since we're about to shift right
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001030 // by at least 9 bits, that carried 1 can be ignored unless the higher 64-bit
1031 // limb's low 9 bits are all on.
Nigel Tao74d4af62020-07-10 11:27:17 +10001032 if (((x_hi & 0x1FF) == 0x1FF) && ((x_lo + man) < man)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001033 // Refine our calculation of (man * e). Before, our approximation of e used
1034 // a "low resolution" 64-bit mantissa. Now use a "high resolution" 128-bit
1035 // mantissa. We've already calculated x = (man * bits_0_to_63_incl_of_e).
1036 // Now calculate y = (man * bits_64_to_127_incl_of_e).
Nigel Taoe722bb02020-07-10 20:33:09 +10001037#if defined(__SIZEOF_INT128__)
Nigel Taobff01d72020-07-10 12:15:09 +10001038 // See commit 18449ad75d582dd015c236abc85a16f333b796f3 "Optimize 128-bit
1039 // muls in parse_number_f64_eisel" for benchmark numbers.
Nigel Tao18449ad2020-07-10 11:48:52 +10001040 __uint128_t y = ((__uint128_t)man) *
1041 (((uint64_t)po10[0]) | (((uint64_t)po10[1]) << 32));
1042 uint64_t y_hi = ((uint64_t)(y >> 64));
1043 uint64_t y_lo = ((uint64_t)(y));
1044#else
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001045 wuffs_base__multiply_u64__output y = wuffs_base__multiply_u64(
1046 man, ((uint64_t)po10[0]) | (((uint64_t)po10[1]) << 32));
Nigel Tao74d4af62020-07-10 11:27:17 +10001047 uint64_t y_hi = y.hi;
1048 uint64_t y_lo = y.lo;
Nigel Tao18449ad2020-07-10 11:48:52 +10001049#endif
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001050
1051 // Merge the 128-bit x and 128-bit y, which overlap by 64 bits, to
1052 // calculate the 192-bit product of the 64-bit man by the 128-bit e.
1053 // As we exit this if-block, we only care about the high 128 bits
1054 // (merged_hi and merged_lo) of that 192-bit product.
Nigel Tao74d4af62020-07-10 11:27:17 +10001055 uint64_t merged_hi = x_hi;
1056 uint64_t merged_lo = x_lo + y_hi;
1057 if (merged_lo < x_lo) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001058 merged_hi++; // Carry the overflow bit.
1059 }
1060
1061 // The "high resolution" approximation of e is still a lower bound. Once
1062 // again, see if the upper bound is large enough to produce a different
1063 // result. This time, if it does, give up instead of reaching for an even
1064 // more precise approximation to e.
1065 //
1066 // This three-part check is similar to the two-part check that guarded the
1067 // if block that we're now in, but it has an extra term for the middle 64
1068 // bits (checking that adding 1 to merged_lo would overflow).
1069 if (((merged_hi & 0x1FF) == 0x1FF) && ((merged_lo + 1) == 0) &&
Nigel Tao74d4af62020-07-10 11:27:17 +10001070 (y_lo + man < man)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001071 return -1;
1072 }
1073
1074 // Replace the 128-bit x with merged.
Nigel Tao74d4af62020-07-10 11:27:17 +10001075 x_hi = merged_hi;
1076 x_lo = merged_lo;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001077 }
1078
Nigel Tao74d4af62020-07-10 11:27:17 +10001079 // As mentioned above, shifting x_hi right by either 9 or 10 bits will leave
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001080 // the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on. If the
1081 // MSB (before shifting) was on, adjust ret_exp2 for the larger shift.
1082 //
1083 // Having bit 53 on (and higher bits off) means that ret_mantissa is a 54-bit
1084 // number.
Nigel Tao74d4af62020-07-10 11:27:17 +10001085 uint64_t msb = x_hi >> 63;
1086 uint64_t ret_mantissa = x_hi >> (msb + 9);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001087 ret_exp2 -= 1 ^ msb;
1088
1089 // IEEE 754 rounds to-nearest with ties rounded to-even. Rounding to-even can
1090 // be tricky. If we're half-way between two exactly representable numbers
1091 // (x's low 73 bits are zero and the next 2 bits that matter are "01"), give
1092 // up instead of trying to pick the winner.
1093 //
1094 // Technically, we could tighten the condition by changing "73" to "73 or 74,
1095 // depending on msb", but a flat "73" is simpler.
Nigel Tao74d4af62020-07-10 11:27:17 +10001096 if ((x_lo == 0) && ((x_hi & 0x1FF) == 0) && ((ret_mantissa & 3) == 1)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001097 return -1;
1098 }
1099
1100 // If we're not halfway then it's rounding to-nearest. Starting with a 54-bit
1101 // number, carry the lowest bit (bit 0) up if it's on. Regardless of whether
1102 // it was on or off, shifting right by one then produces a 53-bit number. If
1103 // carrying up overflowed, shift again.
1104 ret_mantissa += ret_mantissa & 1;
1105 ret_mantissa >>= 1;
1106 if ((ret_mantissa >> 53) > 0) {
1107 ret_mantissa >>= 1;
1108 ret_exp2++;
1109 }
1110
1111 // Starting with a 53-bit number, IEEE 754 double-precision normal numbers
1112 // have an implicit mantissa bit. Mask that away and keep the low 52 bits.
1113 ret_mantissa &= 0x000FFFFFFFFFFFFF;
1114
1115 // IEEE 754 double-precision floating point has 11 exponent bits. All off (0)
1116 // means subnormal numbers. All on (2047) means infinity or NaN.
1117 if ((ret_exp2 <= 0) || (2047 <= ret_exp2)) {
1118 return -1;
1119 }
1120
1121 // Pack the bits and return.
1122 return ((int64_t)(ret_mantissa | (ret_exp2 << 52)));
1123}
1124
1125// --------
1126
1127static wuffs_base__result_f64 //
Nigel Taoe0c5de92020-07-11 11:48:17 +10001128wuffs_base__private_implementation__parse_number_f64_special(
1129 wuffs_base__slice_u8 s,
Nigel Tao4d61a052020-07-11 12:34:40 +10001130 uint32_t options) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001131 do {
Nigel Tao4d61a052020-07-11 12:34:40 +10001132 if (options & WUFFS_BASE__PARSE_NUMBER_FXX__REJECT_INF_AND_NAN) {
1133 goto fail;
1134 }
1135
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001136 uint8_t* p = s.ptr;
1137 uint8_t* q = s.ptr + s.len;
1138
1139 for (; (p < q) && (*p == '_'); p++) {
1140 }
1141 if (p >= q) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001142 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001143 }
1144
1145 // Parse sign.
1146 bool negative = false;
1147 do {
1148 if (*p == '+') {
1149 p++;
1150 } else if (*p == '-') {
1151 negative = true;
1152 p++;
1153 } else {
1154 break;
1155 }
1156 for (; (p < q) && (*p == '_'); p++) {
1157 }
1158 } while (0);
1159 if (p >= q) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001160 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001161 }
1162
1163 bool nan = false;
1164 switch (p[0]) {
1165 case 'I':
1166 case 'i':
1167 if (((q - p) < 3) || //
1168 ((p[1] != 'N') && (p[1] != 'n')) || //
1169 ((p[2] != 'F') && (p[2] != 'f'))) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001170 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001171 }
1172 p += 3;
1173
1174 if ((p >= q) || (*p == '_')) {
1175 break;
1176 } else if (((q - p) < 5) || //
1177 ((p[0] != 'I') && (p[0] != 'i')) || //
1178 ((p[1] != 'N') && (p[1] != 'n')) || //
1179 ((p[2] != 'I') && (p[2] != 'i')) || //
1180 ((p[3] != 'T') && (p[3] != 't')) || //
1181 ((p[4] != 'Y') && (p[4] != 'y'))) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001182 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001183 }
1184 p += 5;
1185
1186 if ((p >= q) || (*p == '_')) {
1187 break;
1188 }
Nigel Tao4d61a052020-07-11 12:34:40 +10001189 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001190
1191 case 'N':
1192 case 'n':
1193 if (((q - p) < 3) || //
1194 ((p[1] != 'A') && (p[1] != 'a')) || //
1195 ((p[2] != 'N') && (p[2] != 'n'))) {
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 nan = true;
1202 break;
1203 }
Nigel Tao4d61a052020-07-11 12:34:40 +10001204 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001205
1206 default:
Nigel Tao4d61a052020-07-11 12:34:40 +10001207 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001208 }
1209
1210 // Finish.
1211 for (; (p < q) && (*p == '_'); p++) {
1212 }
1213 if (p != q) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001214 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001215 }
1216 wuffs_base__result_f64 ret;
1217 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001218 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001219 (nan ? 0x7FFFFFFFFFFFFFFF : 0x7FF0000000000000) |
1220 (negative ? 0x8000000000000000 : 0));
1221 return ret;
1222 } while (0);
1223
Nigel Tao4d61a052020-07-11 12:34:40 +10001224fail:
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001225 do {
1226 wuffs_base__result_f64 ret;
Nigel Tao4d61a052020-07-11 12:34:40 +10001227 ret.status.repr = wuffs_base__error__bad_argument;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001228 ret.value = 0;
1229 return ret;
1230 } while (0);
1231}
1232
1233WUFFS_BASE__MAYBE_STATIC wuffs_base__result_f64 //
Nigel Taoe0c5de92020-07-11 11:48:17 +10001234wuffs_base__private_implementation__high_prec_dec__to_f64(
Nigel Tao4d61a052020-07-11 12:34:40 +10001235 wuffs_base__private_implementation__high_prec_dec* h,
1236 uint32_t options) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001237 do {
1238 // powers converts decimal powers of 10 to binary powers of 2. For example,
1239 // (10000 >> 13) is 1. It stops before the elements exceed 60, also known
1240 // as WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL.
1241 static const uint32_t num_powers = 19;
1242 static const uint8_t powers[19] = {
1243 0, 3, 6, 9, 13, 16, 19, 23, 26, 29, //
1244 33, 36, 39, 43, 46, 49, 53, 56, 59, //
1245 };
1246
1247 // Handle zero and obvious extremes. The largest and smallest positive
1248 // finite f64 values are approximately 1.8e+308 and 4.9e-324.
1249 if ((h->num_digits == 0) || (h->decimal_point < -326)) {
1250 goto zero;
1251 } else if (h->decimal_point > 310) {
1252 goto infinity;
1253 }
1254
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001255 // Try the fast Eisel-Lemire algorithm again. Calculating the (man, exp10)
1256 // pair from the high_prec_dec h is more correct but slower than the
1257 // approach taken in wuffs_base__parse_number_f64. The latter is optimized
1258 // for the common cases (e.g. assuming no underscores or a leading '+'
1259 // sign) rather than the full set of cases allowed by the Wuffs API.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001260 if (h->num_digits <= 19) {
1261 uint64_t man = 0;
1262 uint32_t i;
1263 for (i = 0; i < h->num_digits; i++) {
1264 man = (10 * man) + h->digits[i];
1265 }
1266 int32_t exp10 = h->decimal_point - ((int32_t)(h->num_digits));
1267 if ((man != 0) && (-326 <= exp10) && (exp10 <= 310)) {
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001268 int64_t r =
1269 wuffs_base__private_implementation__parse_number_f64_eisel_lemire(
1270 man, exp10);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001271 if (r >= 0) {
1272 wuffs_base__result_f64 ret;
1273 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001274 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001275 ((uint64_t)r) | (((uint64_t)(h->negative)) << 63));
1276 return ret;
1277 }
1278 }
1279 }
1280
1281 // Scale by powers of 2 until we're in the range [½ .. 1], which gives us
1282 // our exponent (in base-2). First we shift right, possibly a little too
1283 // far, ending with a value certainly below 1 and possibly below ½...
1284 const int32_t f64_bias = -1023;
1285 int32_t exp2 = 0;
1286 while (h->decimal_point > 0) {
1287 uint32_t n = (uint32_t)(+h->decimal_point);
1288 uint32_t shift =
1289 (n < num_powers)
1290 ? powers[n]
1291 : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1292
1293 wuffs_base__private_implementation__high_prec_dec__small_rshift(h, shift);
1294 if (h->decimal_point <
1295 -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
1296 goto zero;
1297 }
1298 exp2 += (int32_t)shift;
1299 }
1300 // ...then we shift left, putting us in [½ .. 1].
1301 while (h->decimal_point <= 0) {
1302 uint32_t shift;
1303 if (h->decimal_point == 0) {
1304 if (h->digits[0] >= 5) {
1305 break;
1306 }
1307 shift = (h->digits[0] <= 2) ? 2 : 1;
1308 } else {
1309 uint32_t n = (uint32_t)(-h->decimal_point);
1310 shift = (n < num_powers)
1311 ? powers[n]
1312 : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1313 }
1314
1315 wuffs_base__private_implementation__high_prec_dec__small_lshift(h, shift);
1316 if (h->decimal_point >
1317 +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
1318 goto infinity;
1319 }
1320 exp2 -= (int32_t)shift;
1321 }
1322
1323 // We're in the range [½ .. 1] but f64 uses [1 .. 2].
1324 exp2--;
1325
1326 // The minimum normal exponent is (f64_bias + 1).
1327 while ((f64_bias + 1) > exp2) {
1328 uint32_t n = (uint32_t)((f64_bias + 1) - exp2);
1329 if (n > WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) {
1330 n = WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1331 }
1332 wuffs_base__private_implementation__high_prec_dec__small_rshift(h, n);
1333 exp2 += (int32_t)n;
1334 }
1335
1336 // Check for overflow.
1337 if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1.
1338 goto infinity;
1339 }
1340
1341 // Extract 53 bits for the mantissa (in base-2).
1342 wuffs_base__private_implementation__high_prec_dec__small_lshift(h, 53);
1343 uint64_t man2 =
1344 wuffs_base__private_implementation__high_prec_dec__rounded_integer(h);
1345
1346 // Rounding might have added one bit. If so, shift and re-check overflow.
1347 if ((man2 >> 53) != 0) {
1348 man2 >>= 1;
1349 exp2++;
1350 if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1.
1351 goto infinity;
1352 }
1353 }
1354
1355 // Handle subnormal numbers.
1356 if ((man2 >> 52) == 0) {
1357 exp2 = f64_bias;
1358 }
1359
1360 // Pack the bits and return.
1361 uint64_t exp2_bits =
1362 (uint64_t)((exp2 - f64_bias) & 0x07FF); // (1 << 11) - 1.
1363 uint64_t bits = (man2 & 0x000FFFFFFFFFFFFF) | // (1 << 52) - 1.
1364 (exp2_bits << 52) | //
1365 (h->negative ? 0x8000000000000000 : 0); // (1 << 63).
1366
1367 wuffs_base__result_f64 ret;
1368 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001369 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001370 return ret;
1371 } while (0);
1372
1373zero:
1374 do {
1375 uint64_t bits = h->negative ? 0x8000000000000000 : 0;
1376
1377 wuffs_base__result_f64 ret;
1378 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001379 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001380 return ret;
1381 } while (0);
1382
1383infinity:
1384 do {
Nigel Tao4d61a052020-07-11 12:34:40 +10001385 if (options & WUFFS_BASE__PARSE_NUMBER_FXX__REJECT_INF_AND_NAN) {
1386 wuffs_base__result_f64 ret;
1387 ret.status.repr = wuffs_base__error__bad_argument;
1388 ret.value = 0;
1389 return ret;
1390 }
1391
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001392 uint64_t bits = h->negative ? 0xFFF0000000000000 : 0x7FF0000000000000;
1393
1394 wuffs_base__result_f64 ret;
1395 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001396 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001397 return ret;
1398 } while (0);
1399}
1400
1401static inline bool //
1402wuffs_base__private_implementation__is_decimal_digit(uint8_t c) {
1403 return ('0' <= c) && (c <= '9');
1404}
1405
1406WUFFS_BASE__MAYBE_STATIC wuffs_base__result_f64 //
1407wuffs_base__parse_number_f64(wuffs_base__slice_u8 s, uint32_t options) {
1408 // In practice, almost all "dd.ddddE±xxx" numbers can be represented
1409 // losslessly by a uint64_t mantissa "dddddd" and an int32_t base-10
1410 // exponent, adjusting "xxx" for the position (if present) of the decimal
1411 // separator '.' or ','.
1412 //
1413 // This (u64 man, i32 exp10) data structure is superficially similar to the
1414 // "Do It Yourself Floating Point" type from Loitsch (†), but the exponent
1415 // here is base-10, not base-2.
1416 //
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001417 // If s's number fits in a (man, exp10), parse that pair with the
1418 // Eisel-Lemire algorithm. If not, or if Eisel-Lemire fails, parsing s with
1419 // the fallback algorithm is slower but comprehensive.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001420 //
1421 // † "Printing Floating-Point Numbers Quickly and Accurately with Integers"
1422 // (https://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf).
1423 // Florian Loitsch is also the primary contributor to
1424 // https://github.com/google/double-conversion
1425 do {
1426 // Calculating that (man, exp10) pair needs to stay within s's bounds.
1427 // Provided that s isn't extremely long, work on a NUL-terminated copy of
1428 // s's contents. The NUL byte isn't a valid part of "±dd.ddddE±xxx".
1429 //
1430 // As the pointer p walks the contents, it's faster to repeatedly check "is
1431 // *p a valid digit" than "is p within bounds and *p a valid digit".
1432 if (s.len >= 256) {
1433 goto fallback;
1434 }
1435 uint8_t z[256];
1436 memcpy(&z[0], s.ptr, s.len);
1437 z[s.len] = 0;
1438 const uint8_t* p = &z[0];
1439
1440 // Look for a leading minus sign. Technically, we could also look for an
1441 // optional plus sign, but the "script/process-json-numbers.c with -p"
1442 // benchmark is noticably slower if we do. It's optional and, in practice,
1443 // usually absent. Let the fallback catch it.
1444 bool negative = (*p == '-');
1445 if (negative) {
1446 p++;
1447 }
1448
1449 // After walking "dd.dddd", comparing p later with p now will produce the
1450 // number of "d"s and "."s.
1451 const uint8_t* const start_of_digits_ptr = p;
1452
1453 // Walk the "d"s before a '.', 'E', NUL byte, etc. If it starts with '0',
1454 // it must be a single '0'. If it starts with a non-zero decimal digit, it
1455 // can be a sequence of decimal digits.
1456 //
1457 // Update the man variable during the walk. It's OK if man overflows now.
1458 // We'll detect that later.
1459 uint64_t man;
1460 if (*p == '0') {
1461 man = 0;
1462 p++;
1463 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1464 goto fallback;
1465 }
1466 } else if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1467 man = ((uint8_t)(*p - '0'));
1468 p++;
1469 for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) {
1470 man = (10 * man) + ((uint8_t)(*p - '0'));
1471 }
1472 } else {
1473 goto fallback;
1474 }
1475
1476 // Walk the "d"s after the optional decimal separator ('.' or ','),
1477 // updating the man and exp10 variables.
1478 int32_t exp10 = 0;
Nigel Taoe0c5de92020-07-11 11:48:17 +10001479 if (*p ==
1480 ((options & WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1481 ? ','
1482 : '.')) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001483 p++;
1484 const uint8_t* first_after_separator_ptr = p;
1485 if (!wuffs_base__private_implementation__is_decimal_digit(*p)) {
1486 goto fallback;
1487 }
1488 man = (10 * man) + ((uint8_t)(*p - '0'));
1489 p++;
1490 for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) {
1491 man = (10 * man) + ((uint8_t)(*p - '0'));
1492 }
1493 exp10 = ((int32_t)(first_after_separator_ptr - p));
1494 }
1495
1496 // Count the number of digits:
1497 // - for an input of "314159", digit_count is 6.
1498 // - for an input of "3.14159", digit_count is 7.
1499 //
1500 // This is off-by-one if there is a decimal separator. That's OK for now.
1501 // We'll correct for that later. The "script/process-json-numbers.c with
1502 // -p" benchmark is noticably slower if we try to correct for that now.
1503 uint32_t digit_count = (uint32_t)(p - start_of_digits_ptr);
1504
1505 // Update exp10 for the optional exponent, starting with 'E' or 'e'.
1506 if ((*p | 0x20) == 'e') {
1507 p++;
1508 int32_t exp_sign = +1;
1509 if (*p == '-') {
1510 p++;
1511 exp_sign = -1;
1512 } else if (*p == '+') {
1513 p++;
1514 }
1515 if (!wuffs_base__private_implementation__is_decimal_digit(*p)) {
1516 goto fallback;
1517 }
1518 int32_t exp_num = ((uint8_t)(*p - '0'));
1519 p++;
1520 // The rest of the exp_num walking has a peculiar control flow but, once
1521 // again, the "script/process-json-numbers.c with -p" benchmark is
1522 // sensitive to alternative formulations.
1523 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1524 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1525 p++;
1526 }
1527 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1528 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1529 p++;
1530 }
1531 while (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1532 if (exp_num > 0x1000000) {
1533 goto fallback;
1534 }
1535 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1536 p++;
1537 }
1538 exp10 += exp_sign * exp_num;
1539 }
1540
1541 // The Wuffs API is that the original slice has no trailing data. It also
1542 // allows underscores, which we don't catch here but the fallback should.
1543 if (p != &z[s.len]) {
1544 goto fallback;
1545 }
1546
1547 // Check that the uint64_t typed man variable has not overflowed, based on
1548 // digit_count.
1549 //
1550 // For reference:
1551 // - (1 << 63) is 9223372036854775808, which has 19 decimal digits.
1552 // - (1 << 64) is 18446744073709551616, which has 20 decimal digits.
1553 // - 19 nines, 9999999999999999999, is 0x8AC7230489E7FFFF, which has 64
1554 // bits and 16 hexadecimal digits.
1555 // - 20 nines, 99999999999999999999, is 0x56BC75E2D630FFFFF, which has 67
1556 // bits and 17 hexadecimal digits.
1557 if (digit_count > 19) {
1558 // Even if we have more than 19 pseudo-digits, it's not yet definitely an
1559 // overflow. Recall that digit_count might be off-by-one (too large) if
1560 // there's a decimal separator. It will also over-report the number of
1561 // meaningful digits if the input looks something like "0.000dddExxx".
1562 //
1563 // We adjust by the number of leading '0's and '.'s and re-compare to 19.
1564 // Once again, technically, we could skip ','s too, but that perturbs the
1565 // "script/process-json-numbers.c with -p" benchmark.
1566 const uint8_t* q = start_of_digits_ptr;
1567 for (; (*q == '0') || (*q == '.'); q++) {
1568 }
1569 digit_count -= (uint32_t)(q - start_of_digits_ptr);
1570 if (digit_count > 19) {
1571 goto fallback;
1572 }
1573 }
1574
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001575 // The wuffs_base__private_implementation__parse_number_f64_eisel_lemire
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001576 // preconditions include that exp10 is in the range -326 ..= 310.
1577 if ((exp10 < -326) || (310 < exp10)) {
1578 goto fallback;
1579 }
1580
1581 // If man and exp10 are small enough, all three of (man), (10 ** exp10) and
1582 // (man ** (10 ** exp10)) are exactly representable by a double. We don't
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001583 // need to run the Eisel-Lemire algorithm.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001584 if ((-22 <= exp10) && (exp10 <= 22) && ((man >> 53) == 0)) {
1585 double d = (double)man;
1586 if (exp10 >= 0) {
1587 d *= wuffs_base__private_implementation__f64_powers_of_10[+exp10];
1588 } else {
1589 d /= wuffs_base__private_implementation__f64_powers_of_10[-exp10];
1590 }
1591 wuffs_base__result_f64 ret;
1592 ret.status.repr = NULL;
1593 ret.value = negative ? -d : +d;
1594 return ret;
1595 }
1596
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001597 // The wuffs_base__private_implementation__parse_number_f64_eisel_lemire
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001598 // preconditions include that man is non-zero. Parsing "0" should be caught
1599 // by the "If man and exp10 are small enough" above, but "0e99" might not.
1600 if (man == 0) {
1601 goto fallback;
1602 }
1603
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001604 // Our man and exp10 are in range. Run the Eisel-Lemire algorithm.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001605 int64_t r =
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001606 wuffs_base__private_implementation__parse_number_f64_eisel_lemire(
1607 man, exp10);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001608 if (r < 0) {
1609 goto fallback;
1610 }
1611 wuffs_base__result_f64 ret;
1612 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001613 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001614 ((uint64_t)r) | (((uint64_t)negative) << 63));
1615 return ret;
1616 } while (0);
1617
1618fallback:
1619 do {
1620 wuffs_base__private_implementation__high_prec_dec h;
1621 wuffs_base__status status =
Nigel Taoe0c5de92020-07-11 11:48:17 +10001622 wuffs_base__private_implementation__high_prec_dec__parse(&h, s,
1623 options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001624 if (status.repr) {
Nigel Taoe0c5de92020-07-11 11:48:17 +10001625 return wuffs_base__private_implementation__parse_number_f64_special(
Nigel Tao4d61a052020-07-11 12:34:40 +10001626 s, options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001627 }
Nigel Tao4d61a052020-07-11 12:34:40 +10001628 return wuffs_base__private_implementation__high_prec_dec__to_f64(&h,
1629 options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001630 } while (0);
1631}
1632
1633// --------
1634
1635static inline size_t //
1636wuffs_base__private_implementation__render_inf(wuffs_base__slice_u8 dst,
1637 bool neg,
1638 uint32_t options) {
1639 if (neg) {
1640 if (dst.len < 4) {
1641 return 0;
1642 }
1643 wuffs_base__store_u32le__no_bounds_check(dst.ptr, 0x666E492D); // '-Inf'le.
1644 return 4;
1645 }
1646
1647 if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1648 if (dst.len < 4) {
1649 return 0;
1650 }
1651 wuffs_base__store_u32le__no_bounds_check(dst.ptr, 0x666E492B); // '+Inf'le.
1652 return 4;
1653 }
1654
1655 if (dst.len < 3) {
1656 return 0;
1657 }
1658 wuffs_base__store_u24le__no_bounds_check(dst.ptr, 0x666E49); // 'Inf'le.
1659 return 3;
1660}
1661
1662static inline size_t //
1663wuffs_base__private_implementation__render_nan(wuffs_base__slice_u8 dst) {
1664 if (dst.len < 3) {
1665 return 0;
1666 }
1667 wuffs_base__store_u24le__no_bounds_check(dst.ptr, 0x4E614E); // 'NaN'le.
1668 return 3;
1669}
1670
1671static size_t //
1672wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1673 wuffs_base__slice_u8 dst,
1674 wuffs_base__private_implementation__high_prec_dec* h,
1675 uint32_t precision,
1676 uint32_t options) {
1677 size_t n = (h->negative ||
1678 (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN))
1679 ? 1
1680 : 0;
1681 if (h->decimal_point <= 0) {
1682 n += 1;
1683 } else {
1684 n += (size_t)(h->decimal_point);
1685 }
1686 if (precision > 0) {
1687 n += precision + 1; // +1 for the '.'.
1688 }
1689
1690 // Don't modify dst if the formatted number won't fit.
1691 if (n > dst.len) {
1692 return 0;
1693 }
1694
1695 // Align-left or align-right.
1696 uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT)
1697 ? &dst.ptr[dst.len - n]
1698 : &dst.ptr[0];
1699
1700 // Leading "±".
1701 if (h->negative) {
1702 *ptr++ = '-';
1703 } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1704 *ptr++ = '+';
1705 }
1706
1707 // Integral digits.
1708 if (h->decimal_point <= 0) {
1709 *ptr++ = '0';
1710 } else {
1711 uint32_t m =
1712 wuffs_base__u32__min(h->num_digits, (uint32_t)(h->decimal_point));
1713 uint32_t i = 0;
1714 for (; i < m; i++) {
1715 *ptr++ = (uint8_t)('0' | h->digits[i]);
1716 }
1717 for (; i < (uint32_t)(h->decimal_point); i++) {
1718 *ptr++ = '0';
1719 }
1720 }
1721
1722 // Separator and then fractional digits.
1723 if (precision > 0) {
1724 *ptr++ =
1725 (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1726 ? ','
1727 : '.';
1728 uint32_t i = 0;
1729 for (; i < precision; i++) {
1730 uint32_t j = ((uint32_t)(h->decimal_point)) + i;
1731 *ptr++ = (uint8_t)('0' | ((j < h->num_digits) ? h->digits[j] : 0));
1732 }
1733 }
1734
1735 return n;
1736}
1737
1738static size_t //
1739wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1740 wuffs_base__slice_u8 dst,
1741 wuffs_base__private_implementation__high_prec_dec* h,
1742 uint32_t precision,
1743 uint32_t options) {
1744 int32_t exp = 0;
1745 if (h->num_digits > 0) {
1746 exp = h->decimal_point - 1;
1747 }
1748 bool negative_exp = exp < 0;
1749 if (negative_exp) {
1750 exp = -exp;
1751 }
1752
1753 size_t n = (h->negative ||
1754 (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN))
1755 ? 4
1756 : 3; // Mininum 3 bytes: first digit and then "e±".
1757 if (precision > 0) {
1758 n += precision + 1; // +1 for the '.'.
1759 }
1760 n += (exp < 100) ? 2 : 3;
1761
1762 // Don't modify dst if the formatted number won't fit.
1763 if (n > dst.len) {
1764 return 0;
1765 }
1766
1767 // Align-left or align-right.
1768 uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT)
1769 ? &dst.ptr[dst.len - n]
1770 : &dst.ptr[0];
1771
1772 // Leading "±".
1773 if (h->negative) {
1774 *ptr++ = '-';
1775 } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1776 *ptr++ = '+';
1777 }
1778
1779 // Integral digit.
1780 if (h->num_digits > 0) {
1781 *ptr++ = (uint8_t)('0' | h->digits[0]);
1782 } else {
1783 *ptr++ = '0';
1784 }
1785
1786 // Separator and then fractional digits.
1787 if (precision > 0) {
1788 *ptr++ =
1789 (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1790 ? ','
1791 : '.';
1792 uint32_t i = 1;
1793 uint32_t j = wuffs_base__u32__min(h->num_digits, precision + 1);
1794 for (; i < j; i++) {
1795 *ptr++ = (uint8_t)('0' | h->digits[i]);
1796 }
1797 for (; i <= precision; i++) {
1798 *ptr++ = '0';
1799 }
1800 }
1801
1802 // Exponent: "e±" and then 2 or 3 digits.
1803 *ptr++ = 'e';
1804 *ptr++ = negative_exp ? '-' : '+';
1805 if (exp < 10) {
1806 *ptr++ = '0';
1807 *ptr++ = (uint8_t)('0' | exp);
1808 } else if (exp < 100) {
1809 *ptr++ = (uint8_t)('0' | (exp / 10));
1810 *ptr++ = (uint8_t)('0' | (exp % 10));
1811 } else {
1812 int32_t e = exp / 100;
1813 exp -= e * 100;
1814 *ptr++ = (uint8_t)('0' | e);
1815 *ptr++ = (uint8_t)('0' | (exp / 10));
1816 *ptr++ = (uint8_t)('0' | (exp % 10));
1817 }
1818
1819 return n;
1820}
1821
1822WUFFS_BASE__MAYBE_STATIC size_t //
1823wuffs_base__render_number_f64(wuffs_base__slice_u8 dst,
1824 double x,
1825 uint32_t precision,
1826 uint32_t options) {
1827 // Decompose x (64 bits) into negativity (1 bit), base-2 exponent (11 bits
1828 // with a -1023 bias) and mantissa (52 bits).
Nigel Tao4d449dc2020-07-12 11:00:47 +10001829 uint64_t bits = wuffs_base__ieee_754_bit_representation__from_f64_to_u64(x);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001830 bool neg = (bits >> 63) != 0;
1831 int32_t exp2 = ((int32_t)(bits >> 52)) & 0x7FF;
1832 uint64_t man = bits & 0x000FFFFFFFFFFFFFul;
1833
1834 // Apply the exponent bias and set the implicit top bit of the mantissa,
1835 // unless x is subnormal. Also take care of Inf and NaN.
1836 if (exp2 == 0x7FF) {
1837 if (man != 0) {
1838 return wuffs_base__private_implementation__render_nan(dst);
1839 }
1840 return wuffs_base__private_implementation__render_inf(dst, neg, options);
1841 } else if (exp2 == 0) {
1842 exp2 = -1022;
1843 } else {
1844 exp2 -= 1023;
1845 man |= 0x0010000000000000ul;
1846 }
1847
1848 // Ensure that precision isn't too large.
1849 if (precision > 4095) {
1850 precision = 4095;
1851 }
1852
1853 // Convert from the (neg, exp2, man) tuple to an HPD.
1854 wuffs_base__private_implementation__high_prec_dec h;
1855 wuffs_base__private_implementation__high_prec_dec__assign(&h, man, neg);
1856 if (h.num_digits > 0) {
1857 wuffs_base__private_implementation__high_prec_dec__lshift(
1858 &h, exp2 - 52); // 52 mantissa bits.
1859 }
1860
1861 // Handle the "%e" and "%f" formats.
1862 switch (options & (WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT |
1863 WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT)) {
1864 case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT: // The "%"f" format.
1865 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1866 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1867 &h, exp2, man);
1868 int32_t p = ((int32_t)(h.num_digits)) - h.decimal_point;
1869 precision = ((uint32_t)(wuffs_base__i32__max(0, p)));
1870 } else {
1871 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1872 &h, ((int32_t)precision) + h.decimal_point);
1873 }
1874 return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1875 dst, &h, precision, options);
1876
1877 case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT: // The "%e" format.
1878 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1879 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1880 &h, exp2, man);
1881 precision = (h.num_digits > 0) ? (h.num_digits - 1) : 0;
1882 } else {
1883 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1884 &h, ((int32_t)precision) + 1);
1885 }
1886 return wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1887 dst, &h, precision, options);
1888 }
1889
1890 // We have the "%g" format and so precision means the number of significant
1891 // digits, not the number of digits after the decimal separator. Perform
1892 // rounding and determine whether to use "%e" or "%f".
1893 int32_t e_threshold = 0;
1894 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1895 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1896 &h, exp2, man);
1897 precision = h.num_digits;
1898 e_threshold = 6;
1899 } else {
1900 if (precision == 0) {
1901 precision = 1;
1902 }
1903 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1904 &h, ((int32_t)precision));
1905 e_threshold = ((int32_t)precision);
1906 int32_t nd = ((int32_t)(h.num_digits));
1907 if ((e_threshold > nd) && (nd >= h.decimal_point)) {
1908 e_threshold = nd;
1909 }
1910 }
1911
1912 // Use the "%e" format if the exponent is large.
1913 int32_t e = h.decimal_point - 1;
1914 if ((e < -4) || (e_threshold <= e)) {
1915 uint32_t p = wuffs_base__u32__min(precision, h.num_digits);
1916 return wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1917 dst, &h, (p > 0) ? (p - 1) : 0, options);
1918 }
1919
1920 // Use the "%f" format otherwise.
1921 int32_t p = ((int32_t)precision);
1922 if (p > h.decimal_point) {
1923 p = ((int32_t)(h.num_digits));
1924 }
1925 precision = ((uint32_t)(wuffs_base__i32__max(0, p - h.decimal_point)));
1926 return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1927 dst, &h, precision, options);
1928}