blob: 603d180817d5aa42aa25448de7a0791ee30c7430 [file] [log] [blame]
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001// After editing this file, run "go generate" in the parent directory.
2
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 }
25 uint16_t neg = ((uint16_t)(u >> 63)) << 15;
26 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;
61 uint32_t shift = 1051 - exp; // 1051 = 0x3F0 + 53 - 10.
62 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;
125 uint32_t shift = 926 - exp; // 926 = 0x380 + 53 - 23.
126 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
963// wuffs_base__private_implementation__parse_number_f64_eisel produces the IEEE
Nigel Taob15a0fc2020-07-08 10:50:14 +1000964// 754 double-precision value for an exact mantissa and base-10 exponent. For
965// example:
966// - 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//
975// The algorithm is based on an original idea by Michael Eisel. See
976// https://lemire.me/blog/2020/03/10/fast-float-parsing-in-practice/
977//
978// Preconditions:
979// - man is non-zero.
980// - exp10 is in the range -326 ..= 310, the same range of the
981// wuffs_base__private_implementation__powers_of_10 array.
982static int64_t //
983wuffs_base__private_implementation__parse_number_f64_eisel(uint64_t man,
984 int32_t exp10) {
985 // Look up the (possibly truncated) base-2 representation of (10 ** exp10).
986 // The look-up table was constructed so that it is already normalized: the
987 // table entry's mantissa's MSB (most significant bit) is on.
988 const uint32_t* po10 =
989 &wuffs_base__private_implementation__powers_of_10[5 * (exp10 + 326)];
990
991 // Normalize the man argument. The (man != 0) precondition means that a
992 // non-zero bit exists.
993 uint32_t clz = wuffs_base__count_leading_zeroes_u64(man);
994 man <<= clz;
995
996 // Calculate the return value's base-2 exponent. We might tweak it by ±1
997 // later, but its initial value comes from the look-up table and clz.
998 uint64_t ret_exp2 = ((uint64_t)po10[4]) - ((uint64_t)clz);
999
1000 // Multiply the two mantissas. Normalization means that both mantissas are at
1001 // least (1<<63), so the 128-bit product must be at least (1<<126). The high
Nigel Tao74d4af62020-07-10 11:27:17 +10001002 // 64 bits of the product, x_hi, must therefore be at least (1<<62).
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001003 //
Nigel Tao74d4af62020-07-10 11:27:17 +10001004 // As a consequence, x_hi has either 0 or 1 leading zeroes. Shifting x_hi
1005 // right by either 9 or 10 bits (depending on x_hi's MSB) will therefore
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001006 // leave the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on.
Nigel Taoe722bb02020-07-10 20:33:09 +10001007#if defined(__SIZEOF_INT128__)
Nigel Taobff01d72020-07-10 12:15:09 +10001008 // See commit 18449ad75d582dd015c236abc85a16f333b796f3 "Optimize 128-bit muls
1009 // in parse_number_f64_eisel" for benchmark numbers.
Nigel Tao18449ad2020-07-10 11:48:52 +10001010 __uint128_t x =
1011 ((__uint128_t)man) * (((uint64_t)po10[2]) | (((uint64_t)po10[3]) << 32));
1012 uint64_t x_hi = ((uint64_t)(x >> 64));
1013 uint64_t x_lo = ((uint64_t)(x));
1014#else
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001015 wuffs_base__multiply_u64__output x = wuffs_base__multiply_u64(
1016 man, ((uint64_t)po10[2]) | (((uint64_t)po10[3]) << 32));
Nigel Tao74d4af62020-07-10 11:27:17 +10001017 uint64_t x_hi = x.hi;
1018 uint64_t x_lo = x.lo;
Nigel Tao18449ad2020-07-10 11:48:52 +10001019#endif
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001020
1021 // Before we shift right by at least 9 bits, recall that the look-up table
1022 // entry was possibly truncated. We have so far only calculated a lower bound
1023 // for the product (man * e), where e is (10 ** exp10). The upper bound would
1024 // add a further (man * 1) to the 128-bit product, which overflows the lower
Nigel Tao74d4af62020-07-10 11:27:17 +10001025 // 64-bit limb if ((x_lo + man) < man).
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001026 //
Nigel Tao74d4af62020-07-10 11:27:17 +10001027 // If overflow occurs, that adds 1 to x_hi. Since we're about to shift right
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001028 // by at least 9 bits, that carried 1 can be ignored unless the higher 64-bit
1029 // limb's low 9 bits are all on.
Nigel Tao74d4af62020-07-10 11:27:17 +10001030 if (((x_hi & 0x1FF) == 0x1FF) && ((x_lo + man) < man)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001031 // Refine our calculation of (man * e). Before, our approximation of e used
1032 // a "low resolution" 64-bit mantissa. Now use a "high resolution" 128-bit
1033 // mantissa. We've already calculated x = (man * bits_0_to_63_incl_of_e).
1034 // Now calculate y = (man * bits_64_to_127_incl_of_e).
Nigel Taoe722bb02020-07-10 20:33:09 +10001035#if defined(__SIZEOF_INT128__)
Nigel Taobff01d72020-07-10 12:15:09 +10001036 // See commit 18449ad75d582dd015c236abc85a16f333b796f3 "Optimize 128-bit
1037 // muls in parse_number_f64_eisel" for benchmark numbers.
Nigel Tao18449ad2020-07-10 11:48:52 +10001038 __uint128_t y = ((__uint128_t)man) *
1039 (((uint64_t)po10[0]) | (((uint64_t)po10[1]) << 32));
1040 uint64_t y_hi = ((uint64_t)(y >> 64));
1041 uint64_t y_lo = ((uint64_t)(y));
1042#else
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001043 wuffs_base__multiply_u64__output y = wuffs_base__multiply_u64(
1044 man, ((uint64_t)po10[0]) | (((uint64_t)po10[1]) << 32));
Nigel Tao74d4af62020-07-10 11:27:17 +10001045 uint64_t y_hi = y.hi;
1046 uint64_t y_lo = y.lo;
Nigel Tao18449ad2020-07-10 11:48:52 +10001047#endif
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001048
1049 // Merge the 128-bit x and 128-bit y, which overlap by 64 bits, to
1050 // calculate the 192-bit product of the 64-bit man by the 128-bit e.
1051 // As we exit this if-block, we only care about the high 128 bits
1052 // (merged_hi and merged_lo) of that 192-bit product.
Nigel Tao74d4af62020-07-10 11:27:17 +10001053 uint64_t merged_hi = x_hi;
1054 uint64_t merged_lo = x_lo + y_hi;
1055 if (merged_lo < x_lo) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001056 merged_hi++; // Carry the overflow bit.
1057 }
1058
1059 // The "high resolution" approximation of e is still a lower bound. Once
1060 // again, see if the upper bound is large enough to produce a different
1061 // result. This time, if it does, give up instead of reaching for an even
1062 // more precise approximation to e.
1063 //
1064 // This three-part check is similar to the two-part check that guarded the
1065 // if block that we're now in, but it has an extra term for the middle 64
1066 // bits (checking that adding 1 to merged_lo would overflow).
1067 if (((merged_hi & 0x1FF) == 0x1FF) && ((merged_lo + 1) == 0) &&
Nigel Tao74d4af62020-07-10 11:27:17 +10001068 (y_lo + man < man)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001069 return -1;
1070 }
1071
1072 // Replace the 128-bit x with merged.
Nigel Tao74d4af62020-07-10 11:27:17 +10001073 x_hi = merged_hi;
1074 x_lo = merged_lo;
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001075 }
1076
Nigel Tao74d4af62020-07-10 11:27:17 +10001077 // As mentioned above, shifting x_hi right by either 9 or 10 bits will leave
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001078 // the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on. If the
1079 // MSB (before shifting) was on, adjust ret_exp2 for the larger shift.
1080 //
1081 // Having bit 53 on (and higher bits off) means that ret_mantissa is a 54-bit
1082 // number.
Nigel Tao74d4af62020-07-10 11:27:17 +10001083 uint64_t msb = x_hi >> 63;
1084 uint64_t ret_mantissa = x_hi >> (msb + 9);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001085 ret_exp2 -= 1 ^ msb;
1086
1087 // IEEE 754 rounds to-nearest with ties rounded to-even. Rounding to-even can
1088 // be tricky. If we're half-way between two exactly representable numbers
1089 // (x's low 73 bits are zero and the next 2 bits that matter are "01"), give
1090 // up instead of trying to pick the winner.
1091 //
1092 // Technically, we could tighten the condition by changing "73" to "73 or 74,
1093 // depending on msb", but a flat "73" is simpler.
Nigel Tao74d4af62020-07-10 11:27:17 +10001094 if ((x_lo == 0) && ((x_hi & 0x1FF) == 0) && ((ret_mantissa & 3) == 1)) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001095 return -1;
1096 }
1097
1098 // If we're not halfway then it's rounding to-nearest. Starting with a 54-bit
1099 // number, carry the lowest bit (bit 0) up if it's on. Regardless of whether
1100 // it was on or off, shifting right by one then produces a 53-bit number. If
1101 // carrying up overflowed, shift again.
1102 ret_mantissa += ret_mantissa & 1;
1103 ret_mantissa >>= 1;
1104 if ((ret_mantissa >> 53) > 0) {
1105 ret_mantissa >>= 1;
1106 ret_exp2++;
1107 }
1108
1109 // Starting with a 53-bit number, IEEE 754 double-precision normal numbers
1110 // have an implicit mantissa bit. Mask that away and keep the low 52 bits.
1111 ret_mantissa &= 0x000FFFFFFFFFFFFF;
1112
1113 // IEEE 754 double-precision floating point has 11 exponent bits. All off (0)
1114 // means subnormal numbers. All on (2047) means infinity or NaN.
1115 if ((ret_exp2 <= 0) || (2047 <= ret_exp2)) {
1116 return -1;
1117 }
1118
1119 // 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
1253 // Try the fast Eisel algorithm again. Calculating the (man, exp10) pair
1254 // from the high_prec_dec h is more correct but slower than the approach
1255 // taken in wuffs_base__parse_number_f64. The latter is optimized for the
1256 // common cases (e.g. assuming no underscores or a leading '+' sign) rather
1257 // than the full set of cases allowed by the Wuffs API.
1258 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));
1265 if ((man != 0) && (-326 <= exp10) && (exp10 <= 310)) {
1266 int64_t r = wuffs_base__private_implementation__parse_number_f64_eisel(
1267 man, exp10);
1268 if (r >= 0) {
1269 wuffs_base__result_f64 ret;
1270 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001271 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001272 ((uint64_t)r) | (((uint64_t)(h->negative)) << 63));
1273 return ret;
1274 }
1275 }
1276 }
1277
1278 // Scale by powers of 2 until we're in the range [½ .. 1], which gives us
1279 // our exponent (in base-2). First we shift right, possibly a little too
1280 // far, ending with a value certainly below 1 and possibly below ½...
1281 const int32_t f64_bias = -1023;
1282 int32_t exp2 = 0;
1283 while (h->decimal_point > 0) {
1284 uint32_t n = (uint32_t)(+h->decimal_point);
1285 uint32_t shift =
1286 (n < num_powers)
1287 ? powers[n]
1288 : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1289
1290 wuffs_base__private_implementation__high_prec_dec__small_rshift(h, shift);
1291 if (h->decimal_point <
1292 -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
1293 goto zero;
1294 }
1295 exp2 += (int32_t)shift;
1296 }
1297 // ...then we shift left, putting us in [½ .. 1].
1298 while (h->decimal_point <= 0) {
1299 uint32_t shift;
1300 if (h->decimal_point == 0) {
1301 if (h->digits[0] >= 5) {
1302 break;
1303 }
1304 shift = (h->digits[0] <= 2) ? 2 : 1;
1305 } else {
1306 uint32_t n = (uint32_t)(-h->decimal_point);
1307 shift = (n < num_powers)
1308 ? powers[n]
1309 : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1310 }
1311
1312 wuffs_base__private_implementation__high_prec_dec__small_lshift(h, shift);
1313 if (h->decimal_point >
1314 +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
1315 goto infinity;
1316 }
1317 exp2 -= (int32_t)shift;
1318 }
1319
1320 // We're in the range [½ .. 1] but f64 uses [1 .. 2].
1321 exp2--;
1322
1323 // The minimum normal exponent is (f64_bias + 1).
1324 while ((f64_bias + 1) > exp2) {
1325 uint32_t n = (uint32_t)((f64_bias + 1) - exp2);
1326 if (n > WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) {
1327 n = WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1328 }
1329 wuffs_base__private_implementation__high_prec_dec__small_rshift(h, n);
1330 exp2 += (int32_t)n;
1331 }
1332
1333 // Check for overflow.
1334 if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1.
1335 goto infinity;
1336 }
1337
1338 // Extract 53 bits for the mantissa (in base-2).
1339 wuffs_base__private_implementation__high_prec_dec__small_lshift(h, 53);
1340 uint64_t man2 =
1341 wuffs_base__private_implementation__high_prec_dec__rounded_integer(h);
1342
1343 // Rounding might have added one bit. If so, shift and re-check overflow.
1344 if ((man2 >> 53) != 0) {
1345 man2 >>= 1;
1346 exp2++;
1347 if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1.
1348 goto infinity;
1349 }
1350 }
1351
1352 // Handle subnormal numbers.
1353 if ((man2 >> 52) == 0) {
1354 exp2 = f64_bias;
1355 }
1356
1357 // Pack the bits and return.
1358 uint64_t exp2_bits =
1359 (uint64_t)((exp2 - f64_bias) & 0x07FF); // (1 << 11) - 1.
1360 uint64_t bits = (man2 & 0x000FFFFFFFFFFFFF) | // (1 << 52) - 1.
1361 (exp2_bits << 52) | //
1362 (h->negative ? 0x8000000000000000 : 0); // (1 << 63).
1363
1364 wuffs_base__result_f64 ret;
1365 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001366 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001367 return ret;
1368 } while (0);
1369
1370zero:
1371 do {
1372 uint64_t bits = h->negative ? 0x8000000000000000 : 0;
1373
1374 wuffs_base__result_f64 ret;
1375 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001376 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001377 return ret;
1378 } while (0);
1379
1380infinity:
1381 do {
Nigel Tao4d61a052020-07-11 12:34:40 +10001382 if (options & WUFFS_BASE__PARSE_NUMBER_FXX__REJECT_INF_AND_NAN) {
1383 wuffs_base__result_f64 ret;
1384 ret.status.repr = wuffs_base__error__bad_argument;
1385 ret.value = 0;
1386 return ret;
1387 }
1388
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001389 uint64_t bits = h->negative ? 0xFFF0000000000000 : 0x7FF0000000000000;
1390
1391 wuffs_base__result_f64 ret;
1392 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001393 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001394 return ret;
1395 } while (0);
1396}
1397
1398static inline bool //
1399wuffs_base__private_implementation__is_decimal_digit(uint8_t c) {
1400 return ('0' <= c) && (c <= '9');
1401}
1402
1403WUFFS_BASE__MAYBE_STATIC wuffs_base__result_f64 //
1404wuffs_base__parse_number_f64(wuffs_base__slice_u8 s, uint32_t options) {
1405 // In practice, almost all "dd.ddddE±xxx" numbers can be represented
1406 // losslessly by a uint64_t mantissa "dddddd" and an int32_t base-10
1407 // exponent, adjusting "xxx" for the position (if present) of the decimal
1408 // separator '.' or ','.
1409 //
1410 // This (u64 man, i32 exp10) data structure is superficially similar to the
1411 // "Do It Yourself Floating Point" type from Loitsch (†), but the exponent
1412 // here is base-10, not base-2.
1413 //
1414 // If s's number fits in a (man, exp10), parse that pair with the Eisel
1415 // algorithm. If not, or if Eisel fails, parsing s with the fallback
1416 // algorithm is slower but comprehensive.
1417 //
1418 // † "Printing Floating-Point Numbers Quickly and Accurately with Integers"
1419 // (https://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf).
1420 // Florian Loitsch is also the primary contributor to
1421 // https://github.com/google/double-conversion
1422 do {
1423 // Calculating that (man, exp10) pair needs to stay within s's bounds.
1424 // Provided that s isn't extremely long, work on a NUL-terminated copy of
1425 // s's contents. The NUL byte isn't a valid part of "±dd.ddddE±xxx".
1426 //
1427 // As the pointer p walks the contents, it's faster to repeatedly check "is
1428 // *p a valid digit" than "is p within bounds and *p a valid digit".
1429 if (s.len >= 256) {
1430 goto fallback;
1431 }
1432 uint8_t z[256];
1433 memcpy(&z[0], s.ptr, s.len);
1434 z[s.len] = 0;
1435 const uint8_t* p = &z[0];
1436
1437 // Look for a leading minus sign. Technically, we could also look for an
1438 // optional plus sign, but the "script/process-json-numbers.c with -p"
1439 // benchmark is noticably slower if we do. It's optional and, in practice,
1440 // usually absent. Let the fallback catch it.
1441 bool negative = (*p == '-');
1442 if (negative) {
1443 p++;
1444 }
1445
1446 // After walking "dd.dddd", comparing p later with p now will produce the
1447 // number of "d"s and "."s.
1448 const uint8_t* const start_of_digits_ptr = p;
1449
1450 // Walk the "d"s before a '.', 'E', NUL byte, etc. If it starts with '0',
1451 // it must be a single '0'. If it starts with a non-zero decimal digit, it
1452 // can be a sequence of decimal digits.
1453 //
1454 // Update the man variable during the walk. It's OK if man overflows now.
1455 // We'll detect that later.
1456 uint64_t man;
1457 if (*p == '0') {
1458 man = 0;
1459 p++;
1460 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1461 goto fallback;
1462 }
1463 } else if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1464 man = ((uint8_t)(*p - '0'));
1465 p++;
1466 for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) {
1467 man = (10 * man) + ((uint8_t)(*p - '0'));
1468 }
1469 } else {
1470 goto fallback;
1471 }
1472
1473 // Walk the "d"s after the optional decimal separator ('.' or ','),
1474 // updating the man and exp10 variables.
1475 int32_t exp10 = 0;
Nigel Taoe0c5de92020-07-11 11:48:17 +10001476 if (*p ==
1477 ((options & WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1478 ? ','
1479 : '.')) {
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001480 p++;
1481 const uint8_t* first_after_separator_ptr = p;
1482 if (!wuffs_base__private_implementation__is_decimal_digit(*p)) {
1483 goto fallback;
1484 }
1485 man = (10 * man) + ((uint8_t)(*p - '0'));
1486 p++;
1487 for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) {
1488 man = (10 * man) + ((uint8_t)(*p - '0'));
1489 }
1490 exp10 = ((int32_t)(first_after_separator_ptr - p));
1491 }
1492
1493 // Count the number of digits:
1494 // - for an input of "314159", digit_count is 6.
1495 // - for an input of "3.14159", digit_count is 7.
1496 //
1497 // This is off-by-one if there is a decimal separator. That's OK for now.
1498 // We'll correct for that later. The "script/process-json-numbers.c with
1499 // -p" benchmark is noticably slower if we try to correct for that now.
1500 uint32_t digit_count = (uint32_t)(p - start_of_digits_ptr);
1501
1502 // Update exp10 for the optional exponent, starting with 'E' or 'e'.
1503 if ((*p | 0x20) == 'e') {
1504 p++;
1505 int32_t exp_sign = +1;
1506 if (*p == '-') {
1507 p++;
1508 exp_sign = -1;
1509 } else if (*p == '+') {
1510 p++;
1511 }
1512 if (!wuffs_base__private_implementation__is_decimal_digit(*p)) {
1513 goto fallback;
1514 }
1515 int32_t exp_num = ((uint8_t)(*p - '0'));
1516 p++;
1517 // The rest of the exp_num walking has a peculiar control flow but, once
1518 // again, the "script/process-json-numbers.c with -p" benchmark is
1519 // sensitive to alternative formulations.
1520 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1521 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1522 p++;
1523 }
1524 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1525 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1526 p++;
1527 }
1528 while (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1529 if (exp_num > 0x1000000) {
1530 goto fallback;
1531 }
1532 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1533 p++;
1534 }
1535 exp10 += exp_sign * exp_num;
1536 }
1537
1538 // The Wuffs API is that the original slice has no trailing data. It also
1539 // allows underscores, which we don't catch here but the fallback should.
1540 if (p != &z[s.len]) {
1541 goto fallback;
1542 }
1543
1544 // Check that the uint64_t typed man variable has not overflowed, based on
1545 // digit_count.
1546 //
1547 // For reference:
1548 // - (1 << 63) is 9223372036854775808, which has 19 decimal digits.
1549 // - (1 << 64) is 18446744073709551616, which has 20 decimal digits.
1550 // - 19 nines, 9999999999999999999, is 0x8AC7230489E7FFFF, which has 64
1551 // bits and 16 hexadecimal digits.
1552 // - 20 nines, 99999999999999999999, is 0x56BC75E2D630FFFFF, which has 67
1553 // bits and 17 hexadecimal digits.
1554 if (digit_count > 19) {
1555 // Even if we have more than 19 pseudo-digits, it's not yet definitely an
1556 // overflow. Recall that digit_count might be off-by-one (too large) if
1557 // there's a decimal separator. It will also over-report the number of
1558 // meaningful digits if the input looks something like "0.000dddExxx".
1559 //
1560 // We adjust by the number of leading '0's and '.'s and re-compare to 19.
1561 // Once again, technically, we could skip ','s too, but that perturbs the
1562 // "script/process-json-numbers.c with -p" benchmark.
1563 const uint8_t* q = start_of_digits_ptr;
1564 for (; (*q == '0') || (*q == '.'); q++) {
1565 }
1566 digit_count -= (uint32_t)(q - start_of_digits_ptr);
1567 if (digit_count > 19) {
1568 goto fallback;
1569 }
1570 }
1571
1572 // The wuffs_base__private_implementation__parse_number_f64_eisel
1573 // preconditions include that exp10 is in the range -326 ..= 310.
1574 if ((exp10 < -326) || (310 < exp10)) {
1575 goto fallback;
1576 }
1577
1578 // If man and exp10 are small enough, all three of (man), (10 ** exp10) and
1579 // (man ** (10 ** exp10)) are exactly representable by a double. We don't
1580 // need to run the Eisel algorithm.
1581 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
1594 // The wuffs_base__private_implementation__parse_number_f64_eisel
1595 // preconditions include that man is non-zero. Parsing "0" should be caught
1596 // by the "If man and exp10 are small enough" above, but "0e99" might not.
1597 if (man == 0) {
1598 goto fallback;
1599 }
1600
1601 // Our man and exp10 are in range. Run the Eisel algorithm.
1602 int64_t r =
1603 wuffs_base__private_implementation__parse_number_f64_eisel(man, exp10);
1604 if (r < 0) {
1605 goto fallback;
1606 }
1607 wuffs_base__result_f64 ret;
1608 ret.status.repr = NULL;
Nigel Tao4d449dc2020-07-12 11:00:47 +10001609 ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001610 ((uint64_t)r) | (((uint64_t)negative) << 63));
1611 return ret;
1612 } while (0);
1613
1614fallback:
1615 do {
1616 wuffs_base__private_implementation__high_prec_dec h;
1617 wuffs_base__status status =
Nigel Taoe0c5de92020-07-11 11:48:17 +10001618 wuffs_base__private_implementation__high_prec_dec__parse(&h, s,
1619 options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001620 if (status.repr) {
Nigel Taoe0c5de92020-07-11 11:48:17 +10001621 return wuffs_base__private_implementation__parse_number_f64_special(
Nigel Tao4d61a052020-07-11 12:34:40 +10001622 s, options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001623 }
Nigel Tao4d61a052020-07-11 12:34:40 +10001624 return wuffs_base__private_implementation__high_prec_dec__to_f64(&h,
1625 options);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001626 } while (0);
1627}
1628
1629// --------
1630
1631static inline size_t //
1632wuffs_base__private_implementation__render_inf(wuffs_base__slice_u8 dst,
1633 bool neg,
1634 uint32_t options) {
1635 if (neg) {
1636 if (dst.len < 4) {
1637 return 0;
1638 }
1639 wuffs_base__store_u32le__no_bounds_check(dst.ptr, 0x666E492D); // '-Inf'le.
1640 return 4;
1641 }
1642
1643 if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1644 if (dst.len < 4) {
1645 return 0;
1646 }
1647 wuffs_base__store_u32le__no_bounds_check(dst.ptr, 0x666E492B); // '+Inf'le.
1648 return 4;
1649 }
1650
1651 if (dst.len < 3) {
1652 return 0;
1653 }
1654 wuffs_base__store_u24le__no_bounds_check(dst.ptr, 0x666E49); // 'Inf'le.
1655 return 3;
1656}
1657
1658static inline size_t //
1659wuffs_base__private_implementation__render_nan(wuffs_base__slice_u8 dst) {
1660 if (dst.len < 3) {
1661 return 0;
1662 }
1663 wuffs_base__store_u24le__no_bounds_check(dst.ptr, 0x4E614E); // 'NaN'le.
1664 return 3;
1665}
1666
1667static size_t //
1668wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1669 wuffs_base__slice_u8 dst,
1670 wuffs_base__private_implementation__high_prec_dec* h,
1671 uint32_t precision,
1672 uint32_t options) {
1673 size_t n = (h->negative ||
1674 (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN))
1675 ? 1
1676 : 0;
1677 if (h->decimal_point <= 0) {
1678 n += 1;
1679 } else {
1680 n += (size_t)(h->decimal_point);
1681 }
1682 if (precision > 0) {
1683 n += precision + 1; // +1 for the '.'.
1684 }
1685
1686 // Don't modify dst if the formatted number won't fit.
1687 if (n > dst.len) {
1688 return 0;
1689 }
1690
1691 // Align-left or align-right.
1692 uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT)
1693 ? &dst.ptr[dst.len - n]
1694 : &dst.ptr[0];
1695
1696 // Leading "±".
1697 if (h->negative) {
1698 *ptr++ = '-';
1699 } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1700 *ptr++ = '+';
1701 }
1702
1703 // Integral digits.
1704 if (h->decimal_point <= 0) {
1705 *ptr++ = '0';
1706 } else {
1707 uint32_t m =
1708 wuffs_base__u32__min(h->num_digits, (uint32_t)(h->decimal_point));
1709 uint32_t i = 0;
1710 for (; i < m; i++) {
1711 *ptr++ = (uint8_t)('0' | h->digits[i]);
1712 }
1713 for (; i < (uint32_t)(h->decimal_point); i++) {
1714 *ptr++ = '0';
1715 }
1716 }
1717
1718 // Separator and then fractional digits.
1719 if (precision > 0) {
1720 *ptr++ =
1721 (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1722 ? ','
1723 : '.';
1724 uint32_t i = 0;
1725 for (; i < precision; i++) {
1726 uint32_t j = ((uint32_t)(h->decimal_point)) + i;
1727 *ptr++ = (uint8_t)('0' | ((j < h->num_digits) ? h->digits[j] : 0));
1728 }
1729 }
1730
1731 return n;
1732}
1733
1734static size_t //
1735wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1736 wuffs_base__slice_u8 dst,
1737 wuffs_base__private_implementation__high_prec_dec* h,
1738 uint32_t precision,
1739 uint32_t options) {
1740 int32_t exp = 0;
1741 if (h->num_digits > 0) {
1742 exp = h->decimal_point - 1;
1743 }
1744 bool negative_exp = exp < 0;
1745 if (negative_exp) {
1746 exp = -exp;
1747 }
1748
1749 size_t n = (h->negative ||
1750 (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN))
1751 ? 4
1752 : 3; // Mininum 3 bytes: first digit and then "e±".
1753 if (precision > 0) {
1754 n += precision + 1; // +1 for the '.'.
1755 }
1756 n += (exp < 100) ? 2 : 3;
1757
1758 // Don't modify dst if the formatted number won't fit.
1759 if (n > dst.len) {
1760 return 0;
1761 }
1762
1763 // Align-left or align-right.
1764 uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT)
1765 ? &dst.ptr[dst.len - n]
1766 : &dst.ptr[0];
1767
1768 // Leading "±".
1769 if (h->negative) {
1770 *ptr++ = '-';
1771 } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1772 *ptr++ = '+';
1773 }
1774
1775 // Integral digit.
1776 if (h->num_digits > 0) {
1777 *ptr++ = (uint8_t)('0' | h->digits[0]);
1778 } else {
1779 *ptr++ = '0';
1780 }
1781
1782 // Separator and then fractional digits.
1783 if (precision > 0) {
1784 *ptr++ =
1785 (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1786 ? ','
1787 : '.';
1788 uint32_t i = 1;
1789 uint32_t j = wuffs_base__u32__min(h->num_digits, precision + 1);
1790 for (; i < j; i++) {
1791 *ptr++ = (uint8_t)('0' | h->digits[i]);
1792 }
1793 for (; i <= precision; i++) {
1794 *ptr++ = '0';
1795 }
1796 }
1797
1798 // Exponent: "e±" and then 2 or 3 digits.
1799 *ptr++ = 'e';
1800 *ptr++ = negative_exp ? '-' : '+';
1801 if (exp < 10) {
1802 *ptr++ = '0';
1803 *ptr++ = (uint8_t)('0' | exp);
1804 } else if (exp < 100) {
1805 *ptr++ = (uint8_t)('0' | (exp / 10));
1806 *ptr++ = (uint8_t)('0' | (exp % 10));
1807 } else {
1808 int32_t e = exp / 100;
1809 exp -= e * 100;
1810 *ptr++ = (uint8_t)('0' | e);
1811 *ptr++ = (uint8_t)('0' | (exp / 10));
1812 *ptr++ = (uint8_t)('0' | (exp % 10));
1813 }
1814
1815 return n;
1816}
1817
1818WUFFS_BASE__MAYBE_STATIC size_t //
1819wuffs_base__render_number_f64(wuffs_base__slice_u8 dst,
1820 double x,
1821 uint32_t precision,
1822 uint32_t options) {
1823 // Decompose x (64 bits) into negativity (1 bit), base-2 exponent (11 bits
1824 // with a -1023 bias) and mantissa (52 bits).
Nigel Tao4d449dc2020-07-12 11:00:47 +10001825 uint64_t bits = wuffs_base__ieee_754_bit_representation__from_f64_to_u64(x);
Nigel Tao2a7e1ed2020-07-07 21:50:06 +10001826 bool neg = (bits >> 63) != 0;
1827 int32_t exp2 = ((int32_t)(bits >> 52)) & 0x7FF;
1828 uint64_t man = bits & 0x000FFFFFFFFFFFFFul;
1829
1830 // Apply the exponent bias and set the implicit top bit of the mantissa,
1831 // unless x is subnormal. Also take care of Inf and NaN.
1832 if (exp2 == 0x7FF) {
1833 if (man != 0) {
1834 return wuffs_base__private_implementation__render_nan(dst);
1835 }
1836 return wuffs_base__private_implementation__render_inf(dst, neg, options);
1837 } else if (exp2 == 0) {
1838 exp2 = -1022;
1839 } else {
1840 exp2 -= 1023;
1841 man |= 0x0010000000000000ul;
1842 }
1843
1844 // Ensure that precision isn't too large.
1845 if (precision > 4095) {
1846 precision = 4095;
1847 }
1848
1849 // Convert from the (neg, exp2, man) tuple to an HPD.
1850 wuffs_base__private_implementation__high_prec_dec h;
1851 wuffs_base__private_implementation__high_prec_dec__assign(&h, man, neg);
1852 if (h.num_digits > 0) {
1853 wuffs_base__private_implementation__high_prec_dec__lshift(
1854 &h, exp2 - 52); // 52 mantissa bits.
1855 }
1856
1857 // Handle the "%e" and "%f" formats.
1858 switch (options & (WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT |
1859 WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT)) {
1860 case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT: // The "%"f" format.
1861 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1862 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1863 &h, exp2, man);
1864 int32_t p = ((int32_t)(h.num_digits)) - h.decimal_point;
1865 precision = ((uint32_t)(wuffs_base__i32__max(0, p)));
1866 } else {
1867 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1868 &h, ((int32_t)precision) + h.decimal_point);
1869 }
1870 return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1871 dst, &h, precision, options);
1872
1873 case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT: // The "%e" format.
1874 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1875 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1876 &h, exp2, man);
1877 precision = (h.num_digits > 0) ? (h.num_digits - 1) : 0;
1878 } else {
1879 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1880 &h, ((int32_t)precision) + 1);
1881 }
1882 return wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1883 dst, &h, precision, options);
1884 }
1885
1886 // We have the "%g" format and so precision means the number of significant
1887 // digits, not the number of digits after the decimal separator. Perform
1888 // rounding and determine whether to use "%e" or "%f".
1889 int32_t e_threshold = 0;
1890 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1891 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1892 &h, exp2, man);
1893 precision = h.num_digits;
1894 e_threshold = 6;
1895 } else {
1896 if (precision == 0) {
1897 precision = 1;
1898 }
1899 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1900 &h, ((int32_t)precision));
1901 e_threshold = ((int32_t)precision);
1902 int32_t nd = ((int32_t)(h.num_digits));
1903 if ((e_threshold > nd) && (nd >= h.decimal_point)) {
1904 e_threshold = nd;
1905 }
1906 }
1907
1908 // Use the "%e" format if the exponent is large.
1909 int32_t e = h.decimal_point - 1;
1910 if ((e < -4) || (e_threshold <= e)) {
1911 uint32_t p = wuffs_base__u32__min(precision, h.num_digits);
1912 return wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1913 dst, &h, (p > 0) ? (p - 1) : 0, options);
1914 }
1915
1916 // Use the "%f" format otherwise.
1917 int32_t p = ((int32_t)precision);
1918 if (p > h.decimal_point) {
1919 p = ((int32_t)(h.num_digits));
1920 }
1921 precision = ((uint32_t)(wuffs_base__i32__max(0, p - h.decimal_point)));
1922 return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1923 dst, &h, precision, options);
1924}