/js/src/dtoa.c
C | 3255 lines | 2761 code | 139 blank | 355 comment | 531 complexity | cf5fabfc007467b921c62685fda7dccf MD5 | raw file
Large files files are truncated, but you can click here to view the full file
1/* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */ 2/**************************************************************** 3 * 4 * The author of this software is David M. Gay. 5 * 6 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. 7 * 8 * Permission to use, copy, modify, and distribute this software for any 9 * purpose without fee is hereby granted, provided that this entire notice 10 * is included in all copies of any software which is or includes a copy 11 * or modification of this software and in all copies of the supporting 12 * documentation for such software. 13 * 14 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 15 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY 16 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 17 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 18 * 19 ***************************************************************/ 20 21/* Please send bug reports to David M. Gay (dmg at acm dot org, 22 * with " at " changed at "@" and " dot " changed to "."). */ 23 24/* On a machine with IEEE extended-precision registers, it is 25 * necessary to specify double-precision (53-bit) rounding precision 26 * before invoking strtod or dtoa. If the machine uses (the equivalent 27 * of) Intel 80x87 arithmetic, the call 28 * _control87(PC_53, MCW_PC); 29 * does this with many compilers. Whether this or another call is 30 * appropriate depends on the compiler; for this to work, it may be 31 * necessary to #include "float.h" or another system-dependent header 32 * file. 33 */ 34 35/* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 36 * 37 * This strtod returns a nearest machine number to the input decimal 38 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 39 * broken by the IEEE round-even rule. Otherwise ties are broken by 40 * biased rounding (add half and chop). 41 * 42 * Inspired loosely by William D. Clinger's paper "How to Read Floating 43 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 44 * 45 * Modifications: 46 * 47 * 1. We only require IEEE, IBM, or VAX double-precision 48 * arithmetic (not IEEE double-extended). 49 * 2. We get by with floating-point arithmetic in a case that 50 * Clinger missed -- when we're computing d * 10^n 51 * for a small integer d and the integer n is not too 52 * much larger than 22 (the maximum integer k for which 53 * we can represent 10^k exactly), we may be able to 54 * compute (d*10^k) * 10^(e-k) with just one roundoff. 55 * 3. Rather than a bit-at-a-time adjustment of the binary 56 * result in the hard case, we use floating-point 57 * arithmetic to determine the adjustment to within 58 * one bit; only in really hard cases do we need to 59 * compute a second residual. 60 * 4. Because of 3., we don't need a large table of powers of 10 61 * for ten-to-e (just some small tables, e.g. of 10^k 62 * for 0 <= k <= 22). 63 */ 64 65/* 66 * #define IEEE_8087 for IEEE-arithmetic machines where the least 67 * significant byte has the lowest address. 68 * #define IEEE_MC68k for IEEE-arithmetic machines where the most 69 * significant byte has the lowest address. 70 * #define Long int on machines with 32-bit ints and 64-bit longs. 71 * #define IBM for IBM mainframe-style floating-point arithmetic. 72 * #define VAX for VAX-style floating-point arithmetic (D_floating). 73 * #define No_leftright to omit left-right logic in fast floating-point 74 * computation of dtoa. 75 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 76 * and strtod and dtoa should round accordingly. 77 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 78 * and Honor_FLT_ROUNDS is not #defined. 79 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 80 * that use extended-precision instructions to compute rounded 81 * products and quotients) with IBM. 82 * #define ROUND_BIASED for IEEE-format with biased rounding. 83 * #define Inaccurate_Divide for IEEE-format with correctly rounded 84 * products but inaccurate quotients, e.g., for Intel i860. 85 * #define NO_LONG_LONG on machines that do not have a "long long" 86 * integer type (of >= 64 bits). On such machines, you can 87 * #define Just_16 to store 16 bits per 32-bit Long when doing 88 * high-precision integer arithmetic. Whether this speeds things 89 * up or slows things down depends on the machine and the number 90 * being converted. If long long is available and the name is 91 * something other than "long long", #define Llong to be the name, 92 * and if "unsigned Llong" does not work as an unsigned version of 93 * Llong, #define #ULLong to be the corresponding unsigned type. 94 * #define KR_headers for old-style C function headers. 95 * #define Bad_float_h if your system lacks a float.h or if it does not 96 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, 97 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. 98 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) 99 * if memory is available and otherwise does something you deem 100 * appropriate. If MALLOC is undefined, malloc will be invoked 101 * directly -- and assumed always to succeed. Similarly, if you 102 * want something other than the system's free() to be called to 103 * recycle memory acquired from MALLOC, #define FREE to be the 104 * name of the alternate routine. (Unless you #define 105 * NO_GLOBAL_STATE and call destroydtoa, FREE or free is only 106 * called in pathological cases, e.g., in a dtoa call after a dtoa 107 * return in mode 3 with thousands of digits requested.) 108 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making 109 * memory allocations from a private pool of memory when possible. 110 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, 111 * unless #defined to be a different length. This default length 112 * suffices to get rid of MALLOC calls except for unusual cases, 113 * such as decimal-to-binary conversion of a very long string of 114 * digits. The longest string dtoa can return is about 751 bytes 115 * long. For conversions by strtod of strings of 800 digits and 116 * all dtoa conversions in single-threaded executions with 8-byte 117 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte 118 * pointers, PRIVATE_MEM >= 7112 appears adequate. 119 * #define MULTIPLE_THREADS if the system offers preemptively scheduled 120 * multiple threads. In this case, you must provide (or suitably 121 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed 122 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed 123 * in pow5mult, ensures lazy evaluation of only one copy of high 124 * powers of 5; omitting this lock would introduce a small 125 * probability of wasting memory, but would otherwise be harmless.) 126 * You must also invoke freedtoa(s) to free the value s returned by 127 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. 128 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that 129 * avoids underflows on inputs whose result does not underflow. 130 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format 131 * floating-point numbers and flushes underflows to zero rather 132 * than implementing gradual underflow, then you must also #define 133 * Sudden_Underflow. 134 * #define USE_LOCALE to use the current locale's decimal_point value. 135 * #define SET_INEXACT if IEEE arithmetic is being used and extra 136 * computation should be done to set the inexact flag when the 137 * result is inexact and avoid setting inexact when the result 138 * is exact. In this case, dtoa.c must be compiled in 139 * an environment, perhaps provided by #include "dtoa.c" in a 140 * suitable wrapper, that defines two functions, 141 * int get_inexact(void); 142 * void clear_inexact(void); 143 * such that get_inexact() returns a nonzero value if the 144 * inexact bit is already set, and clear_inexact() sets the 145 * inexact bit to 0. When SET_INEXACT is #defined, strtod 146 * also does extra computations to set the underflow and overflow 147 * flags when appropriate (i.e., when the result is tiny and 148 * inexact or when it is a numeric value rounded to +-infinity). 149 * #define NO_ERRNO if strtod should not assign errno = ERANGE when 150 * the result overflows to +-Infinity or underflows to 0. 151 * #define NO_GLOBAL_STATE to avoid defining any non-const global or 152 * static variables. Instead the necessary state is stored in an 153 * opaque struct, DtoaState, a pointer to which must be passed to 154 * every entry point. Two new functions are added to the API: 155 * DtoaState *newdtoa(void); 156 * void destroydtoa(DtoaState *); 157 */ 158 159#ifndef Long 160#define Long long 161#endif 162#ifndef ULong 163typedef unsigned Long ULong; 164#endif 165 166#ifdef DEBUG 167#include "stdio.h" 168#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 169#endif 170 171#include "stdlib.h" 172#include "string.h" 173 174#ifdef USE_LOCALE 175#include "locale.h" 176#endif 177 178#ifdef MALLOC 179#ifdef KR_headers 180extern char *MALLOC(); 181#else 182extern void *MALLOC(size_t); 183#endif 184#else 185#define MALLOC malloc 186#endif 187 188#ifndef FREE 189#define FREE free 190#endif 191 192#ifndef Omit_Private_Memory 193#ifndef PRIVATE_MEM 194#define PRIVATE_MEM 2304 195#endif 196#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) 197#endif 198 199#undef IEEE_Arith 200#undef Avoid_Underflow 201#ifdef IEEE_MC68k 202#define IEEE_Arith 203#endif 204#ifdef IEEE_8087 205#define IEEE_Arith 206#endif 207 208#include "errno.h" 209 210#ifdef Bad_float_h 211 212#ifdef IEEE_Arith 213#define DBL_DIG 15 214#define DBL_MAX_10_EXP 308 215#define DBL_MAX_EXP 1024 216#define FLT_RADIX 2 217#endif /*IEEE_Arith*/ 218 219#ifdef IBM 220#define DBL_DIG 16 221#define DBL_MAX_10_EXP 75 222#define DBL_MAX_EXP 63 223#define FLT_RADIX 16 224#define DBL_MAX 7.2370055773322621e+75 225#endif 226 227#ifdef VAX 228#define DBL_DIG 16 229#define DBL_MAX_10_EXP 38 230#define DBL_MAX_EXP 127 231#define FLT_RADIX 2 232#define DBL_MAX 1.7014118346046923e+38 233#endif 234 235#ifndef LONG_MAX 236#define LONG_MAX 2147483647 237#endif 238 239#else /* ifndef Bad_float_h */ 240#include "float.h" 241#endif /* Bad_float_h */ 242 243#ifndef __MATH_H__ 244#include "math.h" 245#endif 246 247#ifdef __cplusplus 248extern "C" { 249#endif 250 251#ifndef CONST 252#ifdef KR_headers 253#define CONST /* blank */ 254#else 255#define CONST const 256#endif 257#endif 258 259#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 260Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. 261#endif 262 263typedef union { double d; ULong L[2]; } U; 264 265#define dval(x) ((x).d) 266#ifdef IEEE_8087 267#define word0(x) ((x).L[1]) 268#define word1(x) ((x).L[0]) 269#else 270#define word0(x) ((x).L[0]) 271#define word1(x) ((x).L[1]) 272#endif 273 274/* The following definition of Storeinc is appropriate for MIPS processors. 275 * An alternative that might be better on some machines is 276 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 277 */ 278#if defined(IEEE_8087) + defined(VAX) 279#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 280((unsigned short *)a)[0] = (unsigned short)c, a++) 281#else 282#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 283((unsigned short *)a)[1] = (unsigned short)c, a++) 284#endif 285 286/* #define P DBL_MANT_DIG */ 287/* Ten_pmax = floor(P*log(2)/log(5)) */ 288/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 289/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 290/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 291 292#ifdef IEEE_Arith 293#define Exp_shift 20 294#define Exp_shift1 20 295#define Exp_msk1 0x100000 296#define Exp_msk11 0x100000 297#define Exp_mask 0x7ff00000 298#define P 53 299#define Bias 1023 300#define Emin (-1022) 301#define Exp_1 0x3ff00000 302#define Exp_11 0x3ff00000 303#define Ebits 11 304#define Frac_mask 0xfffff 305#define Frac_mask1 0xfffff 306#define Ten_pmax 22 307#define Bletch 0x10 308#define Bndry_mask 0xfffff 309#define Bndry_mask1 0xfffff 310#define LSB 1 311#define Sign_bit 0x80000000 312#define Log2P 1 313#define Tiny0 0 314#define Tiny1 1 315#define Quick_max 14 316#define Int_max 14 317#ifndef NO_IEEE_Scale 318#define Avoid_Underflow 319#ifdef Flush_Denorm /* debugging option */ 320#undef Sudden_Underflow 321#endif 322#endif 323 324#ifndef Flt_Rounds 325#ifdef FLT_ROUNDS 326#define Flt_Rounds FLT_ROUNDS 327#else 328#define Flt_Rounds 1 329#endif 330#endif /*Flt_Rounds*/ 331 332#ifdef Honor_FLT_ROUNDS 333#define Rounding rounding 334#undef Check_FLT_ROUNDS 335#define Check_FLT_ROUNDS 336#else 337#define Rounding Flt_Rounds 338#endif 339 340#else /* ifndef IEEE_Arith */ 341#undef Check_FLT_ROUNDS 342#undef Honor_FLT_ROUNDS 343#undef SET_INEXACT 344#undef Sudden_Underflow 345#define Sudden_Underflow 346#ifdef IBM 347#undef Flt_Rounds 348#define Flt_Rounds 0 349#define Exp_shift 24 350#define Exp_shift1 24 351#define Exp_msk1 0x1000000 352#define Exp_msk11 0x1000000 353#define Exp_mask 0x7f000000 354#define P 14 355#define Bias 65 356#define Exp_1 0x41000000 357#define Exp_11 0x41000000 358#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 359#define Frac_mask 0xffffff 360#define Frac_mask1 0xffffff 361#define Bletch 4 362#define Ten_pmax 22 363#define Bndry_mask 0xefffff 364#define Bndry_mask1 0xffffff 365#define LSB 1 366#define Sign_bit 0x80000000 367#define Log2P 4 368#define Tiny0 0x100000 369#define Tiny1 0 370#define Quick_max 14 371#define Int_max 15 372#else /* VAX */ 373#undef Flt_Rounds 374#define Flt_Rounds 1 375#define Exp_shift 23 376#define Exp_shift1 7 377#define Exp_msk1 0x80 378#define Exp_msk11 0x800000 379#define Exp_mask 0x7f80 380#define P 56 381#define Bias 129 382#define Exp_1 0x40800000 383#define Exp_11 0x4080 384#define Ebits 8 385#define Frac_mask 0x7fffff 386#define Frac_mask1 0xffff007f 387#define Ten_pmax 24 388#define Bletch 2 389#define Bndry_mask 0xffff007f 390#define Bndry_mask1 0xffff007f 391#define LSB 0x10000 392#define Sign_bit 0x8000 393#define Log2P 1 394#define Tiny0 0x80 395#define Tiny1 0 396#define Quick_max 15 397#define Int_max 15 398#endif /* IBM, VAX */ 399#endif /* IEEE_Arith */ 400 401#ifndef IEEE_Arith 402#define ROUND_BIASED 403#endif 404 405#ifdef RND_PRODQUOT 406#define rounded_product(a,b) a = rnd_prod(a, b) 407#define rounded_quotient(a,b) a = rnd_quot(a, b) 408#ifdef KR_headers 409extern double rnd_prod(), rnd_quot(); 410#else 411extern double rnd_prod(double, double), rnd_quot(double, double); 412#endif 413#else 414#define rounded_product(a,b) a *= b 415#define rounded_quotient(a,b) a /= b 416#endif 417 418#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 419#define Big1 0xffffffff 420 421#ifndef Pack_32 422#define Pack_32 423#endif 424 425#ifdef KR_headers 426#define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff) 427#else 428#define FFFFFFFF 0xffffffffUL 429#endif 430 431#ifdef NO_LONG_LONG 432#undef ULLong 433#ifdef Just_16 434#undef Pack_32 435/* When Pack_32 is not defined, we store 16 bits per 32-bit Long. 436 * This makes some inner loops simpler and sometimes saves work 437 * during multiplications, but it often seems to make things slightly 438 * slower. Hence the default is now to store 32 bits per Long. 439 */ 440#endif 441#else /* long long available */ 442#ifndef Llong 443#define Llong long long 444#endif 445#ifndef ULLong 446#define ULLong unsigned Llong 447#endif 448#endif /* NO_LONG_LONG */ 449 450#ifndef MULTIPLE_THREADS 451#define ACQUIRE_DTOA_LOCK(n) /*nothing*/ 452#define FREE_DTOA_LOCK(n) /*nothing*/ 453#endif 454 455#define Kmax 7 456 457 struct 458Bigint { 459 struct Bigint *next; 460 int k, maxwds, sign, wds; 461 ULong x[1]; 462 }; 463 464 typedef struct Bigint Bigint; 465 466#ifdef NO_GLOBAL_STATE 467#ifdef MULTIPLE_THREADS 468#error "cannot have both NO_GLOBAL_STATE and MULTIPLE_THREADS" 469#endif 470 struct 471DtoaState { 472#define DECLARE_GLOBAL_STATE /* nothing */ 473#else 474#define DECLARE_GLOBAL_STATE static 475#endif 476 477 DECLARE_GLOBAL_STATE Bigint *freelist[Kmax+1]; 478 DECLARE_GLOBAL_STATE Bigint *p5s; 479#ifndef Omit_Private_Memory 480 DECLARE_GLOBAL_STATE double private_mem[PRIVATE_mem]; 481 DECLARE_GLOBAL_STATE double *pmem_next 482#ifndef NO_GLOBAL_STATE 483 = private_mem 484#endif 485 ; 486#endif 487#ifdef NO_GLOBAL_STATE 488 }; 489 typedef struct DtoaState DtoaState; 490#ifdef KR_headers 491#define STATE_PARAM state, 492#define STATE_PARAM_DECL DtoaState *state; 493#else 494#define STATE_PARAM DtoaState *state, 495#endif 496#define PASS_STATE state, 497#define GET_STATE(field) (state->field) 498 499 static DtoaState * 500newdtoa(void) 501{ 502 DtoaState *state = (DtoaState *) MALLOC(sizeof(DtoaState)); 503 if (state) { 504 memset(state, 0, sizeof(DtoaState)); 505#ifndef Omit_Private_Memory 506 state->pmem_next = state->private_mem; 507#endif 508 } 509 return state; 510} 511 512 static void 513destroydtoa 514#ifdef KR_headers 515 (state) STATE_PARAM_DECL 516#else 517 (DtoaState *state) 518#endif 519{ 520 int i; 521 Bigint *v, *next; 522 523 for (i = 0; i <= Kmax; i++) { 524 for (v = GET_STATE(freelist)[i]; v; v = next) { 525 next = v->next; 526#ifndef Omit_Private_Memory 527 if ((double*)v < GET_STATE(private_mem) || 528 (double*)v >= GET_STATE(private_mem) + PRIVATE_mem) 529#endif 530 FREE((void*)v); 531 } 532 } 533 FREE((void *)state); 534} 535 536#else 537#define STATE_PARAM /* nothing */ 538#define STATE_PARAM_DECL /* nothing */ 539#define PASS_STATE /* nothing */ 540#define GET_STATE(name) name 541#endif 542 543 static Bigint * 544Balloc 545#ifdef KR_headers 546 (STATE_PARAM k) STATE_PARAM_DECL int k; 547#else 548 (STATE_PARAM int k) 549#endif 550{ 551 int x; 552 Bigint *rv; 553#ifndef Omit_Private_Memory 554 size_t len; 555#endif 556 557 ACQUIRE_DTOA_LOCK(0); 558 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */ 559 /* but this case seems very unlikely. */ 560 if (k <= Kmax && (rv = GET_STATE(freelist)[k])) 561 GET_STATE(freelist)[k] = rv->next; 562 else { 563 x = 1 << k; 564#ifdef Omit_Private_Memory 565 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); 566#else 567 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 568 /sizeof(double); 569 if (k <= Kmax && GET_STATE(pmem_next) - GET_STATE(private_mem) + len <= PRIVATE_mem) { 570 rv = (Bigint*)GET_STATE(pmem_next); 571 GET_STATE(pmem_next) += len; 572 } 573 else 574 rv = (Bigint*)MALLOC(len*sizeof(double)); 575#endif 576 rv->k = k; 577 rv->maxwds = x; 578 } 579 FREE_DTOA_LOCK(0); 580 rv->sign = rv->wds = 0; 581 return rv; 582 } 583 584 static void 585Bfree 586#ifdef KR_headers 587 (STATE_PARAM v) STATE_PARAM_DECL Bigint *v; 588#else 589 (STATE_PARAM Bigint *v) 590#endif 591{ 592 if (v) { 593 if (v->k > Kmax) 594 FREE((void*)v); 595 else { 596 ACQUIRE_DTOA_LOCK(0); 597 v->next = GET_STATE(freelist)[v->k]; 598 GET_STATE(freelist)[v->k] = v; 599 FREE_DTOA_LOCK(0); 600 } 601 } 602 } 603 604#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 605y->wds*sizeof(Long) + 2*sizeof(int)) 606 607 static Bigint * 608multadd 609#ifdef KR_headers 610 (STATE_PARAM b, m, a) STATE_PARAM_DECL Bigint *b; int m, a; 611#else 612 (STATE_PARAM Bigint *b, int m, int a) /* multiply by m and add a */ 613#endif 614{ 615 int i, wds; 616#ifdef ULLong 617 ULong *x; 618 ULLong carry, y; 619#else 620 ULong carry, *x, y; 621#ifdef Pack_32 622 ULong xi, z; 623#endif 624#endif 625 Bigint *b1; 626 627 wds = b->wds; 628 x = b->x; 629 i = 0; 630 carry = a; 631 do { 632#ifdef ULLong 633 y = *x * (ULLong)m + carry; 634 carry = y >> 32; 635 *x++ = (ULong) y & FFFFFFFF; 636#else 637#ifdef Pack_32 638 xi = *x; 639 y = (xi & 0xffff) * m + carry; 640 z = (xi >> 16) * m + (y >> 16); 641 carry = z >> 16; 642 *x++ = (z << 16) + (y & 0xffff); 643#else 644 y = *x * m + carry; 645 carry = y >> 16; 646 *x++ = y & 0xffff; 647#endif 648#endif 649 } 650 while(++i < wds); 651 if (carry) { 652 if (wds >= b->maxwds) { 653 b1 = Balloc(PASS_STATE b->k+1); 654 Bcopy(b1, b); 655 Bfree(PASS_STATE b); 656 b = b1; 657 } 658 b->x[wds++] = (ULong) carry; 659 b->wds = wds; 660 } 661 return b; 662 } 663 664 static Bigint * 665s2b 666#ifdef KR_headers 667 (STATE_PARAM s, nd0, nd, y9) STATE_PARAM_DECL CONST char *s; int nd0, nd; ULong y9; 668#else 669 (STATE_PARAM CONST char *s, int nd0, int nd, ULong y9) 670#endif 671{ 672 Bigint *b; 673 int i, k; 674 Long x, y; 675 676 x = (nd + 8) / 9; 677 for(k = 0, y = 1; x > y; y <<= 1, k++) ; 678#ifdef Pack_32 679 b = Balloc(PASS_STATE k); 680 b->x[0] = y9; 681 b->wds = 1; 682#else 683 b = Balloc(PASS_STATE k+1); 684 b->x[0] = y9 & 0xffff; 685 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; 686#endif 687 688 i = 9; 689 if (9 < nd0) { 690 s += 9; 691 do b = multadd(PASS_STATE b, 10, *s++ - '0'); 692 while(++i < nd0); 693 s++; 694 } 695 else 696 s += 10; 697 for(; i < nd; i++) 698 b = multadd(PASS_STATE b, 10, *s++ - '0'); 699 return b; 700 } 701 702 static int 703hi0bits 704#ifdef KR_headers 705 (x) register ULong x; 706#else 707 (register ULong x) 708#endif 709{ 710 register int k = 0; 711 712 if (!(x & 0xffff0000)) { 713 k = 16; 714 x <<= 16; 715 } 716 if (!(x & 0xff000000)) { 717 k += 8; 718 x <<= 8; 719 } 720 if (!(x & 0xf0000000)) { 721 k += 4; 722 x <<= 4; 723 } 724 if (!(x & 0xc0000000)) { 725 k += 2; 726 x <<= 2; 727 } 728 if (!(x & 0x80000000)) { 729 k++; 730 if (!(x & 0x40000000)) 731 return 32; 732 } 733 return k; 734 } 735 736 static int 737lo0bits 738#ifdef KR_headers 739 (y) ULong *y; 740#else 741 (ULong *y) 742#endif 743{ 744 register int k; 745 register ULong x = *y; 746 747 if (x & 7) { 748 if (x & 1) 749 return 0; 750 if (x & 2) { 751 *y = x >> 1; 752 return 1; 753 } 754 *y = x >> 2; 755 return 2; 756 } 757 k = 0; 758 if (!(x & 0xffff)) { 759 k = 16; 760 x >>= 16; 761 } 762 if (!(x & 0xff)) { 763 k += 8; 764 x >>= 8; 765 } 766 if (!(x & 0xf)) { 767 k += 4; 768 x >>= 4; 769 } 770 if (!(x & 0x3)) { 771 k += 2; 772 x >>= 2; 773 } 774 if (!(x & 1)) { 775 k++; 776 x >>= 1; 777 if (!x) 778 return 32; 779 } 780 *y = x; 781 return k; 782 } 783 784 static Bigint * 785i2b 786#ifdef KR_headers 787 (STATE_PARAM i) STATE_PARAM_DECL int i; 788#else 789 (STATE_PARAM int i) 790#endif 791{ 792 Bigint *b; 793 794 b = Balloc(PASS_STATE 1); 795 b->x[0] = i; 796 b->wds = 1; 797 return b; 798 } 799 800 static Bigint * 801mult 802#ifdef KR_headers 803 (STATE_PARAM a, b) STATE_PARAM_DECL Bigint *a, *b; 804#else 805 (STATE_PARAM Bigint *a, Bigint *b) 806#endif 807{ 808 Bigint *c; 809 int k, wa, wb, wc; 810 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 811 ULong y; 812#ifdef ULLong 813 ULLong carry, z; 814#else 815 ULong carry, z; 816#ifdef Pack_32 817 ULong z2; 818#endif 819#endif 820 821 if (a->wds < b->wds) { 822 c = a; 823 a = b; 824 b = c; 825 } 826 k = a->k; 827 wa = a->wds; 828 wb = b->wds; 829 wc = wa + wb; 830 if (wc > a->maxwds) 831 k++; 832 c = Balloc(PASS_STATE k); 833 for(x = c->x, xa = x + wc; x < xa; x++) 834 *x = 0; 835 xa = a->x; 836 xae = xa + wa; 837 xb = b->x; 838 xbe = xb + wb; 839 xc0 = c->x; 840#ifdef ULLong 841 for(; xb < xbe; xc0++) { 842 if ((y = *xb++)) { 843 x = xa; 844 xc = xc0; 845 carry = 0; 846 do { 847 z = *x++ * (ULLong)y + *xc + carry; 848 carry = z >> 32; 849 *xc++ = (ULong) z & FFFFFFFF; 850 } 851 while(x < xae); 852 *xc = (ULong) carry; 853 } 854 } 855#else 856#ifdef Pack_32 857 for(; xb < xbe; xb++, xc0++) { 858 if (y = *xb & 0xffff) { 859 x = xa; 860 xc = xc0; 861 carry = 0; 862 do { 863 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 864 carry = z >> 16; 865 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 866 carry = z2 >> 16; 867 Storeinc(xc, z2, z); 868 } 869 while(x < xae); 870 *xc = carry; 871 } 872 if (y = *xb >> 16) { 873 x = xa; 874 xc = xc0; 875 carry = 0; 876 z2 = *xc; 877 do { 878 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 879 carry = z >> 16; 880 Storeinc(xc, z, z2); 881 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 882 carry = z2 >> 16; 883 } 884 while(x < xae); 885 *xc = z2; 886 } 887 } 888#else 889 for(; xb < xbe; xc0++) { 890 if (y = *xb++) { 891 x = xa; 892 xc = xc0; 893 carry = 0; 894 do { 895 z = *x++ * y + *xc + carry; 896 carry = z >> 16; 897 *xc++ = z & 0xffff; 898 } 899 while(x < xae); 900 *xc = carry; 901 } 902 } 903#endif 904#endif 905 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 906 c->wds = wc; 907 return c; 908 } 909 910 static Bigint * 911pow5mult 912#ifdef KR_headers 913 (STATE_PARAM b, k) STATE_PARAM_DECL Bigint *b; int k; 914#else 915 (STATE_PARAM Bigint *b, int k) 916#endif 917{ 918 Bigint *b1, *p5, *p51; 919 int i; 920 static CONST int p05[3] = { 5, 25, 125 }; 921 922 if ((i = k & 3)) 923 b = multadd(PASS_STATE b, p05[i-1], 0); 924 925 if (!(k >>= 2)) 926 return b; 927 if (!(p5 = GET_STATE(p5s))) { 928 /* first time */ 929#ifdef MULTIPLE_THREADS 930 ACQUIRE_DTOA_LOCK(1); 931 if (!(p5 = p5s)) { 932 p5 = p5s = i2b(625); 933 p5->next = 0; 934 } 935 FREE_DTOA_LOCK(1); 936#else 937 p5 = GET_STATE(p5s) = i2b(PASS_STATE 625); 938 p5->next = 0; 939#endif 940 } 941 for(;;) { 942 if (k & 1) { 943 b1 = mult(PASS_STATE b, p5); 944 Bfree(PASS_STATE b); 945 b = b1; 946 } 947 if (!(k >>= 1)) 948 break; 949 if (!(p51 = p5->next)) { 950#ifdef MULTIPLE_THREADS 951 ACQUIRE_DTOA_LOCK(1); 952 if (!(p51 = p5->next)) { 953 p51 = p5->next = mult(p5,p5); 954 p51->next = 0; 955 } 956 FREE_DTOA_LOCK(1); 957#else 958 p51 = p5->next = mult(PASS_STATE p5,p5); 959 p51->next = 0; 960#endif 961 } 962 p5 = p51; 963 } 964 return b; 965 } 966 967 static Bigint * 968lshift 969#ifdef KR_headers 970 (STATE_PARAM b, k) STATE_PARAM_DECL Bigint *b; int k; 971#else 972 (STATE_PARAM Bigint *b, int k) 973#endif 974{ 975 int i, k1, n, n1; 976 Bigint *b1; 977 ULong *x, *x1, *xe, z; 978 979#ifdef Pack_32 980 n = k >> 5; 981#else 982 n = k >> 4; 983#endif 984 k1 = b->k; 985 n1 = n + b->wds + 1; 986 for(i = b->maxwds; n1 > i; i <<= 1) 987 k1++; 988 b1 = Balloc(PASS_STATE k1); 989 x1 = b1->x; 990 for(i = 0; i < n; i++) 991 *x1++ = 0; 992 x = b->x; 993 xe = x + b->wds; 994#ifdef Pack_32 995 if (k &= 0x1f) { 996 k1 = 32 - k; 997 z = 0; 998 do { 999 *x1++ = *x << k | z; 1000 z = *x++ >> k1; 1001 } 1002 while(x < xe); 1003 if ((*x1 = z)) 1004 ++n1; 1005 } 1006#else 1007 if (k &= 0xf) { 1008 k1 = 16 - k; 1009 z = 0; 1010 do { 1011 *x1++ = *x << k & 0xffff | z; 1012 z = *x++ >> k1; 1013 } 1014 while(x < xe); 1015 if (*x1 = z) 1016 ++n1; 1017 } 1018#endif 1019 else do 1020 *x1++ = *x++; 1021 while(x < xe); 1022 b1->wds = n1 - 1; 1023 Bfree(PASS_STATE b); 1024 return b1; 1025 } 1026 1027 static int 1028cmp 1029#ifdef KR_headers 1030 (a, b) Bigint *a, *b; 1031#else 1032 (Bigint *a, Bigint *b) 1033#endif 1034{ 1035 ULong *xa, *xa0, *xb, *xb0; 1036 int i, j; 1037 1038 i = a->wds; 1039 j = b->wds; 1040#ifdef DEBUG 1041 if (i > 1 && !a->x[i-1]) 1042 Bug("cmp called with a->x[a->wds-1] == 0"); 1043 if (j > 1 && !b->x[j-1]) 1044 Bug("cmp called with b->x[b->wds-1] == 0"); 1045#endif 1046 if (i -= j) 1047 return i; 1048 xa0 = a->x; 1049 xa = xa0 + j; 1050 xb0 = b->x; 1051 xb = xb0 + j; 1052 for(;;) { 1053 if (*--xa != *--xb) 1054 return *xa < *xb ? -1 : 1; 1055 if (xa <= xa0) 1056 break; 1057 } 1058 return 0; 1059 } 1060 1061 static Bigint * 1062diff 1063#ifdef KR_headers 1064 (STATE_PARAM a, b) STATE_PARAM_DECL Bigint *a, *b; 1065#else 1066 (STATE_PARAM Bigint *a, Bigint *b) 1067#endif 1068{ 1069 Bigint *c; 1070 int i, wa, wb; 1071 ULong *xa, *xae, *xb, *xbe, *xc; 1072#ifdef ULLong 1073 ULLong borrow, y; 1074#else 1075 ULong borrow, y; 1076#ifdef Pack_32 1077 ULong z; 1078#endif 1079#endif 1080 1081 i = cmp(a,b); 1082 if (!i) { 1083 c = Balloc(PASS_STATE 0); 1084 c->wds = 1; 1085 c->x[0] = 0; 1086 return c; 1087 } 1088 if (i < 0) { 1089 c = a; 1090 a = b; 1091 b = c; 1092 i = 1; 1093 } 1094 else 1095 i = 0; 1096 c = Balloc(PASS_STATE a->k); 1097 c->sign = i; 1098 wa = a->wds; 1099 xa = a->x; 1100 xae = xa + wa; 1101 wb = b->wds; 1102 xb = b->x; 1103 xbe = xb + wb; 1104 xc = c->x; 1105 borrow = 0; 1106#ifdef ULLong 1107 do { 1108 y = (ULLong)*xa++ - *xb++ - borrow; 1109 borrow = y >> 32 & (ULong)1; 1110 *xc++ = (ULong) y & FFFFFFFF; 1111 } 1112 while(xb < xbe); 1113 while(xa < xae) { 1114 y = *xa++ - borrow; 1115 borrow = y >> 32 & (ULong)1; 1116 *xc++ = (ULong) y & FFFFFFFF; 1117 } 1118#else 1119#ifdef Pack_32 1120 do { 1121 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 1122 borrow = (y & 0x10000) >> 16; 1123 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 1124 borrow = (z & 0x10000) >> 16; 1125 Storeinc(xc, z, y); 1126 } 1127 while(xb < xbe); 1128 while(xa < xae) { 1129 y = (*xa & 0xffff) - borrow; 1130 borrow = (y & 0x10000) >> 16; 1131 z = (*xa++ >> 16) - borrow; 1132 borrow = (z & 0x10000) >> 16; 1133 Storeinc(xc, z, y); 1134 } 1135#else 1136 do { 1137 y = *xa++ - *xb++ - borrow; 1138 borrow = (y & 0x10000) >> 16; 1139 *xc++ = y & 0xffff; 1140 } 1141 while(xb < xbe); 1142 while(xa < xae) { 1143 y = *xa++ - borrow; 1144 borrow = (y & 0x10000) >> 16; 1145 *xc++ = y & 0xffff; 1146 } 1147#endif 1148#endif 1149 while(!*--xc) 1150 wa--; 1151 c->wds = wa; 1152 return c; 1153 } 1154 1155 static double 1156ulp 1157#ifdef KR_headers 1158 (x) U x; 1159#else 1160 (U x) 1161#endif 1162{ 1163 register Long L; 1164 U a; 1165 1166 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 1167#ifndef Avoid_Underflow 1168#ifndef Sudden_Underflow 1169 if (L > 0) { 1170#endif 1171#endif 1172#ifdef IBM 1173 L |= Exp_msk1 >> 4; 1174#endif 1175 word0(a) = L; 1176 word1(a) = 0; 1177#ifndef Avoid_Underflow 1178#ifndef Sudden_Underflow 1179 } 1180 else { 1181 L = -L >> Exp_shift; 1182 if (L < Exp_shift) { 1183 word0(a) = 0x80000 >> L; 1184 word1(a) = 0; 1185 } 1186 else { 1187 word0(a) = 0; 1188 L -= Exp_shift; 1189 word1(a) = L >= 31 ? 1 : 1 << 31 - L; 1190 } 1191 } 1192#endif 1193#endif 1194 return dval(a); 1195 } 1196 1197 static double 1198b2d 1199#ifdef KR_headers 1200 (a, e) Bigint *a; int *e; 1201#else 1202 (Bigint *a, int *e) 1203#endif 1204{ 1205 ULong *xa, *xa0, w, y, z; 1206 int k; 1207 U d; 1208#ifdef VAX 1209 ULong d0, d1; 1210#else 1211#define d0 word0(d) 1212#define d1 word1(d) 1213#endif 1214 1215 xa0 = a->x; 1216 xa = xa0 + a->wds; 1217 y = *--xa; 1218#ifdef DEBUG 1219 if (!y) Bug("zero y in b2d"); 1220#endif 1221 k = hi0bits(y); 1222 *e = 32 - k; 1223#ifdef Pack_32 1224 if (k < Ebits) { 1225 d0 = Exp_1 | y >> (Ebits - k); 1226 w = xa > xa0 ? *--xa : 0; 1227 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k); 1228 goto ret_d; 1229 } 1230 z = xa > xa0 ? *--xa : 0; 1231 if (k -= Ebits) { 1232 d0 = Exp_1 | y << k | z >> (32 - k); 1233 y = xa > xa0 ? *--xa : 0; 1234 d1 = z << k | y >> (32 - k); 1235 } 1236 else { 1237 d0 = Exp_1 | y; 1238 d1 = z; 1239 } 1240#else 1241 if (k < Ebits + 16) { 1242 z = xa > xa0 ? *--xa : 0; 1243 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 1244 w = xa > xa0 ? *--xa : 0; 1245 y = xa > xa0 ? *--xa : 0; 1246 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 1247 goto ret_d; 1248 } 1249 z = xa > xa0 ? *--xa : 0; 1250 w = xa > xa0 ? *--xa : 0; 1251 k -= Ebits + 16; 1252 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 1253 y = xa > xa0 ? *--xa : 0; 1254 d1 = w << k + 16 | y << k; 1255#endif 1256 ret_d: 1257#ifdef VAX 1258 word0(d) = d0 >> 16 | d0 << 16; 1259 word1(d) = d1 >> 16 | d1 << 16; 1260#else 1261#undef d0 1262#undef d1 1263#endif 1264 return dval(d); 1265 } 1266 1267 static Bigint * 1268d2b 1269#ifdef KR_headers 1270 (STATE_PARAM d, e, bits) STATE_PARAM_DECL U d; int *e, *bits; 1271#else 1272 (STATE_PARAM U d, int *e, int *bits) 1273#endif 1274{ 1275 Bigint *b; 1276 int de, k; 1277 ULong *x, y, z; 1278#ifndef Sudden_Underflow 1279 int i; 1280#endif 1281#ifdef VAX 1282 ULong d0, d1; 1283 d0 = word0(d) >> 16 | word0(d) << 16; 1284 d1 = word1(d) >> 16 | word1(d) << 16; 1285#else 1286#define d0 word0(d) 1287#define d1 word1(d) 1288#endif 1289 1290#ifdef Pack_32 1291 b = Balloc(PASS_STATE 1); 1292#else 1293 b = Balloc(PASS_STATE 2); 1294#endif 1295 x = b->x; 1296 1297 z = d0 & Frac_mask; 1298 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 1299#ifdef Sudden_Underflow 1300 de = (int)(d0 >> Exp_shift); 1301#ifndef IBM 1302 z |= Exp_msk11; 1303#endif 1304#else 1305 if ((de = (int)(d0 >> Exp_shift))) 1306 z |= Exp_msk1; 1307#endif 1308#ifdef Pack_32 1309 if ((y = d1)) { 1310 if ((k = lo0bits(&y))) { 1311 x[0] = y | z << (32 - k); 1312 z >>= k; 1313 } 1314 else 1315 x[0] = y; 1316#ifndef Sudden_Underflow 1317 i = 1318#endif 1319 b->wds = (x[1] = z) ? 2 : 1; 1320 } 1321 else { 1322 k = lo0bits(&z); 1323 x[0] = z; 1324#ifndef Sudden_Underflow 1325 i = 1326#endif 1327 b->wds = 1; 1328 k += 32; 1329 } 1330#else 1331 if (y = d1) { 1332 if (k = lo0bits(&y)) 1333 if (k >= 16) { 1334 x[0] = y | z << 32 - k & 0xffff; 1335 x[1] = z >> k - 16 & 0xffff; 1336 x[2] = z >> k; 1337 i = 2; 1338 } 1339 else { 1340 x[0] = y & 0xffff; 1341 x[1] = y >> 16 | z << 16 - k & 0xffff; 1342 x[2] = z >> k & 0xffff; 1343 x[3] = z >> k+16; 1344 i = 3; 1345 } 1346 else { 1347 x[0] = y & 0xffff; 1348 x[1] = y >> 16; 1349 x[2] = z & 0xffff; 1350 x[3] = z >> 16; 1351 i = 3; 1352 } 1353 } 1354 else { 1355#ifdef DEBUG 1356 if (!z) 1357 Bug("Zero passed to d2b"); 1358#endif 1359 k = lo0bits(&z); 1360 if (k >= 16) { 1361 x[0] = z; 1362 i = 0; 1363 } 1364 else { 1365 x[0] = z & 0xffff; 1366 x[1] = z >> 16; 1367 i = 1; 1368 } 1369 k += 32; 1370 } 1371 while(!x[i]) 1372 --i; 1373 b->wds = i + 1; 1374#endif 1375#ifndef Sudden_Underflow 1376 if (de) { 1377#endif 1378#ifdef IBM 1379 *e = (de - Bias - (P-1) << 2) + k; 1380 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 1381#else 1382 *e = de - Bias - (P-1) + k; 1383 *bits = P - k; 1384#endif 1385#ifndef Sudden_Underflow 1386 } 1387 else { 1388 *e = de - Bias - (P-1) + 1 + k; 1389#ifdef Pack_32 1390 *bits = 32*i - hi0bits(x[i-1]); 1391#else 1392 *bits = (i+2)*16 - hi0bits(x[i]); 1393#endif 1394 } 1395#endif 1396 return b; 1397 } 1398#undef d0 1399#undef d1 1400 1401 static double 1402ratio 1403#ifdef KR_headers 1404 (a, b) Bigint *a, *b; 1405#else 1406 (Bigint *a, Bigint *b) 1407#endif 1408{ 1409 U da, db; 1410 int k, ka, kb; 1411 1412 dval(da) = b2d(a, &ka); 1413 dval(db) = b2d(b, &kb); 1414#ifdef Pack_32 1415 k = ka - kb + 32*(a->wds - b->wds); 1416#else 1417 k = ka - kb + 16*(a->wds - b->wds); 1418#endif 1419#ifdef IBM 1420 if (k > 0) { 1421 word0(da) += (k >> 2)*Exp_msk1; 1422 if (k &= 3) 1423 dval(da) *= 1 << k; 1424 } 1425 else { 1426 k = -k; 1427 word0(db) += (k >> 2)*Exp_msk1; 1428 if (k &= 3) 1429 dval(db) *= 1 << k; 1430 } 1431#else 1432 if (k > 0) 1433 word0(da) += k*Exp_msk1; 1434 else { 1435 k = -k; 1436 word0(db) += k*Exp_msk1; 1437 } 1438#endif 1439 return dval(da) / dval(db); 1440 } 1441 1442 static CONST double 1443tens[] = { 1444 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1445 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1446 1e20, 1e21, 1e22 1447#ifdef VAX 1448 , 1e23, 1e24 1449#endif 1450 }; 1451 1452 static CONST double 1453#ifdef IEEE_Arith 1454bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1455static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1456#ifdef Avoid_Underflow 1457 9007199254740992.*9007199254740992.e-256 1458 /* = 2^106 * 1e-53 */ 1459#else 1460 1e-256 1461#endif 1462 }; 1463/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 1464/* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 1465#define Scale_Bit 0x10 1466#define n_bigtens 5 1467#else 1468#ifdef IBM 1469bigtens[] = { 1e16, 1e32, 1e64 }; 1470static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 1471#define n_bigtens 3 1472#else 1473bigtens[] = { 1e16, 1e32 }; 1474static CONST double tinytens[] = { 1e-16, 1e-32 }; 1475#define n_bigtens 2 1476#endif 1477#endif 1478 1479 static double 1480_strtod 1481#ifdef KR_headers 1482 (STATE_PARAM s00, se) STATE_PARAM_DECL CONST char *s00; char **se; 1483#else 1484 (STATE_PARAM CONST char *s00, char **se) 1485#endif 1486{ 1487#ifdef Avoid_Underflow 1488 int scale; 1489#endif 1490 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 1491 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 1492 CONST char *s, *s0, *s1; 1493 double aadj, adj; 1494 U aadj1, rv, rv0; 1495 Long L; 1496 ULong y, z; 1497 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 1498#ifdef SET_INEXACT 1499 int inexact, oldinexact; 1500#endif 1501#ifdef Honor_FLT_ROUNDS 1502 int rounding; 1503#endif 1504#ifdef USE_LOCALE 1505 CONST char *s2; 1506#endif 1507 1508#ifdef __GNUC__ 1509 delta = bb = bd = bs = 0; 1510#endif 1511 1512 sign = nz0 = nz = 0; 1513 dval(rv) = 0.; 1514 for(s = s00;;s++) switch(*s) { 1515 case '-': 1516 sign = 1; 1517 /* no break */ 1518 case '+': 1519 if (*++s) 1520 goto break2; 1521 /* no break */ 1522 case 0: 1523 goto ret0; 1524 case '\t': 1525 case '\n': 1526 case '\v': 1527 case '\f': 1528 case '\r': 1529 case ' ': 1530 continue; 1531 default: 1532 goto break2; 1533 } 1534 break2: 1535 if (*s == '0') { 1536 nz0 = 1; 1537 while(*++s == '0') ; 1538 if (!*s) 1539 goto ret; 1540 } 1541 s0 = s; 1542 y = z = 0; 1543 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 1544 if (nd < 9) 1545 y = 10*y + c - '0'; 1546 else if (nd < 16) 1547 z = 10*z + c - '0'; 1548 nd0 = nd; 1549#ifdef USE_LOCALE 1550 s1 = localeconv()->decimal_point; 1551 if (c == *s1) { 1552 c = '.'; 1553 if (*++s1) { 1554 s2 = s; 1555 for(;;) { 1556 if (*++s2 != *s1) { 1557 c = 0; 1558 break; 1559 } 1560 if (!*++s1) { 1561 s = s2; 1562 break; 1563 } 1564 } 1565 } 1566 } 1567#endif 1568 if (c == '.') { 1569 c = *++s; 1570 if (!nd) { 1571 for(; c == '0'; c = *++s) 1572 nz++; 1573 if (c > '0' && c <= '9') { 1574 s0 = s; 1575 nf += nz; 1576 nz = 0; 1577 goto have_dig; 1578 } 1579 goto dig_done; 1580 } 1581 for(; c >= '0' && c <= '9'; c = *++s) { 1582 have_dig: 1583 nz++; 1584 if (c -= '0') { 1585 nf += nz; 1586 for(i = 1; i < nz; i++) 1587 if (nd++ < 9) 1588 y *= 10; 1589 else if (nd <= DBL_DIG + 1) 1590 z *= 10; 1591 if (nd++ < 9) 1592 y = 10*y + c; 1593 else if (nd <= DBL_DIG + 1) 1594 z = 10*z + c; 1595 nz = 0; 1596 } 1597 } 1598 } 1599 dig_done: 1600 e = 0; 1601 if (c == 'e' || c == 'E') { 1602 if (!nd && !nz && !nz0) { 1603 goto ret0; 1604 } 1605 s00 = s; 1606 esign = 0; 1607 switch(c = *++s) { 1608 case '-': 1609 esign = 1; 1610 case '+': 1611 c = *++s; 1612 } 1613 if (c >= '0' && c <= '9') { 1614 while(c == '0') 1615 c = *++s; 1616 if (c > '0' && c <= '9') { 1617 L = c - '0'; 1618 s1 = s; 1619 while((c = *++s) >= '0' && c <= '9') 1620 L = 10*L + c - '0'; 1621 if (s - s1 > 8 || L > 19999) 1622 /* Avoid confusion from exponents 1623 * so large that e might overflow. 1624 */ 1625 e = 19999; /* safe for 16 bit ints */ 1626 else 1627 e = (int)L; 1628 if (esign) 1629 e = -e; 1630 } 1631 else 1632 e = 0; 1633 } 1634 else 1635 s = s00; 1636 } 1637 if (!nd) { 1638 if (!nz && !nz0) { 1639 ret0: 1640 s = s00; 1641 sign = 0; 1642 } 1643 goto ret; 1644 } 1645 e1 = e -= nf; 1646 1647 /* Now we have nd0 digits, starting at s0, followed by a 1648 * decimal point, followed by nd-nd0 digits. The number we're 1649 * after is the integer represented by those digits times 1650 * 10**e */ 1651 1652 if (!nd0) 1653 nd0 = nd; 1654 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 1655 dval(rv) = y; 1656 if (k > 9) { 1657#ifdef SET_INEXACT 1658 if (k > DBL_DIG) 1659 oldinexact = get_inexact(); 1660#endif 1661 dval(rv) = tens[k - 9] * dval(rv) + z; 1662 } 1663 bd0 = 0; 1664 if (nd <= DBL_DIG 1665#ifndef RND_PRODQUOT 1666#ifndef Honor_FLT_ROUNDS 1667 && Flt_Rounds == 1 1668#endif 1669#endif 1670 ) { 1671 if (!e) 1672 goto ret; 1673 if (e > 0) { 1674 if (e <= Ten_pmax) { 1675#ifdef VAX 1676 goto vax_ovfl_check; 1677#else 1678#ifdef Honor_FLT_ROUNDS 1679 /* round correctly FLT_ROUNDS = 2 or 3 */ 1680 if (sign) { 1681 rv = -rv; 1682 sign = 0; 1683 } 1684#endif 1685 /* rv = */ rounded_product(dval(rv), tens[e]); 1686 goto ret; 1687#endif 1688 } 1689 i = DBL_DIG - nd; 1690 if (e <= Ten_pmax + i) { 1691 /* A fancier test would sometimes let us do 1692 * this for larger i values. 1693 */ 1694#ifdef Honor_FLT_ROUNDS 1695 /* round correctly FLT_ROUNDS = 2 or 3 */ 1696 if (sign) { 1697 rv = -rv; 1698 sign = 0; 1699 } 1700#endif 1701 e -= i; 1702 dval(rv) *= tens[i]; 1703#ifdef VAX 1704 /* VAX exponent range is so narrow we must 1705 * worry about overflow here... 1706 */ 1707 vax_ovfl_check: 1708 word0(rv) -= P*Exp_msk1; 1709 /* rv = */ rounded_product(dval(rv), tens[e]); 1710 if ((word0(rv) & Exp_mask) 1711 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 1712 goto ovfl; 1713 word0(rv) += P*Exp_msk1; 1714#else 1715 /* rv = */ rounded_product(dval(rv), tens[e]); 1716#endif 1717 goto ret; 1718 } 1719 } 1720#ifndef Inaccurate_Divide 1721 else if (e >= -Ten_pmax) { 1722#ifdef Honor_FLT_ROUNDS 1723 /* round correctly FLT_ROUNDS = 2 or 3 */ 1724 if (sign) { 1725 rv = -rv; 1726 sign = 0; 1727 } 1728#endif 1729 /* rv = */ rounded_quotient(dval(rv), tens[-e]); 1730 goto ret; 1731 } 1732#endif 1733 } 1734 e1 += nd - k; 1735 1736#ifdef IEEE_Arith 1737#ifdef SET_INEXACT 1738 inexact = 1; 1739 if (k <= DBL_DIG) 1740 oldinexact = get_inexact(); 1741#endif 1742#ifdef Avoid_Underflow 1743 scale = 0; 1744#endif 1745#ifdef Honor_FLT_ROUNDS 1746 if ((rounding = Flt_Rounds) >= 2) { 1747 if (sign) 1748 rounding = rounding == 2 ? 0 : 2; 1749 else 1750 if (rounding != 2) 1751 rounding = 0; 1752 } 1753#endif 1754#endif /*IEEE_Arith*/ 1755 1756 /* Get starting approximation = rv * 10**e1 */ 1757 1758 if (e1 > 0) { 1759 if ((i = e1 & 15)) 1760 dval(rv) *= tens[i]; 1761 if (e1 &= ~15) { 1762 if (e1 > DBL_MAX_10_EXP) { 1763 ovfl: 1764#ifndef NO_ERRNO 1765 errno = ERANGE; 1766#endif 1767 /* Can't trust HUGE_VAL */ 1768#ifdef IEEE_Arith 1769#ifdef Honor_FLT_ROUNDS 1770 switch(rounding) { 1771 case 0: /* toward 0 */ 1772 case 3: /* toward -infinity */ 1773 word0(rv) = Big0; 1774 word1(rv) = Big1; 1775 break; 1776 default: 1777 word0(rv) = Exp_mask; 1778 word1(rv) = 0; 1779 } 1780#else /*Honor_FLT_ROUNDS*/ 1781 word0(rv) = Exp_mask; 1782 word1(rv) = 0; 1783#endif /*Honor_FLT_ROUNDS*/ 1784#ifdef SET_INEXACT 1785 /* set overflow bit */ 1786 dval(rv0) = 1e300; 1787 dval(rv0) *= dval(rv0); 1788#endif 1789#else /*IEEE_Arith*/ 1790 word0(rv) = Big0; 1791 word1(rv) = Big1; 1792#endif /*IEEE_Arith*/ 1793 if (bd0) 1794 goto retfree; 1795 goto ret; 1796 } 1797 e1 >>= 4; 1798 for(j = 0; e1 > 1; j++, e1 >>= 1) 1799 if (e1 & 1) 1800 dval(rv) *= bigtens[j]; 1801 /* The last multiplication could overflow. */ 1802 word0(rv) -= P*Exp_msk1; 1803 dval(rv) *= bigtens[j]; 1804 if ((z = word0(rv) & Exp_mask) 1805 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 1806 goto ovfl; 1807 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 1808 /* set to largest number */ 1809 /* (Can't trust DBL_MAX) */ 1810 word0(rv) = Big0; 1811 word1(rv) = Big1; 1812 } 1813 else 1814 word0(rv) += P*Exp_msk1; 1815 } 1816 } 1817 else if (e1 < 0) { 1818 e1 = -e1; 1819 if ((i = e1 & 15)) 1820 dval(rv) /= tens[i]; 1821 if (e1 >>= 4) { 1822 if (e1 >= 1 << n_bigtens) 1823 goto undfl; 1824#ifdef Avoid_Underflow 1825 if (e1 & Scale_Bit) 1826 scale = 2*P; 1827 for(j = 0; e1 > 0; j++, e1 >>= 1) 1828 if (e1 & 1) 1829 dval(rv) *= tinytens[j]; 1830 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask) 1831 >> Exp_shift)) > 0) { 1832 /* scaled rv is denormal; zap j low bits */ 1833 if (j >= 32) { 1834 word1(rv) = 0; 1835 if (j >= 53) 1836 word0(rv) = (P+2)*Exp_msk1; 1837 else 1838 word0(rv) &= 0xffffffff << (j-32); 1839 } 1840 else 1841 word1(rv) &= 0xffffffff << j; 1842 } 1843#else 1844 for(j = 0; e1 > 1; j++, e1 >>= 1) 1845 if (e1 & 1) 1846 dval(rv) *= tinytens[j]; 1847 /* The last multiplication could underflow. */ 1848 dval(rv0) = dval(rv); 1849 dval(rv) *= tinytens[j]; 1850 if (!dval(rv)) { 1851 dval(rv) = 2.*dval(rv0); 1852 dval(rv) *= tinytens[j]; 1853#endif 1854 if (!dval(rv)) { 1855 undfl: 1856 dval(rv) = 0.; 1857#ifndef NO_ERRNO 1858 errno = ERANGE; 1859#endif 1860 if (bd0) 1861 goto retfree; 1862 goto ret; 1863 } 1864#ifndef Avoid_Underflow 1865 word0(rv) = Tiny0; 1866 word1(rv) = Tiny1; 1867 /* The refinement below will clean 1868 * this approximation up. 1869 */ 1870 } 1871#endif 1872 } 1873 } 1874 1875 /* Now the hard part -- adjusting rv to the correct value.*/ 1876 1877 /* Put digits into bd: true value = bd * 10^e */ 1878 1879 bd0 = s2b(PASS_STATE s0, nd0, nd, y); 1880 1881 for(;;) { 1882 bd = Balloc(PASS_STATE bd0->k); 1883 Bcopy(bd, bd0); 1884 bb = d2b(PASS_STATE rv, &bbe, &bbbits); /* rv = bb * 2^bbe */ 1885 bs = i2b(PASS_STATE 1); 1886 1887 if (e >= 0) { 1888 bb2 = bb5 = 0; 1889 bd2 = bd5 = e; 1890 } 1891 else { 1892 bb2 = bb5 = -e; 1893 bd2 = bd5 = 0; 1894 } 1895 if (bbe >= 0) 1896 bb2 += bbe; 1897 else 1898 bd2 -= bbe; 1899 bs2 = bb2; 1900#ifdef Honor_FLT_ROUNDS 1901 if (rounding != 1) 1902 bs2++; 1903#endif 1904#ifdef Avoid_Underflow 1905 j = bbe - scale; 1906 i = j + bbbits - 1; /* logb(rv) */ 1907 if (i < Emin) /* denormal */ 1908 j += P - Emin; 1909 else 1910 j = P + 1 - bbbits; 1911#else /*Avoid_Underflow*/ 1912#ifdef Sudden_Underflow 1913#ifdef IBM 1914 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 1915#else 1916 j = P + 1 - bbbits; 1917#endif 1918#else /*Sudden_Underflow*/ 1919 j = bbe; 1920 i = j + bbbits - 1; /* logb(rv) */ 1921 if (i < Emin) /* denormal */ 1922 j += P - Emin; 1923 else 1924 j = P + 1 - bbbits; 1925#endif /*Sudden_Underflow*/ 1926#endif /*Avoid_Underflow*/ 1927 bb2 += j; 1928 bd2 += j; 1929#ifdef Avoid_Underflow 1930 bd2 += scale; 1931#endif 1932 i = bb2 < bd2 ? bb2 : bd2; 1933 if (i > bs2) 1934 i = bs2; 1935 if (i > 0) { 1936 bb2 -= i; 1937 bd2 -= i; 1938 bs2 -= i; 1939 } 1940 if (bb5 > 0) { 1941 bs = pow5mult(PASS_STATE bs, bb5); 1942 bb1 = mult(PASS_STATE bs, bb); 1943 Bfree(PASS_STATE bb); 1944 bb = bb1; 1945 } 1946 if (bb2 > 0) 1947 bb = lshift(PASS_STATE bb, bb2); 1948 if (bd5 > 0) 1949 bd = pow5mult(PASS_STATE bd, bd5); 1950 if (bd2 > 0) 1951 bd = lshift(PASS_STATE bd, bd2); 1952 if (bs2 > 0) 1953 bs = lshift(PASS_STATE bs, bs2); 1954 delta = diff(PASS_STATE bb, bd); 1955 dsign = delta->sign; 1956 delta->sign = 0; 1957 i = cmp(delta, bs); 1958#ifdef Honor_FLT_ROUNDS 1959 if (rounding != 1) { 1960 if (i < 0) { 1961 /* Error is less than an ulp */ 1962 if (!delta->x[0] && delta->wds <= 1) { 1963 /* exact */ 1964#ifdef SET_INEXACT 1965 inexact = 0; 1966#endif 1967 break; 1968 } 1969 if (rounding) { 1970 if (dsign) { 1971 adj = 1.; 1972 goto apply_adj; 1973 } 1974 } 1975 else if (!dsign) { 1976 adj = -1.; 1977 if (!word1(rv) 1978 && !(word0(rv) & Frac_mask)) { 1979 y = word0(rv) & Exp_mask; 1980#ifdef Avoid_Underflow 1981 if (!scale || y > 2*P*Exp_msk1) 1982#else 1983 if (y) 1984#endif 1985 { 1986 delta = lshift(PASS_STATE delta,Log2P); 1987 if (cmp(delta, bs) <= 0) 1988 adj = -0.5; 1989 } 1990 } 1991 apply_adj: 1992#ifdef Avoid_Underflow 1993 if (scale && (y = word0(rv) & Exp_mask) 1994 <= 2*P*Exp_msk1) 1995 word0(adj) += (2*P+1)*Exp_msk1 - y; 1996#else 1997#ifdef Sudden_Underflow 1998 if ((word0(rv) & Exp_mask) <= 1999 P*Exp_msk1) { 2000 word0(rv) += P*Exp_msk1; 2001 dval(rv) += adj*ulp(rv); 2002 word0(rv) -= P*Exp_msk1; 2003 } 2004 else 2005#endif /*Sudden_Underflow*/ 2006#endif /*Avoid_Underflow*/ 2007 dval(rv) += adj*ulp(rv); 2008 } 2009 break; 2010 } 2011 adj = ratio(delta, bs); 2012 if (adj < 1.) 2013 adj = 1.; 2014 if (adj <= 0x7ffffffe) { 2015 /* adj = rounding ? ceil(adj) : floor(adj); */ 2016 y = adj; 2017 if (y != adj) { 2018 if (!((rounding>>1) ^ dsign)) 2019 y++; 2020 adj = y; 2021 } 2022 } 2023#ifdef Avoid_Underflow 2024 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 2025 word0(adj) += (2*P+1)*Exp_msk1 - y; 2026#else 2027#ifdef Sudden_Underflow 2028 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 2029 word0(rv) += P*Exp_msk1; 2030 adj *= ulp(rv); 2031 if (dsign) 2032 dval(rv) += adj; 2033 else 2034 dval(rv) -= adj; 2035 word0(rv) -= P*Exp_msk1; 2036 goto cont; 2037 } 2038#endif /*Sudden_Underflow*/ 2039#endif /*Avoid_Underflow*/ 2040 adj *= ulp(rv); 2041 if (dsign) 2042 dval(rv) += adj; 2043 else 2044 dval(rv) -= adj; 2045 goto cont; 2046 } 2047#endif /*Honor_FLT_ROUNDS*/ 2048 2049 if (i < 0) { 2050 /* Error is less than half an ulp -- check for 2051 * special case of mantissa a power of two. 2052 */ 2053 if (dsign || word1(rv) || word0(rv) & Bndry_mask 2054#ifdef IEEE_Arith 2055#ifdef Avoid_Underflow 2056 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 2057#else 2058 || (word0(rv) & Exp_mask) <= Exp_msk1 2059#endif 2060#endif 2061 ) { 2062#ifdef SET_INEXACT 2063 if (!delta->x[0] && delta->wds <= 1) 2064 inexact = 0; 2065#endif 2066 break; 2067 } 2068 if (!delta->x[0] && delta->wds <= 1) { 2069 /* exact result */ 2070#ifdef SET_INEXACT 2071 inexact = 0; 2072#endif 2073 break; 2074 } 2075 delta = lshift(PASS_STATE delta,Log2P); 2076 if (cmp(delta, bs) > 0) 2077 goto drop_down; 2078 break; 2079 } 2080 if (i == 0) { 2081 /* exactly half-way between */ 2082 if (dsign) { 2083 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 2084 && word1(rv) == ( 2085#ifdef Avoid_Underflow 2086 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 2087 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 2088#endif 2089 0xffffffff)) { 2090 /*boundary case -- increment exponent*/ 2091 word0(rv) = (word0(rv) & Exp_mask) 2092 + Exp_msk1 2093#ifdef IBM 2094 | Exp_msk1 >> 4 2095#endif 2096 ; 2097 word1(rv) = 0; 2098#ifdef Avoid_Underflow 2099 dsign = 0; 2100#endif 2101 break; 2102 } 2103 } 2104 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 2105 drop_down: 2106 /* boundary case -- decrement exponent */ 2107#ifdef Sudden_Underflow /*{{*/ 2108 L = word0(rv) & Exp_mask; 2109#ifdef IBM 2110 if (L < Exp_msk1) 2111#else 2112#ifdef Avoid_Underflow 2113 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 2114#else 2115 if (L <= Exp_msk1) 2116#endif /*Avoid_Underflow*/ 2117#endif /*IBM*/ 2118 goto undfl; 2119 L -= Exp_msk1; 2120#else /*Sudden_Underflow}{*/ 2121#ifdef Avoid_Underflow 2122 if (scale) { 2123 L = word0(rv) & Exp_mask; 2124 if (L <= (2*P+1)*Exp_msk1) { 2125 if (L > (P+2)*Exp_msk1) 2126 /* round even ==> */ 2127 /* accept rv */ 2128 break; 2129 /* rv = smallest denormal */ 2130 goto undfl; 2131 } 2132 } 2133#endif /*Avoid_Underflow*/ 2134 L = (word0(rv) & Exp_mask) - Exp_msk1; 2135#endif /*Sudden_Underflow}}*/ 2136 word0(rv) = L | Bndry_mask1; 2137 word1(rv) = 0xffffffff; 2138#ifdef IBM 2139 goto cont; 2140#else 2141 break; 2142#endif 2143 } 2144#ifndef ROUND_BIASED 2145 if (!(word1(rv) & LSB)) 2146 break; 2147#endif 2148 if (dsign) 2149 dval(rv) += ulp(rv); 2150#ifndef ROUND_BIASED 2151 else { 2152 dval(rv) -= ulp(rv); 2153#ifndef Sudden_Underflow 2154 if (!dval(rv)) 2155 goto undfl; 2156#endif 2157 } 2158#ifdef Avoid_Underflow 2159 dsign = 1 - dsign; 2160#endif 2161#endif 2162 break; 2163 } 2164 if ((aadj = ratio(delta, bs)) <= 2.) { 2165 if (dsign) 2166 aadj = dval(aadj1) = 1.; 2167 else if (word1(rv) || word0(rv) & Bndry_mask) { 2168#ifndef Sudden_Underflow 2169 if (word1(rv) == Tiny1 && !word0(rv)) 2170 goto undfl; 2171#endif 2172 aadj = 1.; 2173 dval(aadj1) = -1.; 2174 } 2175 else { 2176 /* special case -- power of FLT_RADIX to be */ 2177 /* rounded down... */ 2178 2179 if (aadj < 2./FLT_RADIX) 2180 aadj = 1./FLT_RADIX; 2181 else 2182 aadj *= 0.5; 2183 dval(aadj1) = -aadj; 2184 } 2185 } 2186 else { 2187 aadj *= 0.5; 2188 dval(aadj1) = dsign ? aadj : -aadj; 2189#ifdef Check_FLT_ROUNDS 2190 switch(Rounding) { 2191 case 2: /* towards +infinity */ 2192 dval(aadj1) -= 0.5; 2193 break; 2194 case 0: /* towards 0 */ 2195 case 3: /* towards -infinity */ 2196 dval(aadj1) += 0.5; 2197 } 2198#else 2199 if (Flt_Rounds == 0) 2200 dval(aadj1) += 0.5; 2201#endif /*Check_FLT_ROUNDS*/ 2202 } 2203 y = word0(rv) & Exp_mask; 2204 2205 /* Check for overflow */ 2206 2207 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 2208 dval(rv0) = dval(rv); 2209 word0(rv) -= P*Exp_msk1; 2210 adj = dval(aadj1) * ulp(rv); 2211 dval(rv) += adj; 2212 if ((word0(rv) & Exp_mask) >= 2213 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 2214 if (word0(rv0) == Big0 && word1(rv0) == Big1) 2215 goto ovfl; 2216 word0(rv) = Big0; 2217 word1(rv) = Big1; 2218 goto cont; 2219 } 2220 else 2221 word0(rv) += P*Exp_msk1; 2222 } 2223 else { 2224#ifdef Avoid_Underflow 2225 if (scale && y <= 2*P*Exp_msk1) { 2226 if (aadj <= 0x7fffffff) { 2227 if ((z = (ULong) aadj) <= 0) 2228 z = 1; 2229 aadj = z; 2230 dval(aadj1) = dsign ? aadj : -aadj; 2231 } 2232 word0(aadj1) += (2*P+1)*Exp_msk1 - y; 2233 } 2234 adj = dval(aadj1) * ulp(rv); 2235 dval(rv) += adj; 2236#else 2237#ifdef Sudden_Underflow 2238 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 2239 dval(rv0) = dval(rv); 2240 word0(rv) += P*Exp_msk1; 2241 adj = dval(aadj1) * ulp(rv); 2242 dval(rv) += adj; 2243#ifdef IBM 2244 if ((word0(rv) & Exp_mask) < P*Exp_msk1) 2245#else 2246 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 2247#endif 2248 { 2249 if (word0(rv0) == Tiny0 2250 && word1(rv0) == Tiny1) 2251 goto undfl; 2252 word0(rv) = Tiny0; 2253 word1(rv) = Tiny1; 2254 goto cont; 2255 } 2256 else 2257 word0(rv) -= P*Exp_msk1; 2258 } 2259 else { 2260 adj = dval(aadj1) * ulp(rv); 2261 dval(rv) += adj; 2262 } 2263#else /*Sudden_Underflow*/ 2264 /* Compute adj so that the IEEE rounding rules will 2265 * correctly round rv + adj in some half-way cases. 2266 * If rv * ulp(rv) is denormalized (i.e., 2267 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 2268 * trouble from bits lost to denormalization; 2269 * example: 1.2e-307 . 2270 */ 2271 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 2272 dval(aadj1) = (double)(int)(aadj + 0.5); 2273 if (!dsign) 2274 dval(aadj1) = -dval(aadj1); 2275 } 2276 adj = dval(aadj1) * ulp(rv); 2277 dval(rv) += adj; 2278#endif /*Sudden_Underflow*/ 2279#endif /*Avoid_Underflow*/ 2280 } 2281 z = word0(rv) & Exp_mask; 2282#ifndef SET_INEXACT 2283#ifdef Avoid_Underflow 2284 if (!scale) 2285#endif 2286 if (y == z) { 2287 /* Can we stop now? */ 2288 L = (Long)aadj; 2289 aadj -= L; 2290 /* The tolerances below are conservative. */ 2291 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 2292 if (aadj < .4999999 || aadj > .5000001) 2293 break; 2294 } 2295 else if (aadj < .4999999/FLT_RADIX) 2296 break; 2297 } 2298#endif 2299 cont: 2300 Bfree(PASS_STATE bb); 2301 Bfree(PASS_STATE bd); 2302 Bfree(PASS_STATE bs); 2303 Bfree(PASS_STATE delta); 2304 } 2305#ifdef SET_INEXACT 2306 if (inexact) { 2307 if (!oldinexact) { 2308 word0(rv0) = Exp_1 + (70 << Exp_shift); 2309 word1(rv0) = 0; 2310 dval(rv0) += 1.; 2311 } 2312 } 2313 else if (!oldinexact) 2314 clear_inexact(); 2315#endif 2316#ifdef Avoid_Underflow 2317 if (scale) { 2318 word0(rv0) = Exp_1 - 2*P*Exp_msk1; 2319 word1(rv0) = 0; 2320 dval(rv) *= dval(rv0); 2321#ifndef NO_ERRNO 2322 /* try to avoid the bug of testing an 8087 register value */ 2323 if (word0(rv) == 0 && word1(rv) == 0) 2324 errno = ERANGE; 2325#endif 2326 } 2327#endif /* Avoid_Underflow */ 2328#ifdef SET_INEXACT 2329 if (inexact && !(word0(rv) & Exp_mask)) { 2330 /* set underflow bit */ 2331 dval(rv0) = 1e-300; 2332 dval(rv0) *= dval(rv0); 2333 } 2334#endif 2335 retfree: 2336 Bfree(PASS_STATE bb); 2337 Bfree(PASS_STATE bd); 2338 Bfree(PASS_STATE bs); 2339 Bfree(PASS_STATE bd0); 2340 Bfree(PASS_STATE delta); 2341 ret: 2342 if (se) 2343 *se = (char *)s; 2344 return sign ? -dval(rv) : dval(rv); 2345 } 2346 2347 static int 2348quorem 2349#ifdef KR_headers 2350 (b, S) Bigint *b, *S; 2351#else 2352 (Bigint *b, Bigint *S) 2353#endif 2354{ 2355 int n; 2356 ULong *bx, *bxe, q, *sx, *sxe; 2357#ifdef ULLong 2358 ULLong borrow, carry, y, ys; 2359#else 2360 ULong borrow, carry, y, ys; 2361#ifdef Pack_32 2362 ULong si, z, zs; 2363#endif 2364#endif 2365 2366 n = S->wds; 2367#ifdef DEBUG 2368 /*debug*/ if (b->wds > n) 2369 /*debug*/ Bug("oversize b in quorem"); 2370#endif 2371 if (b->wds < n) 2372 return 0; 2373 sx = S->x; 2374 sxe = sx + --n; 2375 bx = b->x; 2376 bxe = bx + n; 2377 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 2378#ifdef DEBUG 2379 /*debug*/ if (q > 9) 2380 /*debug*/ Bug("oversized quotient in quorem"); 2381#endif 2382 if (q) { 2383 borrow = 0; 2384 carry = 0; 2385 do { 2386#ifdef ULLong…
Large files files are truncated, but you can click here to view the full file