blob: 0369e08d7d1363700ac2f4e7b2466316818a4f0b [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
19#define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE 2047
20#define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION 800
21
22// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL is the largest N
23// such that ((10 << N) < (1 << 64)).
24#define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL 60
25
26// wuffs_base__private_implementation__high_prec_dec (abbreviated as HPD) is a
27// fixed precision floating point decimal number, augmented with ±infinity
28// values, but it cannot represent NaN (Not a Number).
29//
30// "High precision" means that the mantissa holds 800 decimal digits. 800 is
31// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION.
32//
33// An HPD isn't for general purpose arithmetic, only for conversions to and
34// from IEEE 754 double-precision floating point, where the largest and
35// smallest positive, finite values are approximately 1.8e+308 and 4.9e-324.
36// HPD exponents above +2047 mean infinity, below -2047 mean zero. The ±2047
37// bounds are further away from zero than ±(324 + 800), where 800 and 2047 is
38// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION and
39// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE.
40//
41// digits[.. num_digits] are the number's digits in big-endian order. The
42// uint8_t values are in the range [0 ..= 9], not ['0' ..= '9'], where e.g. '7'
43// is the ASCII value 0x37.
44//
45// decimal_point is the index (within digits) of the decimal point. It may be
46// negative or be larger than num_digits, in which case the explicit digits are
47// padded with implicit zeroes.
48//
49// For example, if num_digits is 3 and digits is "\x07\x08\x09":
50// - A decimal_point of -2 means ".00789"
51// - A decimal_point of -1 means ".0789"
52// - A decimal_point of +0 means ".789"
53// - A decimal_point of +1 means "7.89"
54// - A decimal_point of +2 means "78.9"
55// - A decimal_point of +3 means "789."
56// - A decimal_point of +4 means "7890."
57// - A decimal_point of +5 means "78900."
58//
59// As above, a decimal_point higher than +2047 means that the overall value is
60// infinity, lower than -2047 means zero.
61//
62// negative is a sign bit. An HPD can distinguish positive and negative zero.
63//
64// truncated is whether there are more than
65// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION digits, and at
66// least one of those extra digits are non-zero. The existence of long-tail
67// digits can affect rounding.
68//
69// The "all fields are zero" value is valid, and represents the number +0.
70typedef struct {
71 uint32_t num_digits;
72 int32_t decimal_point;
73 bool negative;
74 bool truncated;
75 uint8_t digits[WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION];
76} wuffs_base__private_implementation__high_prec_dec;
77
78// wuffs_base__private_implementation__high_prec_dec__trim trims trailing
79// zeroes from the h->digits[.. h->num_digits] slice. They have no benefit,
80// since we explicitly track h->decimal_point.
81//
82// Preconditions:
83// - h is non-NULL.
84static inline void //
85wuffs_base__private_implementation__high_prec_dec__trim(
86 wuffs_base__private_implementation__high_prec_dec* h) {
87 while ((h->num_digits > 0) && (h->digits[h->num_digits - 1] == 0)) {
88 h->num_digits--;
89 }
90}
91
92// wuffs_base__private_implementation__high_prec_dec__assign sets h to
93// represent the number x.
94//
95// Preconditions:
96// - h is non-NULL.
97static void //
98wuffs_base__private_implementation__high_prec_dec__assign(
99 wuffs_base__private_implementation__high_prec_dec* h,
100 uint64_t x,
101 bool negative) {
102 uint32_t n = 0;
103
104 // Set h->digits.
105 if (x > 0) {
106 // Calculate the digits, working right-to-left. After we determine n (how
107 // many digits there are), copy from buf to h->digits.
108 //
109 // UINT64_MAX, 18446744073709551615, is 20 digits long. It can be faster to
110 // copy a constant number of bytes than a variable number (20 instead of
111 // n). Make buf large enough (and start writing to it from the middle) so
112 // that can we always copy 20 bytes: the slice buf[(20-n) .. (40-n)].
113 uint8_t buf[40] = {0};
114 uint8_t* ptr = &buf[20];
115 do {
116 uint64_t remaining = x / 10;
117 x -= remaining * 10;
118 ptr--;
119 *ptr = (uint8_t)x;
120 n++;
121 x = remaining;
122 } while (x > 0);
123 memcpy(h->digits, ptr, 20);
124 }
125
126 // Set h's other fields.
127 h->num_digits = n;
128 h->decimal_point = (int32_t)n;
129 h->negative = negative;
130 h->truncated = false;
131 wuffs_base__private_implementation__high_prec_dec__trim(h);
132}
133
134static wuffs_base__status //
135wuffs_base__private_implementation__high_prec_dec__parse(
136 wuffs_base__private_implementation__high_prec_dec* h,
137 wuffs_base__slice_u8 s) {
138 if (!h) {
139 return wuffs_base__make_status(wuffs_base__error__bad_receiver);
140 }
141 h->num_digits = 0;
142 h->decimal_point = 0;
143 h->negative = false;
144 h->truncated = false;
145
146 uint8_t* p = s.ptr;
147 uint8_t* q = s.ptr + s.len;
148
149 for (;; p++) {
150 if (p >= q) {
151 return wuffs_base__make_status(wuffs_base__error__bad_argument);
152 } else if (*p != '_') {
153 break;
154 }
155 }
156
157 // Parse sign.
158 do {
159 if (*p == '+') {
160 p++;
161 } else if (*p == '-') {
162 h->negative = true;
163 p++;
164 } else {
165 break;
166 }
167 for (;; p++) {
168 if (p >= q) {
169 return wuffs_base__make_status(wuffs_base__error__bad_argument);
170 } else if (*p != '_') {
171 break;
172 }
173 }
174 } while (0);
175
176 // Parse digits, up to (and including) a '.', 'E' or 'e'. Examples for each
177 // limb in this if-else chain:
178 // - "0.789"
179 // - "1002.789"
180 // - ".789"
181 // - Other (invalid input).
182 uint32_t nd = 0;
183 int32_t dp = 0;
184 bool no_digits_before_separator = false;
185 if ('0' == *p) {
186 p++;
187 for (;; p++) {
188 if (p >= q) {
189 goto after_all;
190 } else if ((*p == '.') || (*p == ',')) {
191 p++;
192 goto after_sep;
193 } else if ((*p == 'E') || (*p == 'e')) {
194 p++;
195 goto after_exp;
196 } else if (*p != '_') {
197 return wuffs_base__make_status(wuffs_base__error__bad_argument);
198 }
199 }
200
201 } else if (('0' < *p) && (*p <= '9')) {
202 h->digits[nd++] = (uint8_t)(*p - '0');
203 dp = (int32_t)nd;
204 p++;
205 for (;; p++) {
206 if (p >= q) {
207 goto after_all;
208 } else if (('0' <= *p) && (*p <= '9')) {
209 if (nd < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
210 h->digits[nd++] = (uint8_t)(*p - '0');
211 dp = (int32_t)nd;
212 } else if ('0' != *p) {
213 // Long-tail non-zeroes set the truncated bit.
214 h->truncated = true;
215 }
216 } else if ((*p == '.') || (*p == ',')) {
217 p++;
218 goto after_sep;
219 } else if ((*p == 'E') || (*p == 'e')) {
220 p++;
221 goto after_exp;
222 } else if (*p != '_') {
223 return wuffs_base__make_status(wuffs_base__error__bad_argument);
224 }
225 }
226
227 } else if ((*p == '.') || (*p == ',')) {
228 p++;
229 no_digits_before_separator = true;
230
231 } else {
232 return wuffs_base__make_status(wuffs_base__error__bad_argument);
233 }
234
235after_sep:
236 for (;; p++) {
237 if (p >= q) {
238 goto after_all;
239 } else if ('0' == *p) {
240 if (nd == 0) {
241 // Track leading zeroes implicitly.
242 dp--;
243 } else if (nd <
244 WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
245 h->digits[nd++] = (uint8_t)(*p - '0');
246 }
247 } else if (('0' < *p) && (*p <= '9')) {
248 if (nd < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
249 h->digits[nd++] = (uint8_t)(*p - '0');
250 } else {
251 // Long-tail non-zeroes set the truncated bit.
252 h->truncated = true;
253 }
254 } else if ((*p == 'E') || (*p == 'e')) {
255 p++;
256 goto after_exp;
257 } else if (*p != '_') {
258 return wuffs_base__make_status(wuffs_base__error__bad_argument);
259 }
260 }
261
262after_exp:
263 do {
264 for (;; p++) {
265 if (p >= q) {
266 return wuffs_base__make_status(wuffs_base__error__bad_argument);
267 } else if (*p != '_') {
268 break;
269 }
270 }
271
272 int32_t exp_sign = +1;
273 if (*p == '+') {
274 p++;
275 } else if (*p == '-') {
276 exp_sign = -1;
277 p++;
278 }
279
280 int32_t exp = 0;
281 const int32_t exp_large =
282 WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE +
283 WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION;
284 bool saw_exp_digits = false;
285 for (; p < q; p++) {
286 if (*p == '_') {
287 // No-op.
288 } else if (('0' <= *p) && (*p <= '9')) {
289 saw_exp_digits = true;
290 if (exp < exp_large) {
291 exp = (10 * exp) + ((int32_t)(*p - '0'));
292 }
293 } else {
294 break;
295 }
296 }
297 if (!saw_exp_digits) {
298 return wuffs_base__make_status(wuffs_base__error__bad_argument);
299 }
300 dp += exp_sign * exp;
301 } while (0);
302
303after_all:
304 if (p != q) {
305 return wuffs_base__make_status(wuffs_base__error__bad_argument);
306 }
307 h->num_digits = nd;
308 if (nd == 0) {
309 if (no_digits_before_separator) {
310 return wuffs_base__make_status(wuffs_base__error__bad_argument);
311 }
312 h->decimal_point = 0;
313 } else if (dp <
314 -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
315 h->decimal_point =
316 -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE - 1;
317 } else if (dp >
318 +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
319 h->decimal_point =
320 +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE + 1;
321 } else {
322 h->decimal_point = dp;
323 }
324 wuffs_base__private_implementation__high_prec_dec__trim(h);
325 return wuffs_base__make_status(NULL);
326}
327
328// --------
329
330// wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits
331// returns the number of additional decimal digits when left-shifting by shift.
332//
333// See below for preconditions.
334static uint32_t //
335wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits(
336 wuffs_base__private_implementation__high_prec_dec* h,
337 uint32_t shift) {
338 // Masking with 0x3F should be unnecessary (assuming the preconditions) but
339 // it's cheap and ensures that we don't overflow the
340 // wuffs_base__private_implementation__hpd_left_shift array.
341 shift &= 63;
342
343 uint32_t x_a = wuffs_base__private_implementation__hpd_left_shift[shift];
344 uint32_t x_b = wuffs_base__private_implementation__hpd_left_shift[shift + 1];
345 uint32_t num_new_digits = x_a >> 11;
346 uint32_t pow5_a = 0x7FF & x_a;
347 uint32_t pow5_b = 0x7FF & x_b;
348
349 const uint8_t* pow5 =
350 &wuffs_base__private_implementation__powers_of_5[pow5_a];
351 uint32_t i = 0;
352 uint32_t n = pow5_b - pow5_a;
353 for (; i < n; i++) {
354 if (i >= h->num_digits) {
355 return num_new_digits - 1;
356 } else if (h->digits[i] == pow5[i]) {
357 continue;
358 } else if (h->digits[i] < pow5[i]) {
359 return num_new_digits - 1;
360 } else {
361 return num_new_digits;
362 }
363 }
364 return num_new_digits;
365}
366
367// --------
368
369// wuffs_base__private_implementation__high_prec_dec__rounded_integer returns
370// the integral (non-fractional) part of h, provided that it is 18 or fewer
371// decimal digits. For 19 or more digits, it returns UINT64_MAX. Note that:
372// - (1 << 53) is 9007199254740992, which has 16 decimal digits.
373// - (1 << 56) is 72057594037927936, which has 17 decimal digits.
374// - (1 << 59) is 576460752303423488, which has 18 decimal digits.
375// - (1 << 63) is 9223372036854775808, which has 19 decimal digits.
376// and that IEEE 754 double precision has 52 mantissa bits.
377//
378// That integral part is rounded-to-even: rounding 7.5 or 8.5 both give 8.
379//
380// h's negative bit is ignored: rounding -8.6 returns 9.
381//
382// See below for preconditions.
383static uint64_t //
384wuffs_base__private_implementation__high_prec_dec__rounded_integer(
385 wuffs_base__private_implementation__high_prec_dec* h) {
386 if ((h->num_digits == 0) || (h->decimal_point < 0)) {
387 return 0;
388 } else if (h->decimal_point > 18) {
389 return UINT64_MAX;
390 }
391
392 uint32_t dp = (uint32_t)(h->decimal_point);
393 uint64_t n = 0;
394 uint32_t i = 0;
395 for (; i < dp; i++) {
396 n = (10 * n) + ((i < h->num_digits) ? h->digits[i] : 0);
397 }
398
399 bool round_up = false;
400 if (dp < h->num_digits) {
401 round_up = h->digits[dp] >= 5;
402 if ((h->digits[dp] == 5) && (dp + 1 == h->num_digits)) {
403 // We are exactly halfway. If we're truncated, round up, otherwise round
404 // to even.
405 round_up = h->truncated || //
406 ((dp > 0) && (1 & h->digits[dp - 1]));
407 }
408 }
409 if (round_up) {
410 n++;
411 }
412
413 return n;
414}
415
416// wuffs_base__private_implementation__high_prec_dec__small_xshift shifts h's
417// number (where 'x' is 'l' or 'r' for left or right) by a small shift value.
418//
419// Preconditions:
420// - h is non-NULL.
421// - h->decimal_point is "not extreme".
422// - shift is non-zero.
423// - shift is "a small shift".
424//
425// "Not extreme" means within
426// ±WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE.
427//
428// "A small shift" means not more than
429// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL.
430//
431// wuffs_base__private_implementation__high_prec_dec__rounded_integer and
432// wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits
433// have the same preconditions.
434//
435// wuffs_base__private_implementation__high_prec_dec__lshift keeps the first
436// two preconditions but not the last two. Its shift argument is signed and
437// does not need to be "small": zero is a no-op, positive means left shift and
438// negative means right shift.
439
440static void //
441wuffs_base__private_implementation__high_prec_dec__small_lshift(
442 wuffs_base__private_implementation__high_prec_dec* h,
443 uint32_t shift) {
444 if (h->num_digits == 0) {
445 return;
446 }
447 uint32_t num_new_digits =
448 wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits(
449 h, shift);
450 uint32_t rx = h->num_digits - 1; // Read index.
451 uint32_t wx = h->num_digits - 1 + num_new_digits; // Write index.
452 uint64_t n = 0;
453
454 // Repeat: pick up a digit, put down a digit, right to left.
455 while (((int32_t)rx) >= 0) {
456 n += ((uint64_t)(h->digits[rx])) << shift;
457 uint64_t quo = n / 10;
458 uint64_t rem = n - (10 * quo);
459 if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
460 h->digits[wx] = (uint8_t)rem;
461 } else if (rem > 0) {
462 h->truncated = true;
463 }
464 n = quo;
465 wx--;
466 rx--;
467 }
468
469 // Put down leading digits, right to left.
470 while (n > 0) {
471 uint64_t quo = n / 10;
472 uint64_t rem = n - (10 * quo);
473 if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
474 h->digits[wx] = (uint8_t)rem;
475 } else if (rem > 0) {
476 h->truncated = true;
477 }
478 n = quo;
479 wx--;
480 }
481
482 // Finish.
483 h->num_digits += num_new_digits;
484 if (h->num_digits >
485 WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
486 h->num_digits = WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION;
487 }
488 h->decimal_point += (int32_t)num_new_digits;
489 wuffs_base__private_implementation__high_prec_dec__trim(h);
490}
491
492static void //
493wuffs_base__private_implementation__high_prec_dec__small_rshift(
494 wuffs_base__private_implementation__high_prec_dec* h,
495 uint32_t shift) {
496 uint32_t rx = 0; // Read index.
497 uint32_t wx = 0; // Write index.
498 uint64_t n = 0;
499
500 // Pick up enough leading digits to cover the first shift.
501 while ((n >> shift) == 0) {
502 if (rx < h->num_digits) {
503 // Read a digit.
504 n = (10 * n) + h->digits[rx++];
505 } else if (n == 0) {
506 // h's number used to be zero and remains zero.
507 return;
508 } else {
509 // Read sufficient implicit trailing zeroes.
510 while ((n >> shift) == 0) {
511 n = 10 * n;
512 rx++;
513 }
514 break;
515 }
516 }
517 h->decimal_point -= ((int32_t)(rx - 1));
518 if (h->decimal_point <
519 -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
520 // After the shift, h's number is effectively zero.
521 h->num_digits = 0;
522 h->decimal_point = 0;
523 h->negative = false;
524 h->truncated = false;
525 return;
526 }
527
528 // Repeat: pick up a digit, put down a digit, left to right.
529 uint64_t mask = (((uint64_t)(1)) << shift) - 1;
530 while (rx < h->num_digits) {
531 uint8_t new_digit = ((uint8_t)(n >> shift));
532 n = (10 * (n & mask)) + h->digits[rx++];
533 h->digits[wx++] = new_digit;
534 }
535
536 // Put down trailing digits, left to right.
537 while (n > 0) {
538 uint8_t new_digit = ((uint8_t)(n >> shift));
539 n = 10 * (n & mask);
540 if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
541 h->digits[wx++] = new_digit;
542 } else if (new_digit > 0) {
543 h->truncated = true;
544 }
545 }
546
547 // Finish.
548 h->num_digits = wx;
549 wuffs_base__private_implementation__high_prec_dec__trim(h);
550}
551
552static void //
553wuffs_base__private_implementation__high_prec_dec__lshift(
554 wuffs_base__private_implementation__high_prec_dec* h,
555 int32_t shift) {
556 if (shift > 0) {
557 while (shift > +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) {
558 wuffs_base__private_implementation__high_prec_dec__small_lshift(
559 h, WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL);
560 shift -= WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
561 }
562 wuffs_base__private_implementation__high_prec_dec__small_lshift(
563 h, ((uint32_t)(+shift)));
564 } else if (shift < 0) {
565 while (shift < -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) {
566 wuffs_base__private_implementation__high_prec_dec__small_rshift(
567 h, WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL);
568 shift += WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
569 }
570 wuffs_base__private_implementation__high_prec_dec__small_rshift(
571 h, ((uint32_t)(-shift)));
572 }
573}
574
575// --------
576
577// wuffs_base__private_implementation__high_prec_dec__round_etc rounds h's
578// number. For those functions that take an n argument, rounding produces at
579// most n digits (which is not necessarily at most n decimal places). Negative
580// n values are ignored, as well as any n greater than or equal to h's number
581// of digits. The etc__round_just_enough function implicitly chooses an n to
582// implement WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION.
583//
584// Preconditions:
585// - h is non-NULL.
586// - h->decimal_point is "not extreme".
587//
588// "Not extreme" means within
589// ±WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE.
590
591static void //
592wuffs_base__private_implementation__high_prec_dec__round_down(
593 wuffs_base__private_implementation__high_prec_dec* h,
594 int32_t n) {
595 if ((n < 0) || (h->num_digits <= (uint32_t)n)) {
596 return;
597 }
598 h->num_digits = (uint32_t)(n);
599 wuffs_base__private_implementation__high_prec_dec__trim(h);
600}
601
602static void //
603wuffs_base__private_implementation__high_prec_dec__round_up(
604 wuffs_base__private_implementation__high_prec_dec* h,
605 int32_t n) {
606 if ((n < 0) || (h->num_digits <= (uint32_t)n)) {
607 return;
608 }
609
610 for (n--; n >= 0; n--) {
611 if (h->digits[n] < 9) {
612 h->digits[n]++;
613 h->num_digits = (uint32_t)(n + 1);
614 return;
615 }
616 }
617
618 // The number is all 9s. Change to a single 1 and adjust the decimal point.
619 h->digits[0] = 1;
620 h->num_digits = 1;
621 h->decimal_point++;
622}
623
624static void //
625wuffs_base__private_implementation__high_prec_dec__round_nearest(
626 wuffs_base__private_implementation__high_prec_dec* h,
627 int32_t n) {
628 if ((n < 0) || (h->num_digits <= (uint32_t)n)) {
629 return;
630 }
631 bool up = h->digits[n] >= 5;
632 if ((h->digits[n] == 5) && ((n + 1) == ((int32_t)(h->num_digits)))) {
633 up = h->truncated || //
634 ((n > 0) && ((h->digits[n - 1] & 1) != 0));
635 }
636
637 if (up) {
638 wuffs_base__private_implementation__high_prec_dec__round_up(h, n);
639 } else {
640 wuffs_base__private_implementation__high_prec_dec__round_down(h, n);
641 }
642}
643
644static void //
645wuffs_base__private_implementation__high_prec_dec__round_just_enough(
646 wuffs_base__private_implementation__high_prec_dec* h,
647 int32_t exp2,
648 uint64_t mantissa) {
649 // The magic numbers 52 and 53 in this function are because IEEE 754 double
650 // precision has 52 mantissa bits.
651 //
652 // Let f be the floating point number represented by exp2 and mantissa (and
653 // also the number in h): the number (mantissa * (2 ** (exp2 - 52))).
654 //
655 // If f is zero or a small integer, we can return early.
656 if ((mantissa == 0) ||
657 ((exp2 < 53) && (h->decimal_point >= ((int32_t)(h->num_digits))))) {
658 return;
659 }
660
661 // The smallest normal f has an exp2 of -1022 and a mantissa of (1 << 52).
662 // Subnormal numbers have the same exp2 but a smaller mantissa.
663 static const int32_t min_incl_normal_exp2 = -1022;
664 static const uint64_t min_incl_normal_mantissa = 0x0010000000000000ul;
665
666 // Compute lower and upper bounds such that any number between them (possibly
667 // inclusive) will round to f. First, the lower bound. Our number f is:
668 // ((mantissa + 0) * (2 ** ( exp2 - 52)))
669 //
670 // The next lowest floating point number is:
671 // ((mantissa - 1) * (2 ** ( exp2 - 52)))
672 // unless (mantissa - 1) drops the (1 << 52) bit and exp2 is not the
673 // min_incl_normal_exp2. Either way, call it:
674 // ((l_mantissa) * (2 ** (l_exp2 - 52)))
675 //
676 // The lower bound is halfway between them (noting that 52 became 53):
677 // (((2 * l_mantissa) + 1) * (2 ** (l_exp2 - 53)))
678 int32_t l_exp2 = exp2;
679 uint64_t l_mantissa = mantissa - 1;
680 if ((exp2 > min_incl_normal_exp2) && (mantissa <= min_incl_normal_mantissa)) {
681 l_exp2 = exp2 - 1;
682 l_mantissa = (2 * mantissa) - 1;
683 }
684 wuffs_base__private_implementation__high_prec_dec lower;
685 wuffs_base__private_implementation__high_prec_dec__assign(
686 &lower, (2 * l_mantissa) + 1, false);
687 wuffs_base__private_implementation__high_prec_dec__lshift(&lower,
688 l_exp2 - 53);
689
690 // Next, the upper bound. Our number f is:
691 // ((mantissa + 0) * (2 ** (exp2 - 52)))
692 //
693 // The next highest floating point number is:
694 // ((mantissa + 1) * (2 ** (exp2 - 52)))
695 //
696 // The upper bound is halfway between them (noting that 52 became 53):
697 // (((2 * mantissa) + 1) * (2 ** (exp2 - 53)))
698 wuffs_base__private_implementation__high_prec_dec upper;
699 wuffs_base__private_implementation__high_prec_dec__assign(
700 &upper, (2 * mantissa) + 1, false);
701 wuffs_base__private_implementation__high_prec_dec__lshift(&upper, exp2 - 53);
702
703 // The lower and upper bounds are possible outputs only if the original
704 // mantissa is even, so that IEEE round-to-even would round to the original
705 // mantissa and not its neighbors.
706 bool inclusive = (mantissa & 1) == 0;
707
708 // As we walk the digits, we want to know whether rounding up would fall
709 // within the upper bound. This is tracked by upper_delta:
710 // - When -1, the digits of h and upper are the same so far.
711 // - When +0, we saw a difference of 1 between h and upper on a previous
712 // digit and subsequently only 9s for h and 0s for upper. Thus, rounding
713 // up may fall outside of the bound if !inclusive.
714 // - When +1, the difference is greater than 1 and we know that rounding up
715 // falls within the bound.
716 //
717 // This is a state machine with three states. The numerical value for each
718 // state (-1, +0 or +1) isn't important, other than their order.
719 int upper_delta = -1;
720
721 // We can now figure out the shortest number of digits required. Walk the
722 // digits until h has distinguished itself from lower or upper.
723 //
724 // The zi and zd variables are indexes and digits, for z in l (lower), h (the
725 // number) and u (upper).
726 //
727 // The lower, h and upper numbers may have their decimal points at different
728 // places. In this case, upper is the longest, so we iterate ui starting from
729 // 0 and iterate li and hi starting from either 0 or -1.
730 int32_t ui = 0;
731 for (;; ui++) {
732 // Calculate hd, the middle number's digit.
733 int32_t hi = ui - upper.decimal_point + h->decimal_point;
734 if (hi >= ((int32_t)(h->num_digits))) {
735 break;
736 }
737 uint8_t hd = (((uint32_t)hi) < h->num_digits) ? h->digits[hi] : 0;
738
739 // Calculate ld, the lower bound's digit.
740 int32_t li = ui - upper.decimal_point + lower.decimal_point;
741 uint8_t ld = (((uint32_t)li) < lower.num_digits) ? lower.digits[li] : 0;
742
743 // We can round down (truncate) if lower has a different digit than h or if
744 // lower is inclusive and is exactly the result of rounding down (i.e. we
745 // have reached the final digit of lower).
746 bool can_round_down =
747 (ld != hd) || //
748 (inclusive && ((li + 1) == ((int32_t)(lower.num_digits))));
749
750 // Calculate ud, the upper bound's digit, and update upper_delta.
751 uint8_t ud = (((uint32_t)ui) < upper.num_digits) ? upper.digits[ui] : 0;
752 if (upper_delta < 0) {
753 if ((hd + 1) < ud) {
754 // For example:
755 // h = 12345???
756 // upper = 12347???
757 upper_delta = +1;
758 } else if (hd != ud) {
759 // For example:
760 // h = 12345???
761 // upper = 12346???
762 upper_delta = +0;
763 }
764 } else if (upper_delta == 0) {
765 if ((hd != 9) || (ud != 0)) {
766 // For example:
767 // h = 1234598?
768 // upper = 1234600?
769 upper_delta = +1;
770 }
771 }
772
773 // We can round up if upper has a different digit than h and either upper
774 // is inclusive or upper is bigger than the result of rounding up.
775 bool can_round_up =
776 (upper_delta > 0) || //
777 ((upper_delta == 0) && //
778 (inclusive || ((ui + 1) < ((int32_t)(upper.num_digits)))));
779
780 // If we can round either way, round to nearest. If we can round only one
781 // way, do it. If we can't round, continue the loop.
782 if (can_round_down) {
783 if (can_round_up) {
784 wuffs_base__private_implementation__high_prec_dec__round_nearest(
785 h, hi + 1);
786 return;
787 } else {
788 wuffs_base__private_implementation__high_prec_dec__round_down(h,
789 hi + 1);
790 return;
791 }
792 } else {
793 if (can_round_up) {
794 wuffs_base__private_implementation__high_prec_dec__round_up(h, hi + 1);
795 return;
796 }
797 }
798 }
799}
800
801// --------
802
803// wuffs_base__private_implementation__parse_number_f64_eisel produces the IEEE
Nigel Taob15a0fc2020-07-08 10:50:14 +1000804// 754 double-precision value for an exact mantissa and base-10 exponent. For
805// example:
806// - when parsing "12345.678e+02", man is 12345678 and exp10 is -1.
807// - when parsing "-12", man is 12 and exp10 is 0. Processing the leading
808// minus sign is the responsibility of the caller, not this function.
Nigel Tao2a7e1ed2020-07-07 21:50:06 +1000809//
810// On success, it returns a non-negative int64_t such that the low 63 bits hold
811// the 11-bit exponent and 52-bit mantissa.
812//
813// On failure, it returns a negative value.
814//
815// The algorithm is based on an original idea by Michael Eisel. See
816// https://lemire.me/blog/2020/03/10/fast-float-parsing-in-practice/
817//
818// Preconditions:
819// - man is non-zero.
820// - exp10 is in the range -326 ..= 310, the same range of the
821// wuffs_base__private_implementation__powers_of_10 array.
822static int64_t //
823wuffs_base__private_implementation__parse_number_f64_eisel(uint64_t man,
824 int32_t exp10) {
825 // Look up the (possibly truncated) base-2 representation of (10 ** exp10).
826 // The look-up table was constructed so that it is already normalized: the
827 // table entry's mantissa's MSB (most significant bit) is on.
828 const uint32_t* po10 =
829 &wuffs_base__private_implementation__powers_of_10[5 * (exp10 + 326)];
830
831 // Normalize the man argument. The (man != 0) precondition means that a
832 // non-zero bit exists.
833 uint32_t clz = wuffs_base__count_leading_zeroes_u64(man);
834 man <<= clz;
835
836 // Calculate the return value's base-2 exponent. We might tweak it by ±1
837 // later, but its initial value comes from the look-up table and clz.
838 uint64_t ret_exp2 = ((uint64_t)po10[4]) - ((uint64_t)clz);
839
840 // Multiply the two mantissas. Normalization means that both mantissas are at
841 // least (1<<63), so the 128-bit product must be at least (1<<126). The high
842 // 64 bits of the product, x.hi, must therefore be at least (1<<62).
843 //
844 // As a consequence, x.hi has either 0 or 1 leading zeroes. Shifting x.hi
845 // right by either 9 or 10 bits (depending on x.hi's MSB) will therefore
846 // leave the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on.
847 wuffs_base__multiply_u64__output x = wuffs_base__multiply_u64(
848 man, ((uint64_t)po10[2]) | (((uint64_t)po10[3]) << 32));
849
850 // Before we shift right by at least 9 bits, recall that the look-up table
851 // entry was possibly truncated. We have so far only calculated a lower bound
852 // for the product (man * e), where e is (10 ** exp10). The upper bound would
853 // add a further (man * 1) to the 128-bit product, which overflows the lower
854 // 64-bit limb if ((x.lo + man) < man).
855 //
856 // If overflow occurs, that adds 1 to x.hi. Since we're about to shift right
857 // by at least 9 bits, that carried 1 can be ignored unless the higher 64-bit
858 // limb's low 9 bits are all on.
859 if (((x.hi & 0x1FF) == 0x1FF) && ((x.lo + man) < man)) {
860 // Refine our calculation of (man * e). Before, our approximation of e used
861 // a "low resolution" 64-bit mantissa. Now use a "high resolution" 128-bit
862 // mantissa. We've already calculated x = (man * bits_0_to_63_incl_of_e).
863 // Now calculate y = (man * bits_64_to_127_incl_of_e).
864 wuffs_base__multiply_u64__output y = wuffs_base__multiply_u64(
865 man, ((uint64_t)po10[0]) | (((uint64_t)po10[1]) << 32));
866
867 // Merge the 128-bit x and 128-bit y, which overlap by 64 bits, to
868 // calculate the 192-bit product of the 64-bit man by the 128-bit e.
869 // As we exit this if-block, we only care about the high 128 bits
870 // (merged_hi and merged_lo) of that 192-bit product.
871 uint64_t merged_hi = x.hi;
872 uint64_t merged_lo = x.lo + y.hi;
873 if (merged_lo < x.lo) {
874 merged_hi++; // Carry the overflow bit.
875 }
876
877 // The "high resolution" approximation of e is still a lower bound. Once
878 // again, see if the upper bound is large enough to produce a different
879 // result. This time, if it does, give up instead of reaching for an even
880 // more precise approximation to e.
881 //
882 // This three-part check is similar to the two-part check that guarded the
883 // if block that we're now in, but it has an extra term for the middle 64
884 // bits (checking that adding 1 to merged_lo would overflow).
885 if (((merged_hi & 0x1FF) == 0x1FF) && ((merged_lo + 1) == 0) &&
886 (y.lo + man < man)) {
887 return -1;
888 }
889
890 // Replace the 128-bit x with merged.
891 x.hi = merged_hi;
892 x.lo = merged_lo;
893 }
894
895 // As mentioned above, shifting x.hi right by either 9 or 10 bits will leave
896 // the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on. If the
897 // MSB (before shifting) was on, adjust ret_exp2 for the larger shift.
898 //
899 // Having bit 53 on (and higher bits off) means that ret_mantissa is a 54-bit
900 // number.
901 uint64_t msb = x.hi >> 63;
902 uint64_t ret_mantissa = x.hi >> (msb + 9);
903 ret_exp2 -= 1 ^ msb;
904
905 // IEEE 754 rounds to-nearest with ties rounded to-even. Rounding to-even can
906 // be tricky. If we're half-way between two exactly representable numbers
907 // (x's low 73 bits are zero and the next 2 bits that matter are "01"), give
908 // up instead of trying to pick the winner.
909 //
910 // Technically, we could tighten the condition by changing "73" to "73 or 74,
911 // depending on msb", but a flat "73" is simpler.
912 if ((x.lo == 0) && ((x.hi & 0x1FF) == 0) && ((ret_mantissa & 3) == 1)) {
913 return -1;
914 }
915
916 // If we're not halfway then it's rounding to-nearest. Starting with a 54-bit
917 // number, carry the lowest bit (bit 0) up if it's on. Regardless of whether
918 // it was on or off, shifting right by one then produces a 53-bit number. If
919 // carrying up overflowed, shift again.
920 ret_mantissa += ret_mantissa & 1;
921 ret_mantissa >>= 1;
922 if ((ret_mantissa >> 53) > 0) {
923 ret_mantissa >>= 1;
924 ret_exp2++;
925 }
926
927 // Starting with a 53-bit number, IEEE 754 double-precision normal numbers
928 // have an implicit mantissa bit. Mask that away and keep the low 52 bits.
929 ret_mantissa &= 0x000FFFFFFFFFFFFF;
930
931 // IEEE 754 double-precision floating point has 11 exponent bits. All off (0)
932 // means subnormal numbers. All on (2047) means infinity or NaN.
933 if ((ret_exp2 <= 0) || (2047 <= ret_exp2)) {
934 return -1;
935 }
936
937 // Pack the bits and return.
938 return ((int64_t)(ret_mantissa | (ret_exp2 << 52)));
939}
940
941// --------
942
943static wuffs_base__result_f64 //
944wuffs_base__parse_number_f64_special(wuffs_base__slice_u8 s,
945 const char* fallback_status_repr) {
946 do {
947 uint8_t* p = s.ptr;
948 uint8_t* q = s.ptr + s.len;
949
950 for (; (p < q) && (*p == '_'); p++) {
951 }
952 if (p >= q) {
953 goto fallback;
954 }
955
956 // Parse sign.
957 bool negative = false;
958 do {
959 if (*p == '+') {
960 p++;
961 } else if (*p == '-') {
962 negative = true;
963 p++;
964 } else {
965 break;
966 }
967 for (; (p < q) && (*p == '_'); p++) {
968 }
969 } while (0);
970 if (p >= q) {
971 goto fallback;
972 }
973
974 bool nan = false;
975 switch (p[0]) {
976 case 'I':
977 case 'i':
978 if (((q - p) < 3) || //
979 ((p[1] != 'N') && (p[1] != 'n')) || //
980 ((p[2] != 'F') && (p[2] != 'f'))) {
981 goto fallback;
982 }
983 p += 3;
984
985 if ((p >= q) || (*p == '_')) {
986 break;
987 } else if (((q - p) < 5) || //
988 ((p[0] != 'I') && (p[0] != 'i')) || //
989 ((p[1] != 'N') && (p[1] != 'n')) || //
990 ((p[2] != 'I') && (p[2] != 'i')) || //
991 ((p[3] != 'T') && (p[3] != 't')) || //
992 ((p[4] != 'Y') && (p[4] != 'y'))) {
993 goto fallback;
994 }
995 p += 5;
996
997 if ((p >= q) || (*p == '_')) {
998 break;
999 }
1000 goto fallback;
1001
1002 case 'N':
1003 case 'n':
1004 if (((q - p) < 3) || //
1005 ((p[1] != 'A') && (p[1] != 'a')) || //
1006 ((p[2] != 'N') && (p[2] != 'n'))) {
1007 goto fallback;
1008 }
1009 p += 3;
1010
1011 if ((p >= q) || (*p == '_')) {
1012 nan = true;
1013 break;
1014 }
1015 goto fallback;
1016
1017 default:
1018 goto fallback;
1019 }
1020
1021 // Finish.
1022 for (; (p < q) && (*p == '_'); p++) {
1023 }
1024 if (p != q) {
1025 goto fallback;
1026 }
1027 wuffs_base__result_f64 ret;
1028 ret.status.repr = NULL;
1029 ret.value = wuffs_base__ieee_754_bit_representation__to_f64(
1030 (nan ? 0x7FFFFFFFFFFFFFFF : 0x7FF0000000000000) |
1031 (negative ? 0x8000000000000000 : 0));
1032 return ret;
1033 } while (0);
1034
1035fallback:
1036 do {
1037 wuffs_base__result_f64 ret;
1038 ret.status.repr = fallback_status_repr;
1039 ret.value = 0;
1040 return ret;
1041 } while (0);
1042}
1043
1044WUFFS_BASE__MAYBE_STATIC wuffs_base__result_f64 //
1045wuffs_base__private_implementation__parse_number_f64__fallback(
1046 wuffs_base__private_implementation__high_prec_dec* h) {
1047 do {
1048 // powers converts decimal powers of 10 to binary powers of 2. For example,
1049 // (10000 >> 13) is 1. It stops before the elements exceed 60, also known
1050 // as WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL.
1051 static const uint32_t num_powers = 19;
1052 static const uint8_t powers[19] = {
1053 0, 3, 6, 9, 13, 16, 19, 23, 26, 29, //
1054 33, 36, 39, 43, 46, 49, 53, 56, 59, //
1055 };
1056
1057 // Handle zero and obvious extremes. The largest and smallest positive
1058 // finite f64 values are approximately 1.8e+308 and 4.9e-324.
1059 if ((h->num_digits == 0) || (h->decimal_point < -326)) {
1060 goto zero;
1061 } else if (h->decimal_point > 310) {
1062 goto infinity;
1063 }
1064
1065 // Try the fast Eisel algorithm again. Calculating the (man, exp10) pair
1066 // from the high_prec_dec h is more correct but slower than the approach
1067 // taken in wuffs_base__parse_number_f64. The latter is optimized for the
1068 // common cases (e.g. assuming no underscores or a leading '+' sign) rather
1069 // than the full set of cases allowed by the Wuffs API.
1070 if (h->num_digits <= 19) {
1071 uint64_t man = 0;
1072 uint32_t i;
1073 for (i = 0; i < h->num_digits; i++) {
1074 man = (10 * man) + h->digits[i];
1075 }
1076 int32_t exp10 = h->decimal_point - ((int32_t)(h->num_digits));
1077 if ((man != 0) && (-326 <= exp10) && (exp10 <= 310)) {
1078 int64_t r = wuffs_base__private_implementation__parse_number_f64_eisel(
1079 man, exp10);
1080 if (r >= 0) {
1081 wuffs_base__result_f64 ret;
1082 ret.status.repr = NULL;
1083 ret.value = wuffs_base__ieee_754_bit_representation__to_f64(
1084 ((uint64_t)r) | (((uint64_t)(h->negative)) << 63));
1085 return ret;
1086 }
1087 }
1088 }
1089
1090 // Scale by powers of 2 until we're in the range [½ .. 1], which gives us
1091 // our exponent (in base-2). First we shift right, possibly a little too
1092 // far, ending with a value certainly below 1 and possibly below ½...
1093 const int32_t f64_bias = -1023;
1094 int32_t exp2 = 0;
1095 while (h->decimal_point > 0) {
1096 uint32_t n = (uint32_t)(+h->decimal_point);
1097 uint32_t shift =
1098 (n < num_powers)
1099 ? powers[n]
1100 : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1101
1102 wuffs_base__private_implementation__high_prec_dec__small_rshift(h, shift);
1103 if (h->decimal_point <
1104 -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
1105 goto zero;
1106 }
1107 exp2 += (int32_t)shift;
1108 }
1109 // ...then we shift left, putting us in [½ .. 1].
1110 while (h->decimal_point <= 0) {
1111 uint32_t shift;
1112 if (h->decimal_point == 0) {
1113 if (h->digits[0] >= 5) {
1114 break;
1115 }
1116 shift = (h->digits[0] <= 2) ? 2 : 1;
1117 } else {
1118 uint32_t n = (uint32_t)(-h->decimal_point);
1119 shift = (n < num_powers)
1120 ? powers[n]
1121 : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1122 }
1123
1124 wuffs_base__private_implementation__high_prec_dec__small_lshift(h, shift);
1125 if (h->decimal_point >
1126 +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
1127 goto infinity;
1128 }
1129 exp2 -= (int32_t)shift;
1130 }
1131
1132 // We're in the range [½ .. 1] but f64 uses [1 .. 2].
1133 exp2--;
1134
1135 // The minimum normal exponent is (f64_bias + 1).
1136 while ((f64_bias + 1) > exp2) {
1137 uint32_t n = (uint32_t)((f64_bias + 1) - exp2);
1138 if (n > WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) {
1139 n = WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
1140 }
1141 wuffs_base__private_implementation__high_prec_dec__small_rshift(h, n);
1142 exp2 += (int32_t)n;
1143 }
1144
1145 // Check for overflow.
1146 if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1.
1147 goto infinity;
1148 }
1149
1150 // Extract 53 bits for the mantissa (in base-2).
1151 wuffs_base__private_implementation__high_prec_dec__small_lshift(h, 53);
1152 uint64_t man2 =
1153 wuffs_base__private_implementation__high_prec_dec__rounded_integer(h);
1154
1155 // Rounding might have added one bit. If so, shift and re-check overflow.
1156 if ((man2 >> 53) != 0) {
1157 man2 >>= 1;
1158 exp2++;
1159 if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1.
1160 goto infinity;
1161 }
1162 }
1163
1164 // Handle subnormal numbers.
1165 if ((man2 >> 52) == 0) {
1166 exp2 = f64_bias;
1167 }
1168
1169 // Pack the bits and return.
1170 uint64_t exp2_bits =
1171 (uint64_t)((exp2 - f64_bias) & 0x07FF); // (1 << 11) - 1.
1172 uint64_t bits = (man2 & 0x000FFFFFFFFFFFFF) | // (1 << 52) - 1.
1173 (exp2_bits << 52) | //
1174 (h->negative ? 0x8000000000000000 : 0); // (1 << 63).
1175
1176 wuffs_base__result_f64 ret;
1177 ret.status.repr = NULL;
1178 ret.value = wuffs_base__ieee_754_bit_representation__to_f64(bits);
1179 return ret;
1180 } while (0);
1181
1182zero:
1183 do {
1184 uint64_t bits = h->negative ? 0x8000000000000000 : 0;
1185
1186 wuffs_base__result_f64 ret;
1187 ret.status.repr = NULL;
1188 ret.value = wuffs_base__ieee_754_bit_representation__to_f64(bits);
1189 return ret;
1190 } while (0);
1191
1192infinity:
1193 do {
1194 uint64_t bits = h->negative ? 0xFFF0000000000000 : 0x7FF0000000000000;
1195
1196 wuffs_base__result_f64 ret;
1197 ret.status.repr = NULL;
1198 ret.value = wuffs_base__ieee_754_bit_representation__to_f64(bits);
1199 return ret;
1200 } while (0);
1201}
1202
1203static inline bool //
1204wuffs_base__private_implementation__is_decimal_digit(uint8_t c) {
1205 return ('0' <= c) && (c <= '9');
1206}
1207
1208WUFFS_BASE__MAYBE_STATIC wuffs_base__result_f64 //
1209wuffs_base__parse_number_f64(wuffs_base__slice_u8 s, uint32_t options) {
1210 // In practice, almost all "dd.ddddE±xxx" numbers can be represented
1211 // losslessly by a uint64_t mantissa "dddddd" and an int32_t base-10
1212 // exponent, adjusting "xxx" for the position (if present) of the decimal
1213 // separator '.' or ','.
1214 //
1215 // This (u64 man, i32 exp10) data structure is superficially similar to the
1216 // "Do It Yourself Floating Point" type from Loitsch (†), but the exponent
1217 // here is base-10, not base-2.
1218 //
1219 // If s's number fits in a (man, exp10), parse that pair with the Eisel
1220 // algorithm. If not, or if Eisel fails, parsing s with the fallback
1221 // algorithm is slower but comprehensive.
1222 //
1223 // † "Printing Floating-Point Numbers Quickly and Accurately with Integers"
1224 // (https://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf).
1225 // Florian Loitsch is also the primary contributor to
1226 // https://github.com/google/double-conversion
1227 do {
1228 // Calculating that (man, exp10) pair needs to stay within s's bounds.
1229 // Provided that s isn't extremely long, work on a NUL-terminated copy of
1230 // s's contents. The NUL byte isn't a valid part of "±dd.ddddE±xxx".
1231 //
1232 // As the pointer p walks the contents, it's faster to repeatedly check "is
1233 // *p a valid digit" than "is p within bounds and *p a valid digit".
1234 if (s.len >= 256) {
1235 goto fallback;
1236 }
1237 uint8_t z[256];
1238 memcpy(&z[0], s.ptr, s.len);
1239 z[s.len] = 0;
1240 const uint8_t* p = &z[0];
1241
1242 // Look for a leading minus sign. Technically, we could also look for an
1243 // optional plus sign, but the "script/process-json-numbers.c with -p"
1244 // benchmark is noticably slower if we do. It's optional and, in practice,
1245 // usually absent. Let the fallback catch it.
1246 bool negative = (*p == '-');
1247 if (negative) {
1248 p++;
1249 }
1250
1251 // After walking "dd.dddd", comparing p later with p now will produce the
1252 // number of "d"s and "."s.
1253 const uint8_t* const start_of_digits_ptr = p;
1254
1255 // Walk the "d"s before a '.', 'E', NUL byte, etc. If it starts with '0',
1256 // it must be a single '0'. If it starts with a non-zero decimal digit, it
1257 // can be a sequence of decimal digits.
1258 //
1259 // Update the man variable during the walk. It's OK if man overflows now.
1260 // We'll detect that later.
1261 uint64_t man;
1262 if (*p == '0') {
1263 man = 0;
1264 p++;
1265 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1266 goto fallback;
1267 }
1268 } else if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1269 man = ((uint8_t)(*p - '0'));
1270 p++;
1271 for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) {
1272 man = (10 * man) + ((uint8_t)(*p - '0'));
1273 }
1274 } else {
1275 goto fallback;
1276 }
1277
1278 // Walk the "d"s after the optional decimal separator ('.' or ','),
1279 // updating the man and exp10 variables.
1280 int32_t exp10 = 0;
1281 if ((*p == '.') || (*p == ',')) {
1282 p++;
1283 const uint8_t* first_after_separator_ptr = p;
1284 if (!wuffs_base__private_implementation__is_decimal_digit(*p)) {
1285 goto fallback;
1286 }
1287 man = (10 * man) + ((uint8_t)(*p - '0'));
1288 p++;
1289 for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) {
1290 man = (10 * man) + ((uint8_t)(*p - '0'));
1291 }
1292 exp10 = ((int32_t)(first_after_separator_ptr - p));
1293 }
1294
1295 // Count the number of digits:
1296 // - for an input of "314159", digit_count is 6.
1297 // - for an input of "3.14159", digit_count is 7.
1298 //
1299 // This is off-by-one if there is a decimal separator. That's OK for now.
1300 // We'll correct for that later. The "script/process-json-numbers.c with
1301 // -p" benchmark is noticably slower if we try to correct for that now.
1302 uint32_t digit_count = (uint32_t)(p - start_of_digits_ptr);
1303
1304 // Update exp10 for the optional exponent, starting with 'E' or 'e'.
1305 if ((*p | 0x20) == 'e') {
1306 p++;
1307 int32_t exp_sign = +1;
1308 if (*p == '-') {
1309 p++;
1310 exp_sign = -1;
1311 } else if (*p == '+') {
1312 p++;
1313 }
1314 if (!wuffs_base__private_implementation__is_decimal_digit(*p)) {
1315 goto fallback;
1316 }
1317 int32_t exp_num = ((uint8_t)(*p - '0'));
1318 p++;
1319 // The rest of the exp_num walking has a peculiar control flow but, once
1320 // again, the "script/process-json-numbers.c with -p" benchmark is
1321 // sensitive to alternative formulations.
1322 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1323 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1324 p++;
1325 }
1326 if (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1327 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1328 p++;
1329 }
1330 while (wuffs_base__private_implementation__is_decimal_digit(*p)) {
1331 if (exp_num > 0x1000000) {
1332 goto fallback;
1333 }
1334 exp_num = (10 * exp_num) + ((uint8_t)(*p - '0'));
1335 p++;
1336 }
1337 exp10 += exp_sign * exp_num;
1338 }
1339
1340 // The Wuffs API is that the original slice has no trailing data. It also
1341 // allows underscores, which we don't catch here but the fallback should.
1342 if (p != &z[s.len]) {
1343 goto fallback;
1344 }
1345
1346 // Check that the uint64_t typed man variable has not overflowed, based on
1347 // digit_count.
1348 //
1349 // For reference:
1350 // - (1 << 63) is 9223372036854775808, which has 19 decimal digits.
1351 // - (1 << 64) is 18446744073709551616, which has 20 decimal digits.
1352 // - 19 nines, 9999999999999999999, is 0x8AC7230489E7FFFF, which has 64
1353 // bits and 16 hexadecimal digits.
1354 // - 20 nines, 99999999999999999999, is 0x56BC75E2D630FFFFF, which has 67
1355 // bits and 17 hexadecimal digits.
1356 if (digit_count > 19) {
1357 // Even if we have more than 19 pseudo-digits, it's not yet definitely an
1358 // overflow. Recall that digit_count might be off-by-one (too large) if
1359 // there's a decimal separator. It will also over-report the number of
1360 // meaningful digits if the input looks something like "0.000dddExxx".
1361 //
1362 // We adjust by the number of leading '0's and '.'s and re-compare to 19.
1363 // Once again, technically, we could skip ','s too, but that perturbs the
1364 // "script/process-json-numbers.c with -p" benchmark.
1365 const uint8_t* q = start_of_digits_ptr;
1366 for (; (*q == '0') || (*q == '.'); q++) {
1367 }
1368 digit_count -= (uint32_t)(q - start_of_digits_ptr);
1369 if (digit_count > 19) {
1370 goto fallback;
1371 }
1372 }
1373
1374 // The wuffs_base__private_implementation__parse_number_f64_eisel
1375 // preconditions include that exp10 is in the range -326 ..= 310.
1376 if ((exp10 < -326) || (310 < exp10)) {
1377 goto fallback;
1378 }
1379
1380 // If man and exp10 are small enough, all three of (man), (10 ** exp10) and
1381 // (man ** (10 ** exp10)) are exactly representable by a double. We don't
1382 // need to run the Eisel algorithm.
1383 if ((-22 <= exp10) && (exp10 <= 22) && ((man >> 53) == 0)) {
1384 double d = (double)man;
1385 if (exp10 >= 0) {
1386 d *= wuffs_base__private_implementation__f64_powers_of_10[+exp10];
1387 } else {
1388 d /= wuffs_base__private_implementation__f64_powers_of_10[-exp10];
1389 }
1390 wuffs_base__result_f64 ret;
1391 ret.status.repr = NULL;
1392 ret.value = negative ? -d : +d;
1393 return ret;
1394 }
1395
1396 // The wuffs_base__private_implementation__parse_number_f64_eisel
1397 // preconditions include that man is non-zero. Parsing "0" should be caught
1398 // by the "If man and exp10 are small enough" above, but "0e99" might not.
1399 if (man == 0) {
1400 goto fallback;
1401 }
1402
1403 // Our man and exp10 are in range. Run the Eisel algorithm.
1404 int64_t r =
1405 wuffs_base__private_implementation__parse_number_f64_eisel(man, exp10);
1406 if (r < 0) {
1407 goto fallback;
1408 }
1409 wuffs_base__result_f64 ret;
1410 ret.status.repr = NULL;
1411 ret.value = wuffs_base__ieee_754_bit_representation__to_f64(
1412 ((uint64_t)r) | (((uint64_t)negative) << 63));
1413 return ret;
1414 } while (0);
1415
1416fallback:
1417 do {
1418 wuffs_base__private_implementation__high_prec_dec h;
1419 wuffs_base__status status =
1420 wuffs_base__private_implementation__high_prec_dec__parse(&h, s);
1421 if (status.repr) {
1422 return wuffs_base__parse_number_f64_special(s, status.repr);
1423 }
1424 return wuffs_base__private_implementation__parse_number_f64__fallback(&h);
1425 } while (0);
1426}
1427
1428// --------
1429
1430static inline size_t //
1431wuffs_base__private_implementation__render_inf(wuffs_base__slice_u8 dst,
1432 bool neg,
1433 uint32_t options) {
1434 if (neg) {
1435 if (dst.len < 4) {
1436 return 0;
1437 }
1438 wuffs_base__store_u32le__no_bounds_check(dst.ptr, 0x666E492D); // '-Inf'le.
1439 return 4;
1440 }
1441
1442 if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1443 if (dst.len < 4) {
1444 return 0;
1445 }
1446 wuffs_base__store_u32le__no_bounds_check(dst.ptr, 0x666E492B); // '+Inf'le.
1447 return 4;
1448 }
1449
1450 if (dst.len < 3) {
1451 return 0;
1452 }
1453 wuffs_base__store_u24le__no_bounds_check(dst.ptr, 0x666E49); // 'Inf'le.
1454 return 3;
1455}
1456
1457static inline size_t //
1458wuffs_base__private_implementation__render_nan(wuffs_base__slice_u8 dst) {
1459 if (dst.len < 3) {
1460 return 0;
1461 }
1462 wuffs_base__store_u24le__no_bounds_check(dst.ptr, 0x4E614E); // 'NaN'le.
1463 return 3;
1464}
1465
1466static size_t //
1467wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1468 wuffs_base__slice_u8 dst,
1469 wuffs_base__private_implementation__high_prec_dec* h,
1470 uint32_t precision,
1471 uint32_t options) {
1472 size_t n = (h->negative ||
1473 (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN))
1474 ? 1
1475 : 0;
1476 if (h->decimal_point <= 0) {
1477 n += 1;
1478 } else {
1479 n += (size_t)(h->decimal_point);
1480 }
1481 if (precision > 0) {
1482 n += precision + 1; // +1 for the '.'.
1483 }
1484
1485 // Don't modify dst if the formatted number won't fit.
1486 if (n > dst.len) {
1487 return 0;
1488 }
1489
1490 // Align-left or align-right.
1491 uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT)
1492 ? &dst.ptr[dst.len - n]
1493 : &dst.ptr[0];
1494
1495 // Leading "±".
1496 if (h->negative) {
1497 *ptr++ = '-';
1498 } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1499 *ptr++ = '+';
1500 }
1501
1502 // Integral digits.
1503 if (h->decimal_point <= 0) {
1504 *ptr++ = '0';
1505 } else {
1506 uint32_t m =
1507 wuffs_base__u32__min(h->num_digits, (uint32_t)(h->decimal_point));
1508 uint32_t i = 0;
1509 for (; i < m; i++) {
1510 *ptr++ = (uint8_t)('0' | h->digits[i]);
1511 }
1512 for (; i < (uint32_t)(h->decimal_point); i++) {
1513 *ptr++ = '0';
1514 }
1515 }
1516
1517 // Separator and then fractional digits.
1518 if (precision > 0) {
1519 *ptr++ =
1520 (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1521 ? ','
1522 : '.';
1523 uint32_t i = 0;
1524 for (; i < precision; i++) {
1525 uint32_t j = ((uint32_t)(h->decimal_point)) + i;
1526 *ptr++ = (uint8_t)('0' | ((j < h->num_digits) ? h->digits[j] : 0));
1527 }
1528 }
1529
1530 return n;
1531}
1532
1533static size_t //
1534wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1535 wuffs_base__slice_u8 dst,
1536 wuffs_base__private_implementation__high_prec_dec* h,
1537 uint32_t precision,
1538 uint32_t options) {
1539 int32_t exp = 0;
1540 if (h->num_digits > 0) {
1541 exp = h->decimal_point - 1;
1542 }
1543 bool negative_exp = exp < 0;
1544 if (negative_exp) {
1545 exp = -exp;
1546 }
1547
1548 size_t n = (h->negative ||
1549 (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN))
1550 ? 4
1551 : 3; // Mininum 3 bytes: first digit and then "e±".
1552 if (precision > 0) {
1553 n += precision + 1; // +1 for the '.'.
1554 }
1555 n += (exp < 100) ? 2 : 3;
1556
1557 // Don't modify dst if the formatted number won't fit.
1558 if (n > dst.len) {
1559 return 0;
1560 }
1561
1562 // Align-left or align-right.
1563 uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT)
1564 ? &dst.ptr[dst.len - n]
1565 : &dst.ptr[0];
1566
1567 // Leading "±".
1568 if (h->negative) {
1569 *ptr++ = '-';
1570 } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) {
1571 *ptr++ = '+';
1572 }
1573
1574 // Integral digit.
1575 if (h->num_digits > 0) {
1576 *ptr++ = (uint8_t)('0' | h->digits[0]);
1577 } else {
1578 *ptr++ = '0';
1579 }
1580
1581 // Separator and then fractional digits.
1582 if (precision > 0) {
1583 *ptr++ =
1584 (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
1585 ? ','
1586 : '.';
1587 uint32_t i = 1;
1588 uint32_t j = wuffs_base__u32__min(h->num_digits, precision + 1);
1589 for (; i < j; i++) {
1590 *ptr++ = (uint8_t)('0' | h->digits[i]);
1591 }
1592 for (; i <= precision; i++) {
1593 *ptr++ = '0';
1594 }
1595 }
1596
1597 // Exponent: "e±" and then 2 or 3 digits.
1598 *ptr++ = 'e';
1599 *ptr++ = negative_exp ? '-' : '+';
1600 if (exp < 10) {
1601 *ptr++ = '0';
1602 *ptr++ = (uint8_t)('0' | exp);
1603 } else if (exp < 100) {
1604 *ptr++ = (uint8_t)('0' | (exp / 10));
1605 *ptr++ = (uint8_t)('0' | (exp % 10));
1606 } else {
1607 int32_t e = exp / 100;
1608 exp -= e * 100;
1609 *ptr++ = (uint8_t)('0' | e);
1610 *ptr++ = (uint8_t)('0' | (exp / 10));
1611 *ptr++ = (uint8_t)('0' | (exp % 10));
1612 }
1613
1614 return n;
1615}
1616
1617WUFFS_BASE__MAYBE_STATIC size_t //
1618wuffs_base__render_number_f64(wuffs_base__slice_u8 dst,
1619 double x,
1620 uint32_t precision,
1621 uint32_t options) {
1622 // Decompose x (64 bits) into negativity (1 bit), base-2 exponent (11 bits
1623 // with a -1023 bias) and mantissa (52 bits).
1624 uint64_t bits = wuffs_base__ieee_754_bit_representation__from_f64(x);
1625 bool neg = (bits >> 63) != 0;
1626 int32_t exp2 = ((int32_t)(bits >> 52)) & 0x7FF;
1627 uint64_t man = bits & 0x000FFFFFFFFFFFFFul;
1628
1629 // Apply the exponent bias and set the implicit top bit of the mantissa,
1630 // unless x is subnormal. Also take care of Inf and NaN.
1631 if (exp2 == 0x7FF) {
1632 if (man != 0) {
1633 return wuffs_base__private_implementation__render_nan(dst);
1634 }
1635 return wuffs_base__private_implementation__render_inf(dst, neg, options);
1636 } else if (exp2 == 0) {
1637 exp2 = -1022;
1638 } else {
1639 exp2 -= 1023;
1640 man |= 0x0010000000000000ul;
1641 }
1642
1643 // Ensure that precision isn't too large.
1644 if (precision > 4095) {
1645 precision = 4095;
1646 }
1647
1648 // Convert from the (neg, exp2, man) tuple to an HPD.
1649 wuffs_base__private_implementation__high_prec_dec h;
1650 wuffs_base__private_implementation__high_prec_dec__assign(&h, man, neg);
1651 if (h.num_digits > 0) {
1652 wuffs_base__private_implementation__high_prec_dec__lshift(
1653 &h, exp2 - 52); // 52 mantissa bits.
1654 }
1655
1656 // Handle the "%e" and "%f" formats.
1657 switch (options & (WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT |
1658 WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT)) {
1659 case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT: // The "%"f" format.
1660 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1661 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1662 &h, exp2, man);
1663 int32_t p = ((int32_t)(h.num_digits)) - h.decimal_point;
1664 precision = ((uint32_t)(wuffs_base__i32__max(0, p)));
1665 } else {
1666 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1667 &h, ((int32_t)precision) + h.decimal_point);
1668 }
1669 return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1670 dst, &h, precision, options);
1671
1672 case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT: // The "%e" format.
1673 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1674 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1675 &h, exp2, man);
1676 precision = (h.num_digits > 0) ? (h.num_digits - 1) : 0;
1677 } else {
1678 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1679 &h, ((int32_t)precision) + 1);
1680 }
1681 return wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1682 dst, &h, precision, options);
1683 }
1684
1685 // We have the "%g" format and so precision means the number of significant
1686 // digits, not the number of digits after the decimal separator. Perform
1687 // rounding and determine whether to use "%e" or "%f".
1688 int32_t e_threshold = 0;
1689 if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) {
1690 wuffs_base__private_implementation__high_prec_dec__round_just_enough(
1691 &h, exp2, man);
1692 precision = h.num_digits;
1693 e_threshold = 6;
1694 } else {
1695 if (precision == 0) {
1696 precision = 1;
1697 }
1698 wuffs_base__private_implementation__high_prec_dec__round_nearest(
1699 &h, ((int32_t)precision));
1700 e_threshold = ((int32_t)precision);
1701 int32_t nd = ((int32_t)(h.num_digits));
1702 if ((e_threshold > nd) && (nd >= h.decimal_point)) {
1703 e_threshold = nd;
1704 }
1705 }
1706
1707 // Use the "%e" format if the exponent is large.
1708 int32_t e = h.decimal_point - 1;
1709 if ((e < -4) || (e_threshold <= e)) {
1710 uint32_t p = wuffs_base__u32__min(precision, h.num_digits);
1711 return wuffs_base__private_implementation__high_prec_dec__render_exponent_present(
1712 dst, &h, (p > 0) ? (p - 1) : 0, options);
1713 }
1714
1715 // Use the "%f" format otherwise.
1716 int32_t p = ((int32_t)precision);
1717 if (p > h.decimal_point) {
1718 p = ((int32_t)(h.num_digits));
1719 }
1720 precision = ((uint32_t)(wuffs_base__i32__max(0, p - h.decimal_point)));
1721 return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent(
1722 dst, &h, precision, options);
1723}