blob: df8ff26d2379161691adda7c60ae2f47e4cde7d5 [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/
978//
979// Preconditions:
980// - man is non-zero.
Nigel Tao09872962020-09-15 22:22:51 +1000981// - exp10 is in the range [-307 ..= 288], the same range of the
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000982// wuffs_base__private_implementation__powers_of_10 array.
Nigel Tao8b45db02020-09-15 21:50:32 +1000983//
984// The exp10 range (and the fact that man is in the range [1 ..= UINT64_MAX],
985// approximately [1 ..= 1.85e+19]) means that (man * (10 ** exp10)) is in the
986// range [1e-307 ..= 1.85e+307]. This is entirely within the range of normal
987// (neither subnormal nor non-finite) f64 values: DBL_MIN and DBL_MAX are
988// approximately 2.23e–308 and 1.80e+308.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000989static int64_t //
Nigel Taoc4fa8e22020-07-18 17:35:13 +1000990wuffs_base__private_implementation__parse_number_f64_eisel_lemire(
991 uint64_t man,
992 int32_t exp10) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000993 // Look up the (possibly truncated) base-2 representation of (10 ** exp10).
994 // The look-up table was constructed so that it is already normalized: the
995 // table entry's mantissa's MSB (most significant bit) is on.
Nigel Taoafe7f272020-09-23 15:52:13 +1000996 const uint64_t* po10 =
997 &wuffs_base__private_implementation__powers_of_10[exp10 + 307][0];
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000998
999 // Normalize the man argument. The (man != 0) precondition means that a
1000 // non-zero bit exists.
1001 uint32_t clz = wuffs_base__count_leading_zeroes_u64(man);
1002 man <<= clz;
1003
1004 // Calculate the return value's base-2 exponent. We might tweak it by ±1
Nigel Taob6d85522020-09-23 15:21:47 +10001005 // later, but its initial value comes from a linear scaling of exp10,
1006 // converting from power-of-10 to power-of-2, and adjusting by clz.
1007 //
1008 // The magic constants are:
1009 // - 1087 = 1023 + 64. The 1023 is the f64 exponent bias. The 64 is because
1010 // the look-up table uses 64-bit mantissas.
1011 // - 217706 is such that the ratio 217706 / 65536 ≈ 3.321930 is close enough
1012 // (over the practical range of exp10) to log(10) / log(2) ≈ 3.321928.
1013 // - 65536 = 1<<16 is arbitrary but a power of 2, so division is a shift.
1014 //
1015 // Equality of the linearly-scaled value and the actual power-of-2, over the
1016 // range of exp10 arguments that this function accepts, is confirmed by
1017 // script/print-mpb-powers-of-10.go
1018 uint64_t ret_exp2 =
1019 ((uint64_t)(((217706 * exp10) >> 16) + 1087)) - ((uint64_t)clz);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001020
1021 // Multiply the two mantissas. Normalization means that both mantissas are at
1022 // least (1<<63), so the 128-bit product must be at least (1<<126). The high
Nigel Tao74d4af62020-07-10 11:27:17 +10001023 // 64 bits of the product, x_hi, must therefore be at least (1<<62).
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001024 //
Nigel Tao74d4af62020-07-10 11:27:17 +10001025 // As a consequence, x_hi has either 0 or 1 leading zeroes. Shifting x_hi
1026 // right by either 9 or 10 bits (depending on x_hi's MSB) will therefore
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001027 // leave the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on.
Nigel Taoafe7f272020-09-23 15:52:13 +10001028 wuffs_base__multiply_u64__output x = wuffs_base__multiply_u64(man, po10[1]);
Nigel Tao74d4af62020-07-10 11:27:17 +10001029 uint64_t x_hi = x.hi;
1030 uint64_t x_lo = x.lo;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001031
1032 // Before we shift right by at least 9 bits, recall that the look-up table
1033 // entry was possibly truncated. We have so far only calculated a lower bound
1034 // for the product (man * e), where e is (10 ** exp10). The upper bound would
1035 // add a further (man * 1) to the 128-bit product, which overflows the lower
Nigel Tao74d4af62020-07-10 11:27:17 +10001036 // 64-bit limb if ((x_lo + man) < man).
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001037 //
Nigel Tao74d4af62020-07-10 11:27:17 +10001038 // If overflow occurs, that adds 1 to x_hi. Since we're about to shift right
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001039 // by at least 9 bits, that carried 1 can be ignored unless the higher 64-bit
1040 // limb's low 9 bits are all on.
Nigel Tao74d4af62020-07-10 11:27:17 +10001041 if (((x_hi & 0x1FF) == 0x1FF) && ((x_lo + man) < man)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001042 // Refine our calculation of (man * e). Before, our approximation of e used
1043 // a "low resolution" 64-bit mantissa. Now use a "high resolution" 128-bit
1044 // mantissa. We've already calculated x = (man * bits_0_to_63_incl_of_e).
1045 // Now calculate y = (man * bits_64_to_127_incl_of_e).
Nigel Taoafe7f272020-09-23 15:52:13 +10001046 wuffs_base__multiply_u64__output y = wuffs_base__multiply_u64(man, po10[0]);
Nigel Tao74d4af62020-07-10 11:27:17 +10001047 uint64_t y_hi = y.hi;
1048 uint64_t y_lo = y.lo;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001049
1050 // Merge the 128-bit x and 128-bit y, which overlap by 64 bits, to
1051 // calculate the 192-bit product of the 64-bit man by the 128-bit e.
1052 // As we exit this if-block, we only care about the high 128 bits
1053 // (merged_hi and merged_lo) of that 192-bit product.
Nigel Tao74d4af62020-07-10 11:27:17 +10001054 uint64_t merged_hi = x_hi;
1055 uint64_t merged_lo = x_lo + y_hi;
1056 if (merged_lo < x_lo) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001057 merged_hi++; // Carry the overflow bit.
1058 }
1059
1060 // The "high resolution" approximation of e is still a lower bound. Once
1061 // again, see if the upper bound is large enough to produce a different
1062 // result. This time, if it does, give up instead of reaching for an even
1063 // more precise approximation to e.
1064 //
1065 // This three-part check is similar to the two-part check that guarded the
1066 // if block that we're now in, but it has an extra term for the middle 64
1067 // bits (checking that adding 1 to merged_lo would overflow).
1068 if (((merged_hi & 0x1FF) == 0x1FF) && ((merged_lo + 1) == 0) &&
Nigel Tao74d4af62020-07-10 11:27:17 +10001069 (y_lo + man < man)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001070 return -1;
1071 }
1072
1073 // Replace the 128-bit x with merged.
Nigel Tao74d4af62020-07-10 11:27:17 +10001074 x_hi = merged_hi;
1075 x_lo = merged_lo;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001076 }
1077
Nigel Tao74d4af62020-07-10 11:27:17 +10001078 // As mentioned above, shifting x_hi right by either 9 or 10 bits will leave
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001079 // the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on. If the
1080 // MSB (before shifting) was on, adjust ret_exp2 for the larger shift.
1081 //
1082 // Having bit 53 on (and higher bits off) means that ret_mantissa is a 54-bit
1083 // number.
Nigel Tao74d4af62020-07-10 11:27:17 +10001084 uint64_t msb = x_hi >> 63;
1085 uint64_t ret_mantissa = x_hi >> (msb + 9);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001086 ret_exp2 -= 1 ^ msb;
1087
1088 // IEEE 754 rounds to-nearest with ties rounded to-even. Rounding to-even can
1089 // be tricky. If we're half-way between two exactly representable numbers
1090 // (x's low 73 bits are zero and the next 2 bits that matter are "01"), give
1091 // up instead of trying to pick the winner.
1092 //
1093 // Technically, we could tighten the condition by changing "73" to "73 or 74,
1094 // depending on msb", but a flat "73" is simpler.
Nigel Tao74d4af62020-07-10 11:27:17 +10001095 if ((x_lo == 0) && ((x_hi & 0x1FF) == 0) && ((ret_mantissa & 3) == 1)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001096 return -1;
1097 }
1098
1099 // If we're not halfway then it's rounding to-nearest. Starting with a 54-bit
1100 // number, carry the lowest bit (bit 0) up if it's on. Regardless of whether
1101 // it was on or off, shifting right by one then produces a 53-bit number. If
1102 // carrying up overflowed, shift again.
1103 ret_mantissa += ret_mantissa & 1;
1104 ret_mantissa >>= 1;
Nigel Tao8b45db02020-09-15 21:50:32 +10001105 // This if block is equivalent to (but benchmarks slightly faster than) the
1106 // following branchless form:
1107 // uint64_t overflow_adjustment = ret_mantissa >> 53;
1108 // ret_mantissa >>= overflow_adjustment;
1109 // ret_exp2 += overflow_adjustment;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001110 if ((ret_mantissa >> 53) > 0) {
1111 ret_mantissa >>= 1;
1112 ret_exp2++;
1113 }
1114
1115 // Starting with a 53-bit number, IEEE 754 double-precision normal numbers
1116 // have an implicit mantissa bit. Mask that away and keep the low 52 bits.
1117 ret_mantissa &= 0x000FFFFFFFFFFFFF;
1118
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001119 // Pack the bits and return.
1120 return ((int64_t)(ret_mantissa | (ret_exp2 << 52)));
1121}
1122
1123// --------
1124
1125static wuffs_base__result_f64 //
Nigel Taoe0c5de92020-07-11 11:48:17 +10001126wuffs_base__private_implementation__parse_number_f64_special(
1127 wuffs_base__slice_u8 s,
Nigel Tao4d61a052020-07-11 12:34:40 +10001128 uint32_t options) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001129 do {
Nigel Tao4d61a052020-07-11 12:34:40 +10001130 if (options & WUFFS_BASE__PARSE_NUMBER_FXX__REJECT_INF_AND_NAN) {
1131 goto fail;
1132 }
1133
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001134 uint8_t* p = s.ptr;
1135 uint8_t* q = s.ptr + s.len;
1136
1137 for (; (p < q) && (*p == '_'); p++) {
1138 }
1139 if (p >= q) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001140 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001141 }
1142
1143 // Parse sign.
1144 bool negative = false;
1145 do {
1146 if (*p == '+') {
1147 p++;
1148 } else if (*p == '-') {
1149 negative = true;
1150 p++;
1151 } else {
1152 break;
1153 }
1154 for (; (p < q) && (*p == '_'); p++) {
1155 }
1156 } while (0);
1157 if (p >= q) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001158 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001159 }
1160
1161 bool nan = false;
1162 switch (p[0]) {
1163 case 'I':
1164 case 'i':
1165 if (((q - p) < 3) || //
1166 ((p[1] != 'N') && (p[1] != 'n')) || //
1167 ((p[2] != 'F') && (p[2] != 'f'))) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001168 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001169 }
1170 p += 3;
1171
1172 if ((p >= q) || (*p == '_')) {
1173 break;
1174 } else if (((q - p) < 5) || //
1175 ((p[0] != 'I') && (p[0] != 'i')) || //
1176 ((p[1] != 'N') && (p[1] != 'n')) || //
1177 ((p[2] != 'I') && (p[2] != 'i')) || //
1178 ((p[3] != 'T') && (p[3] != 't')) || //
1179 ((p[4] != 'Y') && (p[4] != 'y'))) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001180 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001181 }
1182 p += 5;
1183
1184 if ((p >= q) || (*p == '_')) {
1185 break;
1186 }
Nigel Tao4d61a052020-07-11 12:34:40 +10001187 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001188
1189 case 'N':
1190 case 'n':
1191 if (((q - p) < 3) || //
1192 ((p[1] != 'A') && (p[1] != 'a')) || //
1193 ((p[2] != 'N') && (p[2] != 'n'))) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001194 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001195 }
1196 p += 3;
1197
1198 if ((p >= q) || (*p == '_')) {
1199 nan = true;
1200 break;
1201 }
Nigel Tao4d61a052020-07-11 12:34:40 +10001202 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001203
1204 default:
Nigel Tao4d61a052020-07-11 12:34:40 +10001205 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001206 }
1207
1208 // Finish.
1209 for (; (p < q) && (*p == '_'); p++) {
1210 }
1211 if (p != q) {
Nigel Tao4d61a052020-07-11 12:34:40 +10001212 goto fail;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001213 }
1214 wuffs_base__result_f64 ret;
1215 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001216 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001217 (nan ? 0x7FFFFFFFFFFFFFFF : 0x7FF0000000000000) |
1218 (negative ? 0x8000000000000000 : 0));
1219 return ret;
1220 } while (0);
1221
Nigel Tao4d61a052020-07-11 12:34:40 +10001222fail:
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001223 do {
1224 wuffs_base__result_f64 ret;
Nigel Tao4d61a052020-07-11 12:34:40 +10001225 ret.status.repr = wuffs_base__error__bad_argument;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001226 ret.value = 0;
1227 return ret;
1228 } while (0);
1229}
1230
1231WUFFS_BASE__MAYBE_STATIC wuffs_base__result_f64 //
Nigel Taoe0c5de92020-07-11 11:48:17 +10001232wuffs_base__private_implementation__high_prec_dec__to_f64(
Nigel Tao4d61a052020-07-11 12:34:40 +10001233 wuffs_base__private_implementation__high_prec_dec* h,
1234 uint32_t options) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001235 do {
1236 // powers converts decimal powers of 10 to binary powers of 2. For example,
1237 // (10000 >> 13) is 1. It stops before the elements exceed 60, also known
1238 // as WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL.
1239 static const uint32_t num_powers = 19;
1240 static const uint8_t powers[19] = {
1241 0, 3, 6, 9, 13, 16, 19, 23, 26, 29, //
1242 33, 36, 39, 43, 46, 49, 53, 56, 59, //
1243 };
1244
1245 // Handle zero and obvious extremes. The largest and smallest positive
1246 // finite f64 values are approximately 1.8e+308 and 4.9e-324.
1247 if ((h->num_digits == 0) || (h->decimal_point < -326)) {
1248 goto zero;
1249 } else if (h->decimal_point > 310) {
1250 goto infinity;
1251 }
1252
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001253 // Try the fast Eisel-Lemire algorithm again. Calculating the (man, exp10)
1254 // pair from the high_prec_dec h is more correct but slower than the
1255 // approach taken in wuffs_base__parse_number_f64. The latter is optimized
1256 // for the common cases (e.g. assuming no underscores or a leading '+'
1257 // sign) rather than the full set of cases allowed by the Wuffs API.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001258 if (h->num_digits <= 19) {
1259 uint64_t man = 0;
1260 uint32_t i;
1261 for (i = 0; i < h->num_digits; i++) {
1262 man = (10 * man) + h->digits[i];
1263 }
1264 int32_t exp10 = h->decimal_point - ((int32_t)(h->num_digits));
Nigel Tao8b45db02020-09-15 21:50:32 +10001265 if ((man != 0) && (-307 <= exp10) && (exp10 <= 288)) {
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001266 int64_t r =
1267 wuffs_base__private_implementation__parse_number_f64_eisel_lemire(
1268 man, exp10);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001269 if (r >= 0) {
1270 wuffs_base__result_f64 ret;
1271 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001272 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001273 ((uint64_t)r) | (((uint64_t)(h->negative)) << 63));
1274 return ret;
1275 }
1276 }
1277 }
1278
1279 // Scale by powers of 2 until we're in the range [½ .. 1], which gives us
1280 // our exponent (in base-2). First we shift right, possibly a little too
1281 // far, ending with a value certainly below 1 and possibly below ½...
1282 const int32_t f64_bias = -1023;
1283 int32_t exp2 = 0;
1284 while (h->decimal_point > 0) {
1285 uint32_t n = (uint32_t)(+h->decimal_point);
1286 uint32_t shift =
1287 (n < num_powers)
1288 ? powers[n]
1289 : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1290
1291 wuffs_base__private_implementation__high_prec_dec__small_rshift(h, shift);
1292 if (h->decimal_point <
1293 -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
1294 goto zero;
1295 }
1296 exp2 += (int32_t)shift;
1297 }
1298 // ...then we shift left, putting us in [½ .. 1].
1299 while (h->decimal_point <= 0) {
1300 uint32_t shift;
1301 if (h->decimal_point == 0) {
1302 if (h->digits[0] >= 5) {
1303 break;
1304 }
Nigel Tao57d47c62020-09-08 16:43:31 +10001305 shift = (h->digits[0] < 2) ? 2 : 1;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001306 } else {
1307 uint32_t n = (uint32_t)(-h->decimal_point);
1308 shift = (n < num_powers)
1309 ? powers[n]
1310 : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1311 }
1312
1313 wuffs_base__private_implementation__high_prec_dec__small_lshift(h, shift);
1314 if (h->decimal_point >
1315 +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
1316 goto infinity;
1317 }
1318 exp2 -= (int32_t)shift;
1319 }
1320
1321 // We're in the range [½ .. 1] but f64 uses [1 .. 2].
1322 exp2--;
1323
1324 // The minimum normal exponent is (f64_bias + 1).
1325 while ((f64_bias + 1) > exp2) {
1326 uint32_t n = (uint32_t)((f64_bias + 1) - exp2);
1327 if (n > WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) {
1328 n = WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1329 }
1330 wuffs_base__private_implementation__high_prec_dec__small_rshift(h, n);
1331 exp2 += (int32_t)n;
1332 }
1333
1334 // Check for overflow.
1335 if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1.
1336 goto infinity;
1337 }
1338
1339 // Extract 53 bits for the mantissa (in base-2).
1340 wuffs_base__private_implementation__high_prec_dec__small_lshift(h, 53);
1341 uint64_t man2 =
1342 wuffs_base__private_implementation__high_prec_dec__rounded_integer(h);
1343
1344 // Rounding might have added one bit. If so, shift and re-check overflow.
1345 if ((man2 >> 53) != 0) {
1346 man2 >>= 1;
1347 exp2++;
1348 if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1.
1349 goto infinity;
1350 }
1351 }
1352
1353 // Handle subnormal numbers.
1354 if ((man2 >> 52) == 0) {
1355 exp2 = f64_bias;
1356 }
1357
1358 // Pack the bits and return.
1359 uint64_t exp2_bits =
1360 (uint64_t)((exp2 - f64_bias) & 0x07FF); // (1 << 11) - 1.
1361 uint64_t bits = (man2 & 0x000FFFFFFFFFFFFF) | // (1 << 52) - 1.
1362 (exp2_bits << 52) | //
1363 (h->negative ? 0x8000000000000000 : 0); // (1 << 63).
1364
1365 wuffs_base__result_f64 ret;
1366 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001367 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001368 return ret;
1369 } while (0);
1370
1371zero:
1372 do {
1373 uint64_t bits = h->negative ? 0x8000000000000000 : 0;
1374
1375 wuffs_base__result_f64 ret;
1376 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001377 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001378 return ret;
1379 } while (0);
1380
1381infinity:
1382 do {
Nigel Tao4d61a052020-07-11 12:34:40 +10001383 if (options & WUFFS_BASE__PARSE_NUMBER_FXX__REJECT_INF_AND_NAN) {
1384 wuffs_base__result_f64 ret;
1385 ret.status.repr = wuffs_base__error__bad_argument;
1386 ret.value = 0;
1387 return ret;
1388 }
1389
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001390 uint64_t bits = h->negative ? 0xFFF0000000000000 : 0x7FF0000000000000;
1391
1392 wuffs_base__result_f64 ret;
1393 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001394 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001395 return ret;
1396 } while (0);
1397}
1398
1399static inline bool //
1400wuffs_base__private_implementation__is_decimal_digit(uint8_t c) {
1401 return ('0' <= c) && (c <= '9');
1402}
1403
1404WUFFS_BASE__MAYBE_STATIC wuffs_base__result_f64 //
1405wuffs_base__parse_number_f64(wuffs_base__slice_u8 s, uint32_t options) {
1406 // In practice, almost all "dd.ddddE±xxx" numbers can be represented
1407 // losslessly by a uint64_t mantissa "dddddd" and an int32_t base-10
1408 // exponent, adjusting "xxx" for the position (if present) of the decimal
1409 // separator '.' or ','.
1410 //
1411 // This (u64 man, i32 exp10) data structure is superficially similar to the
1412 // "Do It Yourself Floating Point" type from Loitsch (†), but the exponent
1413 // here is base-10, not base-2.
1414 //
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001415 // If s's number fits in a (man, exp10), parse that pair with the
1416 // Eisel-Lemire algorithm. If not, or if Eisel-Lemire fails, parsing s with
1417 // the fallback algorithm is slower but comprehensive.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001418 //
1419 // † "Printing Floating-Point Numbers Quickly and Accurately with Integers"
1420 // (https://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf).
1421 // Florian Loitsch is also the primary contributor to
1422 // https://github.com/google/double-conversion
1423 do {
1424 // Calculating that (man, exp10) pair needs to stay within s's bounds.
1425 // Provided that s isn't extremely long, work on a NUL-terminated copy of
1426 // s's contents. The NUL byte isn't a valid part of "±dd.ddddE±xxx".
1427 //
1428 // As the pointer p walks the contents, it's faster to repeatedly check "is
1429 // *p a valid digit" than "is p within bounds and *p a valid digit".
1430 if (s.len >= 256) {
1431 goto fallback;
1432 }
1433 uint8_t z[256];
1434 memcpy(&z[0], s.ptr, s.len);
1435 z[s.len] = 0;
1436 const uint8_t* p = &z[0];
1437
1438 // Look for a leading minus sign. Technically, we could also look for an
1439 // optional plus sign, but the "script/process-json-numbers.c with -p"
1440 // benchmark is noticably slower if we do. It's optional and, in practice,
1441 // usually absent. Let the fallback catch it.
1442 bool negative = (*p == '-');
1443 if (negative) {
1444 p++;
1445 }
1446
1447 // After walking "dd.dddd", comparing p later with p now will produce the
1448 // number of "d"s and "."s.
1449 const uint8_t* const start_of_digits_ptr = p;
1450
1451 // Walk the "d"s before a '.', 'E', NUL byte, etc. If it starts with '0',
1452 // it must be a single '0'. If it starts with a non-zero decimal digit, it
1453 // can be a sequence of decimal digits.
1454 //
1455 // Update the man variable during the walk. It's OK if man overflows now.
1456 // We'll detect that later.
1457 uint64_t man;
1458 if (*p == '0') {
1459 man = 0;
1460 p++;
1461 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1462 goto fallback;
1463 }
1464 } else if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1465 man = ((uint8_t)(*p - '0'));
1466 p++;
1467 for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) {
1468 man = (10 * man) + ((uint8_t)(*p - '0'));
1469 }
1470 } else {
1471 goto fallback;
1472 }
1473
1474 // Walk the "d"s after the optional decimal separator ('.' or ','),
1475 // updating the man and exp10 variables.
1476 int32_t exp10 = 0;
Nigel Taoe0c5de92020-07-11 11:48:17 +10001477 if (*p ==
1478 ((options & WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1479 ? ','
1480 : '.')) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001481 p++;
1482 const uint8_t* first_after_separator_ptr = p;
1483 if (!wuffs_base__private_implementation__is_decimal_digit(*p)) {
1484 goto fallback;
1485 }
1486 man = (10 * man) + ((uint8_t)(*p - '0'));
1487 p++;
1488 for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) {
1489 man = (10 * man) + ((uint8_t)(*p - '0'));
1490 }
1491 exp10 = ((int32_t)(first_after_separator_ptr - p));
1492 }
1493
1494 // Count the number of digits:
1495 // - for an input of "314159", digit_count is 6.
1496 // - for an input of "3.14159", digit_count is 7.
1497 //
1498 // This is off-by-one if there is a decimal separator. That's OK for now.
1499 // We'll correct for that later. The "script/process-json-numbers.c with
1500 // -p" benchmark is noticably slower if we try to correct for that now.
1501 uint32_t digit_count = (uint32_t)(p - start_of_digits_ptr);
1502
1503 // Update exp10 for the optional exponent, starting with 'E' or 'e'.
1504 if ((*p | 0x20) == 'e') {
1505 p++;
1506 int32_t exp_sign = +1;
1507 if (*p == '-') {
1508 p++;
1509 exp_sign = -1;
1510 } else if (*p == '+') {
1511 p++;
1512 }
1513 if (!wuffs_base__private_implementation__is_decimal_digit(*p)) {
1514 goto fallback;
1515 }
1516 int32_t exp_num = ((uint8_t)(*p - '0'));
1517 p++;
1518 // The rest of the exp_num walking has a peculiar control flow but, once
1519 // again, the "script/process-json-numbers.c with -p" benchmark is
1520 // sensitive to alternative formulations.
1521 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1522 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1523 p++;
1524 }
1525 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1526 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1527 p++;
1528 }
1529 while (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1530 if (exp_num > 0x1000000) {
1531 goto fallback;
1532 }
1533 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1534 p++;
1535 }
1536 exp10 += exp_sign * exp_num;
1537 }
1538
1539 // The Wuffs API is that the original slice has no trailing data. It also
1540 // allows underscores, which we don't catch here but the fallback should.
1541 if (p != &z[s.len]) {
1542 goto fallback;
1543 }
1544
1545 // Check that the uint64_t typed man variable has not overflowed, based on
1546 // digit_count.
1547 //
1548 // For reference:
1549 // - (1 << 63) is 9223372036854775808, which has 19 decimal digits.
1550 // - (1 << 64) is 18446744073709551616, which has 20 decimal digits.
1551 // - 19 nines, 9999999999999999999, is 0x8AC7230489E7FFFF, which has 64
1552 // bits and 16 hexadecimal digits.
1553 // - 20 nines, 99999999999999999999, is 0x56BC75E2D630FFFFF, which has 67
1554 // bits and 17 hexadecimal digits.
1555 if (digit_count > 19) {
1556 // Even if we have more than 19 pseudo-digits, it's not yet definitely an
1557 // overflow. Recall that digit_count might be off-by-one (too large) if
1558 // there's a decimal separator. It will also over-report the number of
1559 // meaningful digits if the input looks something like "0.000dddExxx".
1560 //
1561 // We adjust by the number of leading '0's and '.'s and re-compare to 19.
1562 // Once again, technically, we could skip ','s too, but that perturbs the
1563 // "script/process-json-numbers.c with -p" benchmark.
1564 const uint8_t* q = start_of_digits_ptr;
1565 for (; (*q == '0') || (*q == '.'); q++) {
1566 }
1567 digit_count -= (uint32_t)(q - start_of_digits_ptr);
1568 if (digit_count > 19) {
1569 goto fallback;
1570 }
1571 }
1572
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001573 // The wuffs_base__private_implementation__parse_number_f64_eisel_lemire
Nigel Tao8b45db02020-09-15 21:50:32 +10001574 // preconditions include that exp10 is in the range [-307 ..= 288].
1575 if ((exp10 < -307) || (288 < exp10)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001576 goto fallback;
1577 }
1578
Nigel Tao9f22b5e2020-09-11 09:10:08 +10001579 // If both man and (10 ** exp10) are exactly representable by a double, we
1580 // don't need to run the Eisel-Lemire algorithm.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001581 if ((-22 <= exp10) && (exp10 <= 22) && ((man >> 53) == 0)) {
1582 double d = (double)man;
1583 if (exp10 >= 0) {
1584 d *= wuffs_base__private_implementation__f64_powers_of_10[+exp10];
1585 } else {
1586 d /= wuffs_base__private_implementation__f64_powers_of_10[-exp10];
1587 }
1588 wuffs_base__result_f64 ret;
1589 ret.status.repr = NULL;
1590 ret.value = negative ? -d : +d;
1591 return ret;
1592 }
1593
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001594 // The wuffs_base__private_implementation__parse_number_f64_eisel_lemire
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001595 // preconditions include that man is non-zero. Parsing "0" should be caught
Nigel Tao9f22b5e2020-09-11 09:10:08 +10001596 // by the "If both man and (10 ** exp10)" above, but "0e99" might not.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001597 if (man == 0) {
1598 goto fallback;
1599 }
1600
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001601 // Our man and exp10 are in range. Run the Eisel-Lemire algorithm.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001602 int64_t r =
Nigel Taoc4fa8e22020-07-18 17:35:13 +10001603 wuffs_base__private_implementation__parse_number_f64_eisel_lemire(
1604 man, exp10);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001605 if (r < 0) {
1606 goto fallback;
1607 }
1608 wuffs_base__result_f64 ret;
1609 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001610 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001611 ((uint64_t)r) | (((uint64_t)negative) << 63));
1612 return ret;
1613 } while (0);
1614
1615fallback:
1616 do {
1617 wuffs_base__private_implementation__high_prec_dec h;
1618 wuffs_base__status status =
Nigel Taoe0c5de92020-07-11 11:48:17 +10001619 wuffs_base__private_implementation__high_prec_dec__parse(&h, s,
1620 options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001621 if (status.repr) {
Nigel Taoe0c5de92020-07-11 11:48:17 +10001622 return wuffs_base__private_implementation__parse_number_f64_special(
Nigel Tao4d61a052020-07-11 12:34:40 +10001623 s, options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001624 }
Nigel Tao4d61a052020-07-11 12:34:40 +10001625 return wuffs_base__private_implementation__high_prec_dec__to_f64(&h,
1626 options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001627 } while (0);
1628}
1629
1630// --------
1631
1632static inline size_t //
1633wuffs_base__private_implementation__render_inf(wuffs_base__slice_u8 dst,
1634 bool neg,
1635 uint32_t options) {
1636 if (neg) {
1637 if (dst.len < 4) {
1638 return 0;
1639 }
1640 wuffs_base__store_u32le__no_bounds_check(dst.ptr, 0x666E492D); // '-Inf'le.
1641 return 4;
1642 }
1643
1644 if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1645 if (dst.len < 4) {
1646 return 0;
1647 }
1648 wuffs_base__store_u32le__no_bounds_check(dst.ptr, 0x666E492B); // '+Inf'le.
1649 return 4;
1650 }
1651
1652 if (dst.len < 3) {
1653 return 0;
1654 }
1655 wuffs_base__store_u24le__no_bounds_check(dst.ptr, 0x666E49); // 'Inf'le.
1656 return 3;
1657}
1658
1659static inline size_t //
1660wuffs_base__private_implementation__render_nan(wuffs_base__slice_u8 dst) {
1661 if (dst.len < 3) {
1662 return 0;
1663 }
1664 wuffs_base__store_u24le__no_bounds_check(dst.ptr, 0x4E614E); // 'NaN'le.
1665 return 3;
1666}
1667
1668static size_t //
1669wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1670 wuffs_base__slice_u8 dst,
1671 wuffs_base__private_implementation__high_prec_dec* h,
1672 uint32_t precision,
1673 uint32_t options) {
1674 size_t n = (h->negative ||
1675 (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN))
1676 ? 1
1677 : 0;
1678 if (h->decimal_point <= 0) {
1679 n += 1;
1680 } else {
1681 n += (size_t)(h->decimal_point);
1682 }
1683 if (precision > 0) {
1684 n += precision + 1; // +1 for the '.'.
1685 }
1686
1687 // Don't modify dst if the formatted number won't fit.
1688 if (n > dst.len) {
1689 return 0;
1690 }
1691
1692 // Align-left or align-right.
1693 uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT)
1694 ? &dst.ptr[dst.len - n]
1695 : &dst.ptr[0];
1696
1697 // Leading "±".
1698 if (h->negative) {
1699 *ptr++ = '-';
1700 } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1701 *ptr++ = '+';
1702 }
1703
1704 // Integral digits.
1705 if (h->decimal_point <= 0) {
1706 *ptr++ = '0';
1707 } else {
1708 uint32_t m =
1709 wuffs_base__u32__min(h->num_digits, (uint32_t)(h->decimal_point));
1710 uint32_t i = 0;
1711 for (; i < m; i++) {
1712 *ptr++ = (uint8_t)('0' | h->digits[i]);
1713 }
1714 for (; i < (uint32_t)(h->decimal_point); i++) {
1715 *ptr++ = '0';
1716 }
1717 }
1718
1719 // Separator and then fractional digits.
1720 if (precision > 0) {
1721 *ptr++ =
1722 (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1723 ? ','
1724 : '.';
1725 uint32_t i = 0;
1726 for (; i < precision; i++) {
1727 uint32_t j = ((uint32_t)(h->decimal_point)) + i;
1728 *ptr++ = (uint8_t)('0' | ((j < h->num_digits) ? h->digits[j] : 0));
1729 }
1730 }
1731
1732 return n;
1733}
1734
1735static size_t //
1736wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1737 wuffs_base__slice_u8 dst,
1738 wuffs_base__private_implementation__high_prec_dec* h,
1739 uint32_t precision,
1740 uint32_t options) {
1741 int32_t exp = 0;
1742 if (h->num_digits > 0) {
1743 exp = h->decimal_point - 1;
1744 }
1745 bool negative_exp = exp < 0;
1746 if (negative_exp) {
1747 exp = -exp;
1748 }
1749
1750 size_t n = (h->negative ||
1751 (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN))
1752 ? 4
1753 : 3; // Mininum 3 bytes: first digit and then "e±".
1754 if (precision > 0) {
1755 n += precision + 1; // +1 for the '.'.
1756 }
1757 n += (exp < 100) ? 2 : 3;
1758
1759 // Don't modify dst if the formatted number won't fit.
1760 if (n > dst.len) {
1761 return 0;
1762 }
1763
1764 // Align-left or align-right.
1765 uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT)
1766 ? &dst.ptr[dst.len - n]
1767 : &dst.ptr[0];
1768
1769 // Leading "±".
1770 if (h->negative) {
1771 *ptr++ = '-';
1772 } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1773 *ptr++ = '+';
1774 }
1775
1776 // Integral digit.
1777 if (h->num_digits > 0) {
1778 *ptr++ = (uint8_t)('0' | h->digits[0]);
1779 } else {
1780 *ptr++ = '0';
1781 }
1782
1783 // Separator and then fractional digits.
1784 if (precision > 0) {
1785 *ptr++ =
1786 (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1787 ? ','
1788 : '.';
1789 uint32_t i = 1;
1790 uint32_t j = wuffs_base__u32__min(h->num_digits, precision + 1);
1791 for (; i < j; i++) {
1792 *ptr++ = (uint8_t)('0' | h->digits[i]);
1793 }
1794 for (; i <= precision; i++) {
1795 *ptr++ = '0';
1796 }
1797 }
1798
1799 // Exponent: "e±" and then 2 or 3 digits.
1800 *ptr++ = 'e';
1801 *ptr++ = negative_exp ? '-' : '+';
1802 if (exp < 10) {
1803 *ptr++ = '0';
1804 *ptr++ = (uint8_t)('0' | exp);
1805 } else if (exp < 100) {
1806 *ptr++ = (uint8_t)('0' | (exp / 10));
1807 *ptr++ = (uint8_t)('0' | (exp % 10));
1808 } else {
1809 int32_t e = exp / 100;
1810 exp -= e * 100;
1811 *ptr++ = (uint8_t)('0' | e);
1812 *ptr++ = (uint8_t)('0' | (exp / 10));
1813 *ptr++ = (uint8_t)('0' | (exp % 10));
1814 }
1815
1816 return n;
1817}
1818
1819WUFFS_BASE__MAYBE_STATIC size_t //
1820wuffs_base__render_number_f64(wuffs_base__slice_u8 dst,
1821 double x,
1822 uint32_t precision,
1823 uint32_t options) {
1824 // Decompose x (64 bits) into negativity (1 bit), base-2 exponent (11 bits
1825 // with a -1023 bias) and mantissa (52 bits).
Nigel Tao4d449dc2020-07-12 11:00:47 +10001826 uint64_t bits = wuffs_base__ieee_754_bit_representation__from_f64_to_u64(x);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001827 bool neg = (bits >> 63) != 0;
1828 int32_t exp2 = ((int32_t)(bits >> 52)) & 0x7FF;
1829 uint64_t man = bits & 0x000FFFFFFFFFFFFFul;
1830
1831 // Apply the exponent bias and set the implicit top bit of the mantissa,
1832 // unless x is subnormal. Also take care of Inf and NaN.
1833 if (exp2 == 0x7FF) {
1834 if (man != 0) {
1835 return wuffs_base__private_implementation__render_nan(dst);
1836 }
1837 return wuffs_base__private_implementation__render_inf(dst, neg, options);
1838 } else if (exp2 == 0) {
1839 exp2 = -1022;
1840 } else {
1841 exp2 -= 1023;
1842 man |= 0x0010000000000000ul;
1843 }
1844
1845 // Ensure that precision isn't too large.
1846 if (precision > 4095) {
1847 precision = 4095;
1848 }
1849
1850 // Convert from the (neg, exp2, man) tuple to an HPD.
1851 wuffs_base__private_implementation__high_prec_dec h;
1852 wuffs_base__private_implementation__high_prec_dec__assign(&h, man, neg);
1853 if (h.num_digits > 0) {
1854 wuffs_base__private_implementation__high_prec_dec__lshift(
1855 &h, exp2 - 52); // 52 mantissa bits.
1856 }
1857
1858 // Handle the "%e" and "%f" formats.
1859 switch (options & (WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT |
1860 WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT)) {
1861 case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT: // The "%"f" format.
1862 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1863 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1864 &h, exp2, man);
1865 int32_t p = ((int32_t)(h.num_digits)) - h.decimal_point;
1866 precision = ((uint32_t)(wuffs_base__i32__max(0, p)));
1867 } else {
1868 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1869 &h, ((int32_t)precision) + h.decimal_point);
1870 }
1871 return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1872 dst, &h, precision, options);
1873
1874 case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT: // The "%e" format.
1875 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1876 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1877 &h, exp2, man);
1878 precision = (h.num_digits > 0) ? (h.num_digits - 1) : 0;
1879 } else {
1880 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1881 &h, ((int32_t)precision) + 1);
1882 }
1883 return wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1884 dst, &h, precision, options);
1885 }
1886
1887 // We have the "%g" format and so precision means the number of significant
1888 // digits, not the number of digits after the decimal separator. Perform
1889 // rounding and determine whether to use "%e" or "%f".
1890 int32_t e_threshold = 0;
1891 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1892 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1893 &h, exp2, man);
1894 precision = h.num_digits;
1895 e_threshold = 6;
1896 } else {
1897 if (precision == 0) {
1898 precision = 1;
1899 }
1900 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1901 &h, ((int32_t)precision));
1902 e_threshold = ((int32_t)precision);
1903 int32_t nd = ((int32_t)(h.num_digits));
1904 if ((e_threshold > nd) && (nd >= h.decimal_point)) {
1905 e_threshold = nd;
1906 }
1907 }
1908
1909 // Use the "%e" format if the exponent is large.
1910 int32_t e = h.decimal_point - 1;
1911 if ((e < -4) || (e_threshold <= e)) {
1912 uint32_t p = wuffs_base__u32__min(precision, h.num_digits);
1913 return wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1914 dst, &h, (p > 0) ? (p - 1) : 0, options);
1915 }
1916
1917 // Use the "%f" format otherwise.
1918 int32_t p = ((int32_t)precision);
1919 if (p > h.decimal_point) {
1920 p = ((int32_t)(h.num_digits));
1921 }
1922 precision = ((uint32_t)(wuffs_base__i32__max(0, p - h.decimal_point)));
1923 return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1924 dst, &h, precision, options);
1925}