1/*
2 * Tiny float64 printing and parsing library
3 *
4 * Copyright (c) 2024 Fabrice Bellard
5 *
6 * Permission is hereby granted, free of charge, to any person obtaining a copy
7 * of this software and associated documentation files (the "Software"), to deal
8 * in the Software without restriction, including without limitation the rights
9 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10 * copies of the Software, and to permit persons to whom the Software is
11 * furnished to do so, subject to the following conditions:
12 *
13 * The above copyright notice and this permission notice shall be included in
14 * all copies or substantial portions of the Software.
15 *
16 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
19 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
22 * THE SOFTWARE.
23 */
24#include <stdlib.h>
25#include <stdio.h>
26#include <stdarg.h>
27#include <inttypes.h>
28#include <string.h>
29#include <assert.h>
30#include <ctype.h>
31#include <sys/time.h>
32#include <math.h>
33#include <setjmp.h>
34
35#include "cutils.h"
36#include "dtoa.h"
37
38/*
39 TODO:
40 - test n_digits=101 instead of 100
41 - simplify subnormal handling
42 - reduce max memory usage
43 - free format: could add shortcut if exact result
44 - use 64 bit limb_t when possible
45 - use another algorithm for free format dtoa in base 10 (ryu ?)
46*/
47
48#define USE_POW5_TABLE
49/* use fast path to print small integers in free format */
50#define USE_FAST_INT
51
52#define LIMB_LOG2_BITS 5
53
54#define LIMB_BITS (1 << LIMB_LOG2_BITS)
55
56typedef int32_t slimb_t;
57typedef uint32_t limb_t;
58typedef uint64_t dlimb_t;
59
60#define LIMB_DIGITS 9
61
62#define JS_RADIX_MAX 36
63
64#define DBIGNUM_LEN_MAX 52 /* ~ 2^(1072+53)*36^100 (dtoa) */
65#define MANT_LEN_MAX 18 /* < 36^100 */
66
67typedef intptr_t mp_size_t;
68
69/* the represented number is sum(i, tab[i]*2^(LIMB_BITS * i)) */
70typedef struct {
71 int len; /* >= 1 */
72 limb_t tab[];
73} mpb_t;
74
75static limb_t mp_add_ui(limb_t *tab, limb_t b, size_t n)
76{
77 size_t i;
78 limb_t k, a;
79
80 k=b;
81 for(i=0;i<n;i++) {
82 if (k == 0)
83 break;
84 a = tab[i] + k;
85 k = (a < k);
86 tab[i] = a;
87 }
88 return k;
89}
90
91/* tabr[] = taba[] * b + l. Return the high carry */
92static limb_t mp_mul1(limb_t *tabr, const limb_t *taba, limb_t n,
93 limb_t b, limb_t l)
94{
95 limb_t i;
96 dlimb_t t;
97
98 for(i = 0; i < n; i++) {
99 t = (dlimb_t)taba[i] * (dlimb_t)b + l;
100 tabr[i] = t;
101 l = t >> LIMB_BITS;
102 }
103 return l;
104}
105
106/* WARNING: d must be >= 2^(LIMB_BITS-1) */
107static inline limb_t udiv1norm_init(limb_t d)
108{
109 limb_t a0, a1;
110 a1 = -d - 1;
111 a0 = -1;
112 return (((dlimb_t)a1 << LIMB_BITS) | a0) / d;
113}
114
115/* return the quotient and the remainder in '*pr'of 'a1*2^LIMB_BITS+a0
116 / d' with 0 <= a1 < d. */
117static inline limb_t udiv1norm(limb_t *pr, limb_t a1, limb_t a0,
118 limb_t d, limb_t d_inv)
119{
120 limb_t n1m, n_adj, q, r, ah;
121 dlimb_t a;
122 n1m = ((slimb_t)a0 >> (LIMB_BITS - 1));
123 n_adj = a0 + (n1m & d);
124 a = (dlimb_t)d_inv * (a1 - n1m) + n_adj;
125 q = (a >> LIMB_BITS) + a1;
126 /* compute a - q * r and update q so that the remainder is between
127 0 and d - 1 */
128 a = ((dlimb_t)a1 << LIMB_BITS) | a0;
129 a = a - (dlimb_t)q * d - d;
130 ah = a >> LIMB_BITS;
131 q += 1 + ah;
132 r = (limb_t)a + (ah & d);
133 *pr = r;
134 return q;
135}
136
137static limb_t mp_div1(limb_t *tabr, const limb_t *taba, limb_t n,
138 limb_t b, limb_t r)
139{
140 slimb_t i;
141 dlimb_t a1;
142 for(i = n - 1; i >= 0; i--) {
143 a1 = ((dlimb_t)r << LIMB_BITS) | taba[i];
144 tabr[i] = a1 / b;
145 r = a1 % b;
146 }
147 return r;
148}
149
150/* r = (a + high*B^n) >> shift. Return the remainder r (0 <= r < 2^shift).
151 1 <= shift <= LIMB_BITS - 1 */
152static limb_t mp_shr(limb_t *tab_r, const limb_t *tab, mp_size_t n,
153 int shift, limb_t high)
154{
155 mp_size_t i;
156 limb_t l, a;
157
158 assert(shift >= 1 && shift < LIMB_BITS);
159 l = high;
160 for(i = n - 1; i >= 0; i--) {
161 a = tab[i];
162 tab_r[i] = (a >> shift) | (l << (LIMB_BITS - shift));
163 l = a;
164 }
165 return l & (((limb_t)1 << shift) - 1);
166}
167
168/* r = (a << shift) + low. 1 <= shift <= LIMB_BITS - 1, 0 <= low <
169 2^shift. */
170static limb_t mp_shl(limb_t *tab_r, const limb_t *tab, mp_size_t n,
171 int shift, limb_t low)
172{
173 mp_size_t i;
174 limb_t l, a;
175
176 assert(shift >= 1 && shift < LIMB_BITS);
177 l = low;
178 for(i = 0; i < n; i++) {
179 a = tab[i];
180 tab_r[i] = (a << shift) | l;
181 l = (a >> (LIMB_BITS - shift));
182 }
183 return l;
184}
185
186static no_inline limb_t mp_div1norm(limb_t *tabr, const limb_t *taba, limb_t n,
187 limb_t b, limb_t r, limb_t b_inv, int shift)
188{
189 slimb_t i;
190
191 if (shift != 0) {
192 r = (r << shift) | mp_shl(tabr, taba, n, shift, 0);
193 }
194 for(i = n - 1; i >= 0; i--) {
195 tabr[i] = udiv1norm(&r, r, taba[i], b, b_inv);
196 }
197 r >>= shift;
198 return r;
199}
200
201static __maybe_unused void mpb_dump(const char *str, const mpb_t *a)
202{
203 int i;
204
205 printf("%s= 0x", str);
206 for(i = a->len - 1; i >= 0; i--) {
207 printf("%08x", a->tab[i]);
208 if (i != 0)
209 printf("_");
210 }
211 printf("\n");
212}
213
214static void mpb_renorm(mpb_t *r)
215{
216 while (r->len > 1 && r->tab[r->len - 1] == 0)
217 r->len--;
218}
219
220#ifdef USE_POW5_TABLE
221static const uint32_t pow5_table[17] = {
222 0x00000005, 0x00000019, 0x0000007d, 0x00000271,
223 0x00000c35, 0x00003d09, 0x0001312d, 0x0005f5e1,
224 0x001dcd65, 0x009502f9, 0x02e90edd, 0x0e8d4a51,
225 0x48c27395, 0x6bcc41e9, 0x1afd498d, 0x86f26fc1,
226 0xa2bc2ec5,
227};
228
229static const uint8_t pow5h_table[4] = {
230 0x00000001, 0x00000007, 0x00000023, 0x000000b1,
231};
232
233static const uint32_t pow5_inv_table[13] = {
234 0x99999999, 0x47ae147a, 0x0624dd2f, 0xa36e2eb1,
235 0x4f8b588e, 0x0c6f7a0b, 0xad7f29ab, 0x5798ee23,
236 0x12e0be82, 0xb7cdfd9d, 0x5fd7fe17, 0x19799812,
237 0xc25c2684,
238};
239#endif
240
241/* return a^b */
242static uint64_t pow_ui(uint32_t a, uint32_t b)
243{
244 int i, n_bits;
245 uint64_t r;
246 if (b == 0)
247 return 1;
248 if (b == 1)
249 return a;
250#ifdef USE_POW5_TABLE
251 if ((a == 5 || a == 10) && b <= 17) {
252 r = pow5_table[b - 1];
253 if (b >= 14) {
254 r |= (uint64_t)pow5h_table[b - 14] << 32;
255 }
256 if (a == 10)
257 r <<= b;
258 return r;
259 }
260#endif
261 r = a;
262 n_bits = 32 - clz32(b);
263 for(i = n_bits - 2; i >= 0; i--) {
264 r *= r;
265 if ((b >> i) & 1)
266 r *= a;
267 }
268 return r;
269}
270
271static uint32_t pow_ui_inv(uint32_t *pr_inv, int *pshift, uint32_t a, uint32_t b)
272{
273 uint32_t r_inv, r;
274 int shift;
275#ifdef USE_POW5_TABLE
276 if (a == 5 && b >= 1 && b <= 13) {
277 r = pow5_table[b - 1];
278 shift = clz32(r);
279 r <<= shift;
280 r_inv = pow5_inv_table[b - 1];
281 } else
282#endif
283 {
284 r = pow_ui(a, b);
285 shift = clz32(r);
286 r <<= shift;
287 r_inv = udiv1norm_init(r);
288 }
289 *pshift = shift;
290 *pr_inv = r_inv;
291 return r;
292}
293
294enum {
295 JS_RNDN, /* round to nearest, ties to even */
296 JS_RNDNA, /* round to nearest, ties away from zero */
297 JS_RNDZ,
298};
299
300static int mpb_get_bit(const mpb_t *r, int k)
301{
302 int l;
303
304 l = (unsigned)k / LIMB_BITS;
305 k = k & (LIMB_BITS - 1);
306 if (l >= r->len)
307 return 0;
308 else
309 return (r->tab[l] >> k) & 1;
310}
311
312/* compute round(r / 2^shift). 'shift' can be negative */
313static void mpb_shr_round(mpb_t *r, int shift, int rnd_mode)
314{
315 int l, i;
316
317 if (shift == 0)
318 return;
319 if (shift < 0) {
320 shift = -shift;
321 l = (unsigned)shift / LIMB_BITS;
322 shift = shift & (LIMB_BITS - 1);
323 if (shift != 0) {
324 r->tab[r->len] = mp_shl(r->tab, r->tab, r->len, shift, 0);
325 r->len++;
326 mpb_renorm(r);
327 }
328 if (l > 0) {
329 for(i = r->len - 1; i >= 0; i--)
330 r->tab[i + l] = r->tab[i];
331 for(i = 0; i < l; i++)
332 r->tab[i] = 0;
333 r->len += l;
334 }
335 } else {
336 limb_t bit1, bit2;
337 int k, add_one;
338
339 switch(rnd_mode) {
340 default:
341 case JS_RNDZ:
342 add_one = 0;
343 break;
344 case JS_RNDN:
345 case JS_RNDNA:
346 bit1 = mpb_get_bit(r, shift - 1);
347 if (bit1) {
348 if (rnd_mode == JS_RNDNA) {
349 bit2 = 1;
350 } else {
351 /* bit2 = oring of all the bits after bit1 */
352 bit2 = 0;
353 if (shift >= 2) {
354 k = shift - 1;
355 l = (unsigned)k / LIMB_BITS;
356 k = k & (LIMB_BITS - 1);
357 for(i = 0; i < min_int(l, r->len); i++)
358 bit2 |= r->tab[i];
359 if (l < r->len)
360 bit2 |= r->tab[l] & (((limb_t)1 << k) - 1);
361 }
362 }
363 if (bit2) {
364 add_one = 1;
365 } else {
366 /* round to even */
367 add_one = mpb_get_bit(r, shift);
368 }
369 } else {
370 add_one = 0;
371 }
372 break;
373 }
374
375 l = (unsigned)shift / LIMB_BITS;
376 shift = shift & (LIMB_BITS - 1);
377 if (l >= r->len) {
378 r->len = 1;
379 r->tab[0] = add_one;
380 } else {
381 if (l > 0) {
382 r->len -= l;
383 for(i = 0; i < r->len; i++)
384 r->tab[i] = r->tab[i + l];
385 }
386 if (shift != 0) {
387 mp_shr(r->tab, r->tab, r->len, shift, 0);
388 mpb_renorm(r);
389 }
390 if (add_one) {
391 limb_t a;
392 a = mp_add_ui(r->tab, 1, r->len);
393 if (a)
394 r->tab[r->len++] = a;
395 }
396 }
397 }
398}
399
400/* return -1, 0 or 1 */
401static int mpb_cmp(const mpb_t *a, const mpb_t *b)
402{
403 mp_size_t i;
404 if (a->len < b->len)
405 return -1;
406 else if (a->len > b->len)
407 return 1;
408 for(i = a->len - 1; i >= 0; i--) {
409 if (a->tab[i] != b->tab[i]) {
410 if (a->tab[i] < b->tab[i])
411 return -1;
412 else
413 return 1;
414 }
415 }
416 return 0;
417}
418
419static void mpb_set_u64(mpb_t *r, uint64_t m)
420{
421#if LIMB_BITS == 64
422 r->tab[0] = m;
423 r->len = 1;
424#else
425 r->tab[0] = m;
426 r->tab[1] = m >> LIMB_BITS;
427 if (r->tab[1] == 0)
428 r->len = 1;
429 else
430 r->len = 2;
431#endif
432}
433
434static uint64_t mpb_get_u64(mpb_t *r)
435{
436#if LIMB_BITS == 64
437 return r->tab[0];
438#else
439 if (r->len == 1) {
440 return r->tab[0];
441 } else {
442 return r->tab[0] | ((uint64_t)r->tab[1] << LIMB_BITS);
443 }
444#endif
445}
446
447/* floor_log2() = position of the first non zero bit or -1 if zero. */
448static int mpb_floor_log2(mpb_t *a)
449{
450 limb_t v;
451 v = a->tab[a->len - 1];
452 if (v == 0)
453 return -1;
454 else
455 return a->len * LIMB_BITS - 1 - clz32(v);
456}
457
458#define MUL_LOG2_RADIX_BASE_LOG2 24
459
460/* round((1 << MUL_LOG2_RADIX_BASE_LOG2)/log2(i + 2)) */
461static const uint32_t mul_log2_radix_table[JS_RADIX_MAX - 1] = {
462 0x000000, 0xa1849d, 0x000000, 0x6e40d2,
463 0x6308c9, 0x5b3065, 0x000000, 0x50c24e,
464 0x4d104d, 0x4a0027, 0x4768ce, 0x452e54,
465 0x433d00, 0x418677, 0x000000, 0x3ea16b,
466 0x3d645a, 0x3c43c2, 0x3b3b9a, 0x3a4899,
467 0x39680b, 0x3897b3, 0x37d5af, 0x372069,
468 0x367686, 0x35d6df, 0x354072, 0x34b261,
469 0x342bea, 0x33ac62, 0x000000, 0x32bfd9,
470 0x3251dd, 0x31e8d6, 0x318465,
471};
472
473/* return floor(a / log2(radix)) for -2048 <= a <= 2047 */
474static int mul_log2_radix(int a, int radix)
475{
476 int radix_bits, mult;
477
478 if ((radix & (radix - 1)) == 0) {
479 /* if the radix is a power of two better to do it exactly */
480 radix_bits = 31 - clz32(radix);
481 if (a < 0)
482 a -= radix_bits - 1;
483 return a / radix_bits;
484 } else {
485 mult = mul_log2_radix_table[radix - 2];
486 return ((int64_t)a * mult) >> MUL_LOG2_RADIX_BASE_LOG2;
487 }
488}
489
490#if 0
491static void build_mul_log2_radix_table(void)
492{
493 int base, radix, mult, col, base_log2;
494
495 base_log2 = 24;
496 base = 1 << base_log2;
497 col = 0;
498 for(radix = 2; radix <= 36; radix++) {
499 if ((radix & (radix - 1)) == 0)
500 mult = 0;
501 else
502 mult = lrint((double)base / log2(radix));
503 printf("0x%06x, ", mult);
504 if (++col == 4) {
505 printf("\n");
506 col = 0;
507 }
508 }
509 printf("\n");
510}
511
512static void mul_log2_radix_test(void)
513{
514 int radix, i, ref, r;
515
516 for(radix = 2; radix <= 36; radix++) {
517 for(i = -2048; i <= 2047; i++) {
518 ref = (int)floor((double)i / log2(radix));
519 r = mul_log2_radix(i, radix);
520 if (ref != r) {
521 printf("ERROR: radix=%d i=%d r=%d ref=%d\n",
522 radix, i, r, ref);
523 exit(1);
524 }
525 }
526 }
527 if (0)
528 build_mul_log2_radix_table();
529}
530#endif
531
532static void u32toa_len(char *buf, uint32_t n, size_t len)
533{
534 int digit, i;
535 for(i = len - 1; i >= 0; i--) {
536 digit = n % 10;
537 n = n / 10;
538 buf[i] = digit + '0';
539 }
540}
541
542/* for power of 2 radixes. len >= 1 */
543static void u64toa_bin_len(char *buf, uint64_t n, unsigned int radix_bits, int len)
544{
545 int digit, i;
546 unsigned int mask;
547
548 mask = (1 << radix_bits) - 1;
549 for(i = len - 1; i >= 0; i--) {
550 digit = n & mask;
551 n >>= radix_bits;
552 if (digit < 10)
553 digit += '0';
554 else
555 digit += 'a' - 10;
556 buf[i] = digit;
557 }
558}
559
560/* len >= 1. 2 <= radix <= 36 */
561static void limb_to_a(char *buf, limb_t n, unsigned int radix, int len)
562{
563 int digit, i;
564
565 if (radix == 10) {
566 /* specific case with constant divisor */
567#if LIMB_BITS == 32
568 u32toa_len(buf, n, len);
569#else
570 /* XXX: optimize */
571 for(i = len - 1; i >= 0; i--) {
572 digit = (limb_t)n % 10;
573 n = (limb_t)n / 10;
574 buf[i] = digit + '0';
575 }
576#endif
577 } else {
578 for(i = len - 1; i >= 0; i--) {
579 digit = (limb_t)n % radix;
580 n = (limb_t)n / radix;
581 if (digit < 10)
582 digit += '0';
583 else
584 digit += 'a' - 10;
585 buf[i] = digit;
586 }
587 }
588}
589
590size_t u32toa(char *buf, uint32_t n)
591{
592 char buf1[10], *q;
593 size_t len;
594
595 q = buf1 + sizeof(buf1);
596 do {
597 *--q = n % 10 + '0';
598 n /= 10;
599 } while (n != 0);
600 len = buf1 + sizeof(buf1) - q;
601 memcpy(buf, q, len);
602 return len;
603}
604
605size_t i32toa(char *buf, int32_t n)
606{
607 if (n >= 0) {
608 return u32toa(buf, n);
609 } else {
610 buf[0] = '-';
611 return u32toa(buf + 1, -(uint32_t)n) + 1;
612 }
613}
614
615#ifdef USE_FAST_INT
616size_t u64toa(char *buf, uint64_t n)
617{
618 if (n < 0x100000000) {
619 return u32toa(buf, n);
620 } else {
621 uint64_t n1;
622 char *q = buf;
623 uint32_t n2;
624
625 n1 = n / 1000000000;
626 n %= 1000000000;
627 if (n1 >= 0x100000000) {
628 n2 = n1 / 1000000000;
629 n1 = n1 % 1000000000;
630 /* at most two digits */
631 if (n2 >= 10) {
632 *q++ = n2 / 10 + '0';
633 n2 %= 10;
634 }
635 *q++ = n2 + '0';
636 u32toa_len(q, n1, 9);
637 q += 9;
638 } else {
639 q += u32toa(q, n1);
640 }
641 u32toa_len(q, n, 9);
642 q += 9;
643 return q - buf;
644 }
645}
646
647size_t i64toa(char *buf, int64_t n)
648{
649 if (n >= 0) {
650 return u64toa(buf, n);
651 } else {
652 buf[0] = '-';
653 return u64toa(buf + 1, -(uint64_t)n) + 1;
654 }
655}
656
657/* XXX: only tested for 1 <= n < 2^53 */
658size_t u64toa_radix(char *buf, uint64_t n, unsigned int radix)
659{
660 int radix_bits, l;
661 if (likely(radix == 10))
662 return u64toa(buf, n);
663 if ((radix & (radix - 1)) == 0) {
664 radix_bits = 31 - clz32(radix);
665 if (n == 0)
666 l = 1;
667 else
668 l = (64 - clz64(n) + radix_bits - 1) / radix_bits;
669 u64toa_bin_len(buf, n, radix_bits, l);
670 return l;
671 } else {
672 char buf1[41], *q; /* maximum length for radix = 3 */
673 size_t len;
674 int digit;
675 q = buf1 + sizeof(buf1);
676 do {
677 digit = n % radix;
678 n /= radix;
679 if (digit < 10)
680 digit += '0';
681 else
682 digit += 'a' - 10;
683 *--q = digit;
684 } while (n != 0);
685 len = buf1 + sizeof(buf1) - q;
686 memcpy(buf, q, len);
687 return len;
688 }
689}
690
691size_t i64toa_radix(char *buf, int64_t n, unsigned int radix)
692{
693 if (n >= 0) {
694 return u64toa_radix(buf, n, radix);
695 } else {
696 buf[0] = '-';
697 return u64toa_radix(buf + 1, -(uint64_t)n, radix) + 1;
698 }
699}
700#endif /* USE_FAST_INT */
701
702static const uint8_t digits_per_limb_table[JS_RADIX_MAX - 1] = {
703#if LIMB_BITS == 32
70432,20,16,13,12,11,10,10, 9, 9, 8, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
705#else
70664,40,32,27,24,22,21,20,19,18,17,17,16,16,16,15,15,15,14,14,14,14,13,13,13,13,13,13,13,12,12,12,12,12,12,
707#endif
708};
709
710static const uint32_t radix_base_table[JS_RADIX_MAX - 1] = {
711 0x00000000, 0xcfd41b91, 0x00000000, 0x48c27395,
712 0x81bf1000, 0x75db9c97, 0x40000000, 0xcfd41b91,
713 0x3b9aca00, 0x8c8b6d2b, 0x19a10000, 0x309f1021,
714 0x57f6c100, 0x98c29b81, 0x00000000, 0x18754571,
715 0x247dbc80, 0x3547667b, 0x4c4b4000, 0x6b5a6e1d,
716 0x94ace180, 0xcaf18367, 0x0b640000, 0x0e8d4a51,
717 0x1269ae40, 0x17179149, 0x1cb91000, 0x23744899,
718 0x2b73a840, 0x34e63b41, 0x40000000, 0x4cfa3cc1,
719 0x5c13d840, 0x6d91b519, 0x81bf1000,
720};
721
722/* XXX: remove the table ? */
723static uint8_t dtoa_max_digits_table[JS_RADIX_MAX - 1] = {
724 54, 35, 28, 24, 22, 20, 19, 18, 17, 17, 16, 16, 15, 15, 15, 14, 14, 14, 14, 14, 13, 13, 13, 13, 13, 13, 13, 12, 12, 12, 12, 12, 12, 12, 12,
725};
726
727/* we limit the maximum number of significant digits for atod to about
728 128 bits of precision for non power of two bases. The only
729 requirement for Javascript is at least 20 digits in base 10. For
730 power of two bases, we do an exact rounding in all the cases. */
731static uint8_t atod_max_digits_table[JS_RADIX_MAX - 1] = {
732 64, 80, 32, 55, 49, 45, 21, 40, 38, 37, 35, 34, 33, 32, 16, 31, 30, 30, 29, 29, 28, 28, 27, 27, 27, 26, 26, 26, 26, 25, 12, 25, 25, 24, 24,
733};
734
735/* if abs(d) >= B^max_exponent, it is an overflow */
736static const int16_t max_exponent[JS_RADIX_MAX - 1] = {
737 1024, 647, 512, 442, 397, 365, 342, 324,
738 309, 297, 286, 277, 269, 263, 256, 251,
739 246, 242, 237, 234, 230, 227, 224, 221,
740 218, 216, 214, 211, 209, 207, 205, 203,
741 202, 200, 199,
742};
743
744/* if abs(d) <= B^min_exponent, it is an underflow */
745static const int16_t min_exponent[JS_RADIX_MAX - 1] = {
746-1075, -679, -538, -463, -416, -383, -359, -340,
747 -324, -311, -300, -291, -283, -276, -269, -263,
748 -258, -254, -249, -245, -242, -238, -235, -232,
749 -229, -227, -224, -222, -220, -217, -215, -214,
750 -212, -210, -208,
751};
752
753#if 0
754void build_tables(void)
755{
756 int r, j, radix, n, col, i;
757
758 /* radix_base_table */
759 for(radix = 2; radix <= 36; radix++) {
760 r = 1;
761 for(j = 0; j < digits_per_limb_table[radix - 2]; j++) {
762 r *= radix;
763 }
764 printf(" 0x%08x,", r);
765 if ((radix % 4) == 1)
766 printf("\n");
767 }
768 printf("\n");
769
770 /* dtoa_max_digits_table */
771 for(radix = 2; radix <= 36; radix++) {
772 /* Note: over estimated when the radix is a power of two */
773 printf(" %d,", 1 + (int)ceil(53.0 / log2(radix)));
774 }
775 printf("\n");
776
777 /* atod_max_digits_table */
778 for(radix = 2; radix <= 36; radix++) {
779 if ((radix & (radix - 1)) == 0) {
780 /* 64 bits is more than enough */
781 n = (int)floor(64.0 / log2(radix));
782 } else {
783 n = (int)floor(128.0 / log2(radix));
784 }
785 printf(" %d,", n);
786 }
787 printf("\n");
788
789 printf("static const int16_t max_exponent[JS_RADIX_MAX - 1] = {\n");
790 col = 0;
791 for(radix = 2; radix <= 36; radix++) {
792 printf("%5d, ", (int)ceil(1024 / log2(radix)));
793 if (++col == 8) {
794 col = 0;
795 printf("\n");
796 }
797 }
798 printf("\n};\n\n");
799
800 printf("static const int16_t min_exponent[JS_RADIX_MAX - 1] = {\n");
801 col = 0;
802 for(radix = 2; radix <= 36; radix++) {
803 printf("%5d, ", (int)floor(-1075 / log2(radix)));
804 if (++col == 8) {
805 col = 0;
806 printf("\n");
807 }
808 }
809 printf("\n};\n\n");
810
811 printf("static const uint32_t pow5_table[16] = {\n");
812 col = 0;
813 for(i = 2; i <= 17; i++) {
814 r = 1;
815 for(j = 0; j < i; j++) {
816 r *= 5;
817 }
818 printf("0x%08x, ", r);
819 if (++col == 4) {
820 col = 0;
821 printf("\n");
822 }
823 }
824 printf("\n};\n\n");
825
826 /* high part */
827 printf("static const uint8_t pow5h_table[4] = {\n");
828 col = 0;
829 for(i = 14; i <= 17; i++) {
830 uint64_t r1;
831 r1 = 1;
832 for(j = 0; j < i; j++) {
833 r1 *= 5;
834 }
835 printf("0x%08x, ", (uint32_t)(r1 >> 32));
836 if (++col == 4) {
837 col = 0;
838 printf("\n");
839 }
840 }
841 printf("\n};\n\n");
842}
843#endif
844
845/* n_digits >= 1. 0 <= dot_pos <= n_digits. If dot_pos == n_digits,
846 the dot is not displayed. 'a' is modified. */
847static int output_digits(char *buf,
848 mpb_t *a, int radix, int n_digits1,
849 int dot_pos)
850{
851 int n_digits, digits_per_limb, radix_bits, n, len;
852
853 n_digits = n_digits1;
854 if ((radix & (radix - 1)) == 0) {
855 /* radix = 2^radix_bits */
856 radix_bits = 31 - clz32(radix);
857 } else {
858 radix_bits = 0;
859 }
860 digits_per_limb = digits_per_limb_table[radix - 2];
861 if (radix_bits != 0) {
862 for(;;) {
863 n = min_int(n_digits, digits_per_limb);
864 n_digits -= n;
865 u64toa_bin_len(buf + n_digits, a->tab[0], radix_bits, n);
866 if (n_digits == 0)
867 break;
868 mpb_shr_round(a, digits_per_limb * radix_bits, JS_RNDZ);
869 }
870 } else {
871 limb_t r;
872 while (n_digits != 0) {
873 n = min_int(n_digits, digits_per_limb);
874 n_digits -= n;
875 r = mp_div1(a->tab, a->tab, a->len, radix_base_table[radix - 2], 0);
876 mpb_renorm(a);
877 limb_to_a(buf + n_digits, r, radix, n);
878 }
879 }
880
881 /* add the dot */
882 len = n_digits1;
883 if (dot_pos != n_digits1) {
884 memmove(buf + dot_pos + 1, buf + dot_pos, n_digits1 - dot_pos);
885 buf[dot_pos] = '.';
886 len++;
887 }
888 return len;
889}
890
891/* return (a, e_offset) such that a = a * (radix1*2^radix_shift)^f *
892 2^-e_offset. 'f' can be negative. */
893static int mul_pow(mpb_t *a, int radix1, int radix_shift, int f, BOOL is_int, int e)
894{
895 int e_offset, d, n, n0;
896
897 e_offset = -f * radix_shift;
898 if (radix1 != 1) {
899 d = digits_per_limb_table[radix1 - 2];
900 if (f >= 0) {
901 limb_t h, b;
902
903 b = 0;
904 n0 = 0;
905 while (f != 0) {
906 n = min_int(f, d);
907 if (n != n0) {
908 b = pow_ui(radix1, n);
909 n0 = n;
910 }
911 h = mp_mul1(a->tab, a->tab, a->len, b, 0);
912 if (h != 0) {
913 a->tab[a->len++] = h;
914 }
915 f -= n;
916 }
917 } else {
918 int extra_bits, l, shift;
919 limb_t r, rem, b, b_inv;
920
921 f = -f;
922 l = (f + d - 1) / d; /* high bound for the number of limbs (XXX: make it better) */
923 e_offset += l * LIMB_BITS;
924 if (!is_int) {
925 /* at least 'e' bits are needed in the final result for rounding */
926 extra_bits = max_int(e - mpb_floor_log2(a), 0);
927 } else {
928 /* at least two extra bits are needed in the final result
929 for rounding */
930 extra_bits = max_int(2 + e - e_offset, 0);
931 }
932 e_offset += extra_bits;
933 mpb_shr_round(a, -(l * LIMB_BITS + extra_bits), JS_RNDZ);
934
935 b = 0;
936 b_inv = 0;
937 shift = 0;
938 n0 = 0;
939 rem = 0;
940 while (f != 0) {
941 n = min_int(f, d);
942 if (n != n0) {
943 b = pow_ui_inv(&b_inv, &shift, radix1, n);
944 n0 = n;
945 }
946 r = mp_div1norm(a->tab, a->tab, a->len, b, 0, b_inv, shift);
947 rem |= r;
948 mpb_renorm(a);
949 f -= n;
950 }
951 /* if the remainder is non zero, use it for rounding */
952 a->tab[0] |= (rem != 0);
953 }
954 }
955 return e_offset;
956}
957
958/* tmp1 = round(m*2^e*radix^f). 'tmp0' is a temporary storage */
959static void mul_pow_round(mpb_t *tmp1, uint64_t m, int e, int radix1, int radix_shift, int f,
960 int rnd_mode)
961{
962 int e_offset;
963
964 mpb_set_u64(tmp1, m);
965 e_offset = mul_pow(tmp1, radix1, radix_shift, f, TRUE, e);
966 mpb_shr_round(tmp1, -e + e_offset, rnd_mode);
967}
968
969/* return round(a*2^e_offset) rounded as a float64. 'a' is modified */
970static uint64_t round_to_d(int *pe, mpb_t *a, int e_offset, int rnd_mode)
971{
972 int e;
973 uint64_t m;
974
975 if (a->tab[0] == 0 && a->len == 1) {
976 /* zero result */
977 m = 0;
978 e = 0; /* don't care */
979 } else {
980 int prec, prec1, e_min;
981 e = mpb_floor_log2(a) + 1 - e_offset;
982 prec1 = 53;
983 e_min = -1021;
984 if (e < e_min) {
985 /* subnormal result or zero */
986 prec = prec1 - (e_min - e);
987 } else {
988 prec = prec1;
989 }
990 mpb_shr_round(a, e + e_offset - prec, rnd_mode);
991 m = mpb_get_u64(a);
992 m <<= (53 - prec);
993 /* mantissa overflow due to rounding */
994 if (m >= (uint64_t)1 << 53) {
995 m >>= 1;
996 e++;
997 }
998 }
999 *pe = e;
1000 return m;
1001}
1002
1003/* return (m, e) such that m*2^(e-53) = round(a * radix^f) with 2^52
1004 <= m < 2^53 or m = 0.
1005 'a' is modified. */
1006static uint64_t mul_pow_round_to_d(int *pe, mpb_t *a,
1007 int radix1, int radix_shift, int f, int rnd_mode)
1008{
1009 int e_offset;
1010
1011 e_offset = mul_pow(a, radix1, radix_shift, f, FALSE, 55);
1012 return round_to_d(pe, a, e_offset, rnd_mode);
1013}
1014
1015#ifdef JS_DTOA_DUMP_STATS
1016static int out_len_count[17];
1017
1018void js_dtoa_dump_stats(void)
1019{
1020 int i, sum;
1021 sum = 0;
1022 for(i = 0; i < 17; i++)
1023 sum += out_len_count[i];
1024 for(i = 0; i < 17; i++) {
1025 printf("%2d %8d %5.2f%%\n",
1026 i + 1, out_len_count[i], (double)out_len_count[i] / sum * 100);
1027 }
1028}
1029#endif
1030
1031/* return a maximum bound of the string length. The bound depends on
1032 'd' only if format = JS_DTOA_FORMAT_FRAC or if JS_DTOA_EXP_DISABLED
1033 is enabled. */
1034int js_dtoa_max_len(double d, int radix, int n_digits, int flags)
1035{
1036 int fmt = flags & JS_DTOA_FORMAT_MASK;
1037 int n, e;
1038 uint64_t a;
1039
1040 if (fmt != JS_DTOA_FORMAT_FRAC) {
1041 if (fmt == JS_DTOA_FORMAT_FREE) {
1042 n = dtoa_max_digits_table[radix - 2];
1043 } else {
1044 n = n_digits;
1045 }
1046 if ((flags & JS_DTOA_EXP_MASK) == JS_DTOA_EXP_DISABLED) {
1047 /* no exponential */
1048 a = float64_as_uint64(d);
1049 e = (a >> 52) & 0x7ff;
1050 if (e == 0x7ff) {
1051 /* NaN, Infinity */
1052 n = 0;
1053 } else {
1054 e -= 1023;
1055 /* XXX: adjust */
1056 n += 10 + abs(mul_log2_radix(e - 1, radix));
1057 }
1058 } else {
1059 /* extra: sign, 1 dot and exponent "e-1000" */
1060 n += 1 + 1 + 6;
1061 }
1062 } else {
1063 a = float64_as_uint64(d);
1064 e = (a >> 52) & 0x7ff;
1065 if (e == 0x7ff) {
1066 /* NaN, Infinity */
1067 n = 0;
1068 } else {
1069 /* high bound for the integer part */
1070 e -= 1023;
1071 /* x < 2^(e + 1) */
1072 if (e < 0) {
1073 n = 1;
1074 } else {
1075 n = 2 + mul_log2_radix(e - 1, radix);
1076 }
1077 /* sign, extra digit, 1 dot */
1078 n += 1 + 1 + 1 + n_digits;
1079 }
1080 }
1081 return max_int(n, 9); /* also include NaN and [-]Infinity */
1082}
1083
1084#if defined(__SANITIZE_ADDRESS__) && 0
1085static void *dtoa_malloc(uint64_t **pptr, size_t size)
1086{
1087 return malloc(size);
1088}
1089static void dtoa_free(void *ptr)
1090{
1091 free(ptr);
1092}
1093#else
1094static void *dtoa_malloc(uint64_t **pptr, size_t size)
1095{
1096 void *ret;
1097 ret = *pptr;
1098 *pptr += (size + 7) / 8;
1099 return ret;
1100}
1101
1102static void dtoa_free(void *ptr)
1103{
1104}
1105#endif
1106
1107/* return the length */
1108int js_dtoa(char *buf, double d, int radix, int n_digits, int flags,
1109 JSDTOATempMem *tmp_mem)
1110{
1111 uint64_t a, m, *mptr = tmp_mem->mem;
1112 int e, sgn, l, E, P, i, E_max, radix1, radix_shift;
1113 char *q;
1114 mpb_t *tmp1, *mant_max;
1115 int fmt = flags & JS_DTOA_FORMAT_MASK;
1116
1117 tmp1 = dtoa_malloc(&mptr, sizeof(mpb_t) + sizeof(limb_t) * DBIGNUM_LEN_MAX);
1118 mant_max = dtoa_malloc(&mptr, sizeof(mpb_t) + sizeof(limb_t) * MANT_LEN_MAX);
1119 assert((mptr - tmp_mem->mem) <= sizeof(JSDTOATempMem) / sizeof(mptr[0]));
1120
1121 radix_shift = ctz32(radix);
1122 radix1 = radix >> radix_shift;
1123 a = float64_as_uint64(d);
1124 sgn = a >> 63;
1125 e = (a >> 52) & 0x7ff;
1126 m = a & (((uint64_t)1 << 52) - 1);
1127 q = buf;
1128 if (e == 0x7ff) {
1129 if (m == 0) {
1130 if (sgn)
1131 *q++ = '-';
1132 memcpy(q, "Infinity", 8);
1133 q += 8;
1134 } else {
1135 memcpy(q, "NaN", 3);
1136 q += 3;
1137 }
1138 goto done;
1139 } else if (e == 0) {
1140 if (m == 0) {
1141 tmp1->len = 1;
1142 tmp1->tab[0] = 0;
1143 E = 1;
1144 if (fmt == JS_DTOA_FORMAT_FREE)
1145 P = 1;
1146 else if (fmt == JS_DTOA_FORMAT_FRAC)
1147 P = n_digits + 1;
1148 else
1149 P = n_digits;
1150 /* "-0" is displayed as "0" if JS_DTOA_MINUS_ZERO is not present */
1151 if (sgn && (flags & JS_DTOA_MINUS_ZERO))
1152 *q++ = '-';
1153 goto output;
1154 }
1155 /* denormal number: convert to a normal number */
1156 l = clz64(m) - 11;
1157 e -= l - 1;
1158 m <<= l;
1159 } else {
1160 m |= (uint64_t)1 << 52;
1161 }
1162 if (sgn)
1163 *q++ = '-';
1164 /* remove the bias */
1165 e -= 1022;
1166 /* d = 2^(e-53)*m */
1167 // printf("m=0x%016" PRIx64 " e=%d\n", m, e);
1168#ifdef USE_FAST_INT
1169 if (fmt == JS_DTOA_FORMAT_FREE &&
1170 e >= 1 && e <= 53 &&
1171 (m & (((uint64_t)1 << (53 - e)) - 1)) == 0 &&
1172 (flags & JS_DTOA_EXP_MASK) != JS_DTOA_EXP_ENABLED) {
1173 m >>= 53 - e;
1174 /* 'm' is never zero */
1175 q += u64toa_radix(q, m, radix);
1176 goto done;
1177 }
1178#endif
1179
1180 /* this choice of E implies F=round(x*B^(P-E) is such as:
1181 B^(P-1) <= F < 2.B^P. */
1182 E = 1 + mul_log2_radix(e - 1, radix);
1183
1184 if (fmt == JS_DTOA_FORMAT_FREE) {
1185 int P_max, E0, e1, E_found, P_found;
1186 uint64_t m1, mant_found, mant, mant_max1;
1187 /* P_max is guarranteed to work by construction */
1188 P_max = dtoa_max_digits_table[radix - 2];
1189 E0 = E;
1190 E_found = 0;
1191 P_found = 0;
1192 mant_found = 0;
1193 /* find the minimum number of digits by successive tries */
1194 P = P_max; /* P_max is guarateed to work */
1195 for(;;) {
1196 /* mant_max always fits on 64 bits */
1197 mant_max1 = pow_ui(radix, P);
1198 /* compute the mantissa in base B */
1199 E = E0;
1200 for(;;) {
1201 /* XXX: add inexact flag */
1202 mul_pow_round(tmp1, m, e - 53, radix1, radix_shift, P - E, JS_RNDN);
1203 mant = mpb_get_u64(tmp1);
1204 if (mant < mant_max1)
1205 break;
1206 E++; /* at most one iteration is possible */
1207 }
1208 /* remove useless trailing zero digits */
1209 while ((mant % radix) == 0) {
1210 mant /= radix;
1211 P--;
1212 }
1213 /* garanteed to work for P = P_max */
1214 if (P_found == 0)
1215 goto prec_found;
1216 /* convert back to base 2 */
1217 mpb_set_u64(tmp1, mant);
1218 m1 = mul_pow_round_to_d(&e1, tmp1, radix1, radix_shift, E - P, JS_RNDN);
1219 // printf("P=%2d: m=0x%016" PRIx64 " e=%d m1=0x%016" PRIx64 " e1=%d\n", P, m, e, m1, e1);
1220 /* Note: (m, e) is never zero here, so the exponent for m1
1221 = 0 does not matter */
1222 if (m1 == m && e1 == e) {
1223 prec_found:
1224 P_found = P;
1225 E_found = E;
1226 mant_found = mant;
1227 if (P == 1)
1228 break;
1229 P--; /* try lower exponent */
1230 } else {
1231 break;
1232 }
1233 }
1234 P = P_found;
1235 E = E_found;
1236 mpb_set_u64(tmp1, mant_found);
1237#ifdef JS_DTOA_DUMP_STATS
1238 if (radix == 10) {
1239 out_len_count[P - 1]++;
1240 }
1241#endif
1242 } else if (fmt == JS_DTOA_FORMAT_FRAC) {
1243 int len;
1244
1245 assert(n_digits >= 0 && n_digits <= JS_DTOA_MAX_DIGITS);
1246 /* P = max_int(E, 1) + n_digits; */
1247 /* frac is rounded using RNDNA */
1248 mul_pow_round(tmp1, m, e - 53, radix1, radix_shift, n_digits, JS_RNDNA);
1249
1250 /* we add one extra digit on the left and remove it if needed
1251 to avoid testing if the result is < radix^P */
1252 len = output_digits(q, tmp1, radix, max_int(E + 1, 1) + n_digits,
1253 max_int(E + 1, 1));
1254 if (q[0] == '0' && len >= 2 && q[1] != '.') {
1255 len--;
1256 memmove(q, q + 1, len);
1257 }
1258 q += len;
1259 goto done;
1260 } else {
1261 int pow_shift;
1262 assert(n_digits >= 1 && n_digits <= JS_DTOA_MAX_DIGITS);
1263 P = n_digits;
1264 /* mant_max = radix^P */
1265 mant_max->len = 1;
1266 mant_max->tab[0] = 1;
1267 pow_shift = mul_pow(mant_max, radix1, radix_shift, P, FALSE, 0);
1268 mpb_shr_round(mant_max, pow_shift, JS_RNDZ);
1269
1270 for(;;) {
1271 /* fixed and frac are rounded using RNDNA */
1272 mul_pow_round(tmp1, m, e - 53, radix1, radix_shift, P - E, JS_RNDNA);
1273 if (mpb_cmp(tmp1, mant_max) < 0)
1274 break;
1275 E++; /* at most one iteration is possible */
1276 }
1277 }
1278 output:
1279 if (fmt == JS_DTOA_FORMAT_FIXED)
1280 E_max = n_digits;
1281 else
1282 E_max = dtoa_max_digits_table[radix - 2] + 4;
1283 if ((flags & JS_DTOA_EXP_MASK) == JS_DTOA_EXP_ENABLED ||
1284 ((flags & JS_DTOA_EXP_MASK) == JS_DTOA_EXP_AUTO && (E <= -6 || E > E_max))) {
1285 q += output_digits(q, tmp1, radix, P, 1);
1286 E--;
1287 if (radix == 10) {
1288 *q++ = 'e';
1289 } else if (radix1 == 1 && radix_shift <= 4) {
1290 E *= radix_shift;
1291 *q++ = 'p';
1292 } else {
1293 *q++ = '@';
1294 }
1295 if (E < 0) {
1296 *q++ = '-';
1297 E = -E;
1298 } else {
1299 *q++ = '+';
1300 }
1301 q += u32toa(q, E);
1302 } else if (E <= 0) {
1303 *q++ = '0';
1304 *q++ = '.';
1305 for(i = 0; i < -E; i++)
1306 *q++ = '0';
1307 q += output_digits(q, tmp1, radix, P, P);
1308 } else {
1309 q += output_digits(q, tmp1, radix, P, min_int(P, E));
1310 for(i = 0; i < E - P; i++)
1311 *q++ = '0';
1312 }
1313 done:
1314 *q = '\0';
1315 dtoa_free(mant_max);
1316 dtoa_free(tmp1);
1317 return q - buf;
1318}
1319
1320static inline int to_digit(int c)
1321{
1322 if (c >= '0' && c <= '9')
1323 return c - '0';
1324 else if (c >= 'A' && c <= 'Z')
1325 return c - 'A' + 10;
1326 else if (c >= 'a' && c <= 'z')
1327 return c - 'a' + 10;
1328 else
1329 return 36;
1330}
1331
1332/* r = r * radix_base + a. radix_base = 0 means radix_base = 2^32 */
1333static void mpb_mul1_base(mpb_t *r, limb_t radix_base, limb_t a)
1334{
1335 int i;
1336 if (r->tab[0] == 0 && r->len == 1) {
1337 r->tab[0] = a;
1338 } else {
1339 if (radix_base == 0) {
1340 for(i = r->len; i >= 0; i--) {
1341 r->tab[i + 1] = r->tab[i];
1342 }
1343 r->tab[0] = a;
1344 } else {
1345 r->tab[r->len] = mp_mul1(r->tab, r->tab, r->len,
1346 radix_base, a);
1347 }
1348 r->len++;
1349 mpb_renorm(r);
1350 }
1351}
1352
1353/* XXX: add fast path for small integers */
1354double js_atod(const char *str, const char **pnext, int radix, int flags,
1355 JSATODTempMem *tmp_mem)
1356{
1357 uint64_t *mptr = tmp_mem->mem;
1358 const char *p, *p_start;
1359 limb_t cur_limb, radix_base, extra_digits;
1360 int is_neg, digit_count, limb_digit_count, digits_per_limb, sep, radix1, radix_shift;
1361 int radix_bits, expn, e, max_digits, expn_offset, dot_pos, sig_pos, pos;
1362 mpb_t *tmp0;
1363 double dval;
1364 BOOL is_bin_exp, is_zero, expn_overflow;
1365 uint64_t m, a;
1366
1367 tmp0 = dtoa_malloc(&mptr, sizeof(mpb_t) + sizeof(limb_t) * DBIGNUM_LEN_MAX);
1368 assert((mptr - tmp_mem->mem) <= sizeof(JSATODTempMem) / sizeof(mptr[0]));
1369 /* optional separator between digits */
1370 sep = (flags & JS_ATOD_ACCEPT_UNDERSCORES) ? '_' : 256;
1371
1372 p = str;
1373 is_neg = 0;
1374 if (p[0] == '+') {
1375 p++;
1376 p_start = p;
1377 } else if (p[0] == '-') {
1378 is_neg = 1;
1379 p++;
1380 p_start = p;
1381 } else {
1382 p_start = p;
1383 }
1384
1385 if (p[0] == '0') {
1386 if ((p[1] == 'x' || p[1] == 'X') &&
1387 (radix == 0 || radix == 16)) {
1388 p += 2;
1389 radix = 16;
1390 } else if ((p[1] == 'o' || p[1] == 'O') &&
1391 radix == 0 && (flags & JS_ATOD_ACCEPT_BIN_OCT)) {
1392 p += 2;
1393 radix = 8;
1394 } else if ((p[1] == 'b' || p[1] == 'B') &&
1395 radix == 0 && (flags & JS_ATOD_ACCEPT_BIN_OCT)) {
1396 p += 2;
1397 radix = 2;
1398 } else if ((p[1] >= '0' && p[1] <= '9') &&
1399 radix == 0 && (flags & JS_ATOD_ACCEPT_LEGACY_OCTAL)) {
1400 int i;
1401 sep = 256;
1402 for (i = 1; (p[i] >= '0' && p[i] <= '7'); i++)
1403 continue;
1404 if (p[i] == '8' || p[i] == '9')
1405 goto no_prefix;
1406 p += 1;
1407 radix = 8;
1408 } else {
1409 goto no_prefix;
1410 }
1411 /* there must be a digit after the prefix */
1412 if (to_digit((uint8_t)*p) >= radix)
1413 goto fail;
1414 no_prefix: ;
1415 } else {
1416 if (!(flags & JS_ATOD_INT_ONLY) && strstart(p, "Infinity", &p))
1417 goto overflow;
1418 }
1419 if (radix == 0)
1420 radix = 10;
1421
1422 cur_limb = 0;
1423 expn_offset = 0;
1424 digit_count = 0;
1425 limb_digit_count = 0;
1426 max_digits = atod_max_digits_table[radix - 2];
1427 digits_per_limb = digits_per_limb_table[radix - 2];
1428 radix_base = radix_base_table[radix - 2];
1429 radix_shift = ctz32(radix);
1430 radix1 = radix >> radix_shift;
1431 if (radix1 == 1) {
1432 /* radix = 2^radix_bits */
1433 radix_bits = radix_shift;
1434 } else {
1435 radix_bits = 0;
1436 }
1437 tmp0->len = 1;
1438 tmp0->tab[0] = 0;
1439 extra_digits = 0;
1440 pos = 0;
1441 dot_pos = -1;
1442 /* skip leading zeros */
1443 for(;;) {
1444 if (*p == '.' && (p > p_start || to_digit(p[1]) < radix) &&
1445 !(flags & JS_ATOD_INT_ONLY)) {
1446 if (*p == sep)
1447 goto fail;
1448 if (dot_pos >= 0)
1449 break;
1450 dot_pos = pos;
1451 p++;
1452 }
1453 if (*p == sep && p > p_start && p[1] == '0')
1454 p++;
1455 if (*p != '0')
1456 break;
1457 p++;
1458 pos++;
1459 }
1460
1461 sig_pos = pos;
1462 for(;;) {
1463 limb_t c;
1464 if (*p == '.' && (p > p_start || to_digit(p[1]) < radix) &&
1465 !(flags & JS_ATOD_INT_ONLY)) {
1466 if (*p == sep)
1467 goto fail;
1468 if (dot_pos >= 0)
1469 break;
1470 dot_pos = pos;
1471 p++;
1472 }
1473 if (*p == sep && p > p_start && to_digit(p[1]) < radix)
1474 p++;
1475 c = to_digit(*p);
1476 if (c >= radix)
1477 break;
1478 p++;
1479 pos++;
1480 if (digit_count < max_digits) {
1481 /* XXX: could be faster when radix_bits != 0 */
1482 cur_limb = cur_limb * radix + c;
1483 limb_digit_count++;
1484 if (limb_digit_count == digits_per_limb) {
1485 mpb_mul1_base(tmp0, radix_base, cur_limb);
1486 cur_limb = 0;
1487 limb_digit_count = 0;
1488 }
1489 digit_count++;
1490 } else {
1491 extra_digits |= c;
1492 }
1493 }
1494 if (limb_digit_count != 0) {
1495 mpb_mul1_base(tmp0, pow_ui(radix, limb_digit_count), cur_limb);
1496 }
1497 if (digit_count == 0) {
1498 is_zero = TRUE;
1499 expn_offset = 0;
1500 } else {
1501 is_zero = FALSE;
1502 if (dot_pos < 0)
1503 dot_pos = pos;
1504 expn_offset = sig_pos + digit_count - dot_pos;
1505 }
1506
1507 /* Use the extra digits for rounding if the base is a power of
1508 two. Otherwise they are just truncated. */
1509 if (radix_bits != 0 && extra_digits != 0) {
1510 tmp0->tab[0] |= 1;
1511 }
1512
1513 /* parse the exponent, if any */
1514 expn = 0;
1515 expn_overflow = FALSE;
1516 is_bin_exp = FALSE;
1517 if (!(flags & JS_ATOD_INT_ONLY) &&
1518 ((radix == 10 && (*p == 'e' || *p == 'E')) ||
1519 (radix != 10 && (*p == '@' ||
1520 (radix_bits >= 1 && radix_bits <= 4 && (*p == 'p' || *p == 'P'))))) &&
1521 p > p_start) {
1522 BOOL exp_is_neg;
1523 int c;
1524 is_bin_exp = (*p == 'p' || *p == 'P');
1525 p++;
1526 exp_is_neg = 0;
1527 if (*p == '+') {
1528 p++;
1529 } else if (*p == '-') {
1530 exp_is_neg = 1;
1531 p++;
1532 }
1533 c = to_digit(*p);
1534 if (c >= 10)
1535 goto fail; /* XXX: could stop before the exponent part */
1536 expn = c;
1537 p++;
1538 for(;;) {
1539 if (*p == sep && to_digit(p[1]) < 10)
1540 p++;
1541 c = to_digit(*p);
1542 if (c >= 10)
1543 break;
1544 if (!expn_overflow) {
1545 if (unlikely(expn > ((INT32_MAX - 2 - 9) / 10))) {
1546 expn_overflow = TRUE;
1547 } else {
1548 expn = expn * 10 + c;
1549 }
1550 }
1551 p++;
1552 }
1553 if (exp_is_neg)
1554 expn = -expn;
1555 /* if zero result, the exponent can be arbitrarily large */
1556 if (!is_zero && expn_overflow) {
1557 if (exp_is_neg)
1558 a = 0;
1559 else
1560 a = (uint64_t)0x7ff << 52; /* infinity */
1561 goto done;
1562 }
1563 }
1564
1565 if (p == p_start)
1566 goto fail;
1567
1568 if (is_zero) {
1569 a = 0;
1570 } else {
1571 int expn1;
1572 if (radix_bits != 0) {
1573 if (!is_bin_exp)
1574 expn *= radix_bits;
1575 expn -= expn_offset * radix_bits;
1576 expn1 = expn + digit_count * radix_bits;
1577 if (expn1 >= 1024 + radix_bits)
1578 goto overflow;
1579 else if (expn1 <= -1075)
1580 goto underflow;
1581 m = round_to_d(&e, tmp0, -expn, JS_RNDN);
1582 } else {
1583 expn -= expn_offset;
1584 expn1 = expn + digit_count;
1585 if (expn1 >= max_exponent[radix - 2] + 1)
1586 goto overflow;
1587 else if (expn1 <= min_exponent[radix - 2])
1588 goto underflow;
1589 m = mul_pow_round_to_d(&e, tmp0, radix1, radix_shift, expn, JS_RNDN);
1590 }
1591 if (m == 0) {
1592 underflow:
1593 a = 0;
1594 } else if (e > 1024) {
1595 overflow:
1596 /* overflow */
1597 a = (uint64_t)0x7ff << 52;
1598 } else if (e < -1073) {
1599 /* underflow */
1600 /* XXX: check rounding */
1601 a = 0;
1602 } else if (e < -1021) {
1603 /* subnormal */
1604 a = m >> (-e - 1021);
1605 } else {
1606 a = ((uint64_t)(e + 1022) << 52) | (m & (((uint64_t)1 << 52) - 1));
1607 }
1608 }
1609 done:
1610 a |= (uint64_t)is_neg << 63;
1611 dval = uint64_as_float64(a);
1612 done1:
1613 if (pnext)
1614 *pnext = p;
1615 dtoa_free(tmp0);
1616 return dval;
1617 fail:
1618 dval = NAN;
1619 goto done1;
1620}
1621