dtoa.c (44875B)
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 56 typedef int32_t slimb_t; 57 typedef uint32_t limb_t; 58 typedef 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 67 typedef intptr_t mp_size_t; 68 69 /* the represented number is sum(i, tab[i]*2^(LIMB_BITS * i)) */ 70 typedef struct { 71 int len; /* >= 1 */ 72 limb_t tab[]; 73 } mpb_t; 74 75 static 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 */ 92 static 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) */ 107 static 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. */ 117 static 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 137 static 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 */ 152 static 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. */ 170 static 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 186 static 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 201 static __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 214 static 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 221 static 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 229 static const uint8_t pow5h_table[4] = { 230 0x00000001, 0x00000007, 0x00000023, 0x000000b1, 231 }; 232 233 static 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 */ 242 static 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 271 static 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 294 enum { 295 JS_RNDN, /* round to nearest, ties to even */ 296 JS_RNDNA, /* round to nearest, ties away from zero */ 297 JS_RNDZ, 298 }; 299 300 static 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 */ 313 static 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 */ 401 static 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 419 static 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 434 static 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. */ 448 static 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)) */ 461 static 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 */ 474 static 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 491 static 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 512 static 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 532 static 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 */ 543 static 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 */ 561 static 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 590 size_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 605 size_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 616 size_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 647 size_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 */ 658 size_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 691 size_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 702 static const uint8_t digits_per_limb_table[JS_RADIX_MAX - 1] = { 703 #if LIMB_BITS == 32 704 32,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 706 64,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 710 static 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 ? */ 723 static 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. */ 731 static 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 */ 736 static 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 */ 745 static 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 754 void 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. */ 847 static 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. */ 893 static 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 */ 959 static 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 */ 970 static 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. */ 1006 static 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 1016 static int out_len_count[17]; 1017 1018 void 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. */ 1034 int 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 1085 static void *dtoa_malloc(uint64_t **pptr, size_t size) 1086 { 1087 return malloc(size); 1088 } 1089 static void dtoa_free(void *ptr) 1090 { 1091 free(ptr); 1092 } 1093 #else 1094 static 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 1102 static void dtoa_free(void *ptr) 1103 { 1104 } 1105 #endif 1106 1107 /* return the length */ 1108 int 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 1320 static 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 */ 1333 static 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 */ 1354 double 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 }