1/* 2 Note: Although strtod and dtoa seem to be present on some systems 3 they are not always included in the headers or in the libraries. 4 DCJM 6/4/00 5 6 To simplify all this strtod, dtoa and free_dtoa all have 7 a poly_ prefix. 8*/ 9/**************************************************************** 10 * 11 * The author of this software is David M. Gay. 12 * 13 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. 14 * 15 * Permission to use, copy, modify, and distribute this software for any 16 * purpose without fee is hereby granted, provided that this entire notice 17 * is included in all copies of any software which is or includes a copy 18 * or modification of this software and in all copies of the supporting 19 * documentation for such software. 20 * 21 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 22 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY 23 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 24 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 25 * 26 ***************************************************************/ 27 28/* Please send bug reports to David M. Gay (dmg at acm dot org, 29 * with " at " changed at "@" and " dot " changed to "."). */ 30 31/* On a machine with IEEE extended-precision registers, it is 32 * necessary to specify double-precision (53-bit) rounding precision 33 * before invoking strtod or dtoa. If the machine uses (the equivalent 34 * of) Intel 80x87 arithmetic, the call 35 * _control87(PC_53, MCW_PC); 36 * does this with many compilers. Whether this or another call is 37 * appropriate depends on the compiler; for this to work, it may be 38 * necessary to #include "float.h" or another system-dependent header 39 * file. 40 */ 41 42/* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 43 * (Note that IEEE arithmetic is disabled by gcc's -ffast-math flag.) 44 * 45 * This strtod returns a nearest machine number to the input decimal 46 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 47 * broken by the IEEE round-even rule. Otherwise ties are broken by 48 * biased rounding (add half and chop). 49 * 50 * Inspired loosely by William D. Clinger's paper "How to Read Floating 51 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 52 * 53 * Modifications: 54 * 55 * 1. We only require IEEE, IBM, or VAX double-precision 56 * arithmetic (not IEEE double-extended). 57 * 2. We get by with floating-point arithmetic in a case that 58 * Clinger missed -- when we're computing d * 10^n 59 * for a small integer d and the integer n is not too 60 * much larger than 22 (the maximum integer k for which 61 * we can represent 10^k exactly), we may be able to 62 * compute (d*10^k) * 10^(e-k) with just one roundoff. 63 * 3. Rather than a bit-at-a-time adjustment of the binary 64 * result in the hard case, we use floating-point 65 * arithmetic to determine the adjustment to within 66 * one bit; only in really hard cases do we need to 67 * compute a second residual. 68 * 4. Because of 3., we don't need a large table of powers of 10 69 * for ten-to-e (just some small tables, e.g. of 10^k 70 * for 0 <= k <= 22). 71 */ 72 73/* 74 * #define IEEE_8087 for IEEE-arithmetic machines where the least 75 * significant byte has the lowest address. 76 * #define IEEE_MC68k for IEEE-arithmetic machines where the most 77 * significant byte has the lowest address. 78 * #define Long int on machines with 32-bit ints and 64-bit longs. 79 * #define IBM for IBM mainframe-style floating-point arithmetic. 80 * #define VAX for VAX-style floating-point arithmetic (D_floating). 81 * #define No_leftright to omit left-right logic in fast floating-point 82 * computation of dtoa. This will cause dtoa modes 4 and 5 to be 83 * treated the same as modes 2 and 3 for some inputs. 84 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 85 * and strtod and dtoa should round accordingly. Unless Trust_FLT_ROUNDS 86 * is also #defined, fegetround() will be queried for the rounding mode. 87 * Note that both FLT_ROUNDS and fegetround() are specified by the C99 88 * standard (and are specified to be consistent, with fesetround() 89 * affecting the value of FLT_ROUNDS), but that some (Linux) systems 90 * do not work correctly in this regard, so using fegetround() is more 91 * portable than using FLT_ROUNDS directly. 92 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 93 * and Honor_FLT_ROUNDS is not #defined. 94 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 95 * that use extended-precision instructions to compute rounded 96 * products and quotients) with IBM. 97 * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic 98 * that rounds toward +Infinity. 99 * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased 100 * rounding when the underlying floating-point arithmetic uses 101 * unbiased rounding. This prevent using ordinary floating-point 102 * arithmetic when the result could be computed with one rounding error. 103 * #define Inaccurate_Divide for IEEE-format with correctly rounded 104 * products but inaccurate quotients, e.g., for Intel i860. 105 * #define NO_LONG_LONG on machines that do not have a "long long" 106 * integer type (of >= 64 bits). On such machines, you can 107 * #define Just_16 to store 16 bits per 32-bit Long when doing 108 * high-precision integer arithmetic. Whether this speeds things 109 * up or slows things down depends on the machine and the number 110 * being converted. If long long is available and the name is 111 * something other than "long long", #define Llong to be the name, 112 * and if "unsigned Llong" does not work as an unsigned version of 113 * Llong, #define #ULLong to be the corresponding unsigned type. 114 * #define KR_headers for old-style C function headers. 115 * #define Bad_float_h if your system lacks a float.h or if it does not 116 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, 117 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. 118 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) 119 * if memory is available and otherwise does something you deem 120 * appropriate. If MALLOC is undefined, malloc will be invoked 121 * directly -- and assumed always to succeed. Similarly, if you 122 * want something other than the system's free() to be called to 123 * recycle memory acquired from MALLOC, #define FREE to be the 124 * name of the alternate routine. (FREE or free is only called in 125 * pathological cases, e.g., in a dtoa call after a dtoa return in 126 * mode 3 with thousands of digits requested.) 127 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making 128 * memory allocations from a private pool of memory when possible. 129 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, 130 * unless #defined to be a different length. This default length 131 * suffices to get rid of MALLOC calls except for unusual cases, 132 * such as decimal-to-binary conversion of a very long string of 133 * digits. The longest string dtoa can return is about 751 bytes 134 * long. For conversions by strtod of strings of 800 digits and 135 * all dtoa conversions in single-threaded executions with 8-byte 136 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte 137 * pointers, PRIVATE_MEM >= 7112 appears adequate. 138 * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK 139 * #defined automatically on IEEE systems. On such systems, 140 * when INFNAN_CHECK is #defined, strtod checks 141 * for Infinity and NaN (case insensitively). On some systems 142 * (e.g., some HP systems), it may be necessary to #define NAN_WORD0 143 * appropriately -- to the most significant word of a quiet NaN. 144 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.) 145 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined, 146 * strtod also accepts (case insensitively) strings of the form 147 * NaN(x), where x is a string of hexadecimal digits and spaces; 148 * if there is only one string of hexadecimal digits, it is taken 149 * for the 52 fraction bits of the resulting NaN; if there are two 150 * or more strings of hex digits, the first is for the high 20 bits, 151 * the second and subsequent for the low 32 bits, with intervening 152 * white space ignored; but if this results in none of the 52 153 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0 154 * and NAN_WORD1 are used instead. 155 * #define MULTIPLE_THREADS if the system offers preemptively scheduled 156 * multiple threads. In this case, you must provide (or suitably 157 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed 158 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed 159 * in pow5mult, ensures lazy evaluation of only one copy of high 160 * powers of 5; omitting this lock would introduce a small 161 * probability of wasting memory, but would otherwise be harmless.) 162 * You must also invoke freedtoa(s) to free the value s returned by 163 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. 164 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that 165 * avoids underflows on inputs whose result does not underflow. 166 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format 167 * floating-point numbers and flushes underflows to zero rather 168 * than implementing gradual underflow, then you must also #define 169 * Sudden_Underflow. 170 * #define USE_LOCALE to use the current locale's decimal_point value. 171 * #define SET_INEXACT if IEEE arithmetic is being used and extra 172 * computation should be done to set the inexact flag when the 173 * result is inexact and avoid setting inexact when the result 174 * is exact. In this case, dtoa.c must be compiled in 175 * an environment, perhaps provided by #include "dtoa.c" in a 176 * suitable wrapper, that defines two functions, 177 * int get_inexact(void); 178 * void clear_inexact(void); 179 * such that get_inexact() returns a nonzero value if the 180 * inexact bit is already set, and clear_inexact() sets the 181 * inexact bit to 0. When SET_INEXACT is #defined, strtod 182 * also does extra computations to set the underflow and overflow 183 * flags when appropriate (i.e., when the result is tiny and 184 * inexact or when it is a numeric value rounded to +-infinity). 185 * #define NO_ERRNO if strtod should not assign errno = ERANGE when 186 * the result overflows to +-Infinity or underflows to 0. 187 * #define NO_HEX_FP to omit recognition of hexadecimal floating-point 188 * values by strtod. 189 * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now) 190 * to disable logic for "fast" testing of very long input strings 191 * to strtod. This testing proceeds by initially truncating the 192 * input string, then if necessary comparing the whole string with 193 * a decimal expansion to decide close cases. This logic is only 194 * used for input more than STRTOD_DIGLIM digits long (default 40). 195 */ 196 197#ifdef HAVE_CONFIG_H 198#include "config.h" 199#elif defined(_WIN32) 200#include "winconfig.h" 201#else 202#error "No configuration file" 203#endif 204 205#ifdef HAVE_SYS_PARAM_H 206#include <sys/param.h> 207#endif 208 209#include "realconv.h" 210#include "locking.h" 211 212// Try to determine the byte order to set either IEEE_8087 or IEEE_MC68k 213#ifdef __BYTE_ORDER 214#if __BYTE_ORDER == __LITTLE_ENDIAN 215#define IEEE_8087 // Little endian 216#elif __BYTE_ORDER == __BIG_ENDIAN 217#define IEEE_MC68k // Big endian 218#endif 219#endif 220 221#if !defined(IEEE_8087) && ! defined(IEEE_MC68k) 222#if defined(_WIN32) || defined(HOSTARCHITECTURE_X86) || defined (__i386__) || defined (_M_IX86) || \ 223 defined (vax) || defined (__alpha) || defined(HOSTARCHITECTURE_X86_64) || \ 224 defined(HOSTARCHITECTURE_X32) 225#define IEEE_8087 226#else 227#define IEEE_MC68k 228#endif 229#endif 230 231#if (SIZEOF_LONG == 8) 232// If "long" is the same size as "double" we need to define this. 233#define Long int 234#define ULong unsigned 235#endif 236 237#ifndef HAVE_LONG_LONG 238#define NO_LONG_LONG 239#endif 240 241#ifndef Long 242#define Long long 243#endif 244#ifndef ULong 245typedef unsigned Long ULong; 246#endif 247 248#ifdef DEBUG 249#include "stdio.h" 250#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 251#endif 252 253#include "stdlib.h" 254#include "string.h" 255 256#ifdef USE_LOCALE 257#include "locale.h" 258#endif 259 260#ifdef Honor_FLT_ROUNDS 261#ifndef Trust_FLT_ROUNDS 262#include <fenv.h> 263#endif 264#endif 265 266#ifdef MALLOC 267#ifdef KR_headers 268extern char *MALLOC(); 269#else 270extern void *MALLOC(size_t); 271#endif 272#else 273#define MALLOC malloc 274#endif 275 276#ifndef Omit_Private_Memory 277#ifndef PRIVATE_MEM 278#define PRIVATE_MEM 2304 279#endif 280#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) 281static double private_mem[PRIVATE_mem], *pmem_next = private_mem; 282#endif 283 284#undef IEEE_Arith 285#undef Avoid_Underflow 286#ifdef IEEE_MC68k 287#define IEEE_Arith 288#endif 289#ifdef IEEE_8087 290#define IEEE_Arith 291#endif 292 293#ifdef IEEE_Arith 294#ifndef NO_INFNAN_CHECK 295#undef INFNAN_CHECK 296#define INFNAN_CHECK 297#endif 298#else 299#undef INFNAN_CHECK 300#define NO_STRTOD_BIGCOMP 301#endif 302 303#include "errno.h" 304 305#ifdef Bad_float_h 306 307#ifdef IEEE_Arith 308#define DBL_DIG 15 309#define DBL_MAX_10_EXP 308 310#define DBL_MAX_EXP 1024 311#define FLT_RADIX 2 312#endif /*IEEE_Arith*/ 313 314#ifdef IBM 315#define DBL_DIG 16 316#define DBL_MAX_10_EXP 75 317#define DBL_MAX_EXP 63 318#define FLT_RADIX 16 319#define DBL_MAX 7.2370055773322621e+75 320#endif 321 322#ifdef VAX 323#define DBL_DIG 16 324#define DBL_MAX_10_EXP 38 325#define DBL_MAX_EXP 127 326#define FLT_RADIX 2 327#define DBL_MAX 1.7014118346046923e+38 328#endif 329 330#ifndef LONG_MAX 331#define LONG_MAX 2147483647 332#endif 333 334#else /* ifndef Bad_float_h */ 335#include "float.h" 336#endif /* Bad_float_h */ 337 338#ifndef __MATH_H__ 339#include "math.h" 340#endif 341 342#ifdef __cplusplus 343extern "C" { 344#endif 345 346#ifndef CONST 347#ifdef KR_headers 348#define CONST /* blank */ 349#else 350#define CONST const 351#endif 352#endif 353 354#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 355Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. 356#endif 357 358typedef union { double d; ULong L[2]; } U; 359 360#ifdef IEEE_8087 361#define word0(x) (x)->L[1] 362#define word1(x) (x)->L[0] 363#else 364#define word0(x) (x)->L[0] 365#define word1(x) (x)->L[1] 366#endif 367#define dval(x) (x)->d 368 369#ifndef STRTOD_DIGLIM 370#define STRTOD_DIGLIM 40 371#endif 372 373#ifdef DIGLIM_DEBUG 374extern int strtod_diglim; 375#else 376#define strtod_diglim STRTOD_DIGLIM 377#endif 378 379/* The following definition of Storeinc is appropriate for MIPS processors. 380 * An alternative that might be better on some machines is 381 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 382 */ 383#if defined(IEEE_8087) + defined(VAX) 384#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 385((unsigned short *)a)[0] = (unsigned short)c, a++) 386#else 387#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 388((unsigned short *)a)[1] = (unsigned short)c, a++) 389#endif 390 391/* #define P DBL_MANT_DIG */ 392/* Ten_pmax = floor(P*log(2)/log(5)) */ 393/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 394/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 395/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 396 397#ifdef IEEE_Arith 398#define Exp_shift 20 399#define Exp_shift1 20 400#define Exp_msk1 0x100000 401#define Exp_msk11 0x100000 402#define Exp_mask 0x7ff00000 403#define P 53 404#define Nbits 53 405#define Bias 1023 406#define Emax 1023 407#define Emin (-1022) 408#define Exp_1 0x3ff00000 409#define Exp_11 0x3ff00000 410#define Ebits 11 411#define Frac_mask 0xfffff 412#define Frac_mask1 0xfffff 413#define Ten_pmax 22 414#define Bletch 0x10 415#define Bndry_mask 0xfffff 416#define Bndry_mask1 0xfffff 417#define LSB 1 418#define Sign_bit 0x80000000 419#define Log2P 1 420#define Tiny0 0 421#define Tiny1 1 422#define Quick_max 14 423#define Int_max 14 424#ifndef NO_IEEE_Scale 425#define Avoid_Underflow 426#ifdef Flush_Denorm /* debugging option */ 427#undef Sudden_Underflow 428#endif 429#endif 430 431#ifndef Flt_Rounds 432#ifdef FLT_ROUNDS 433#define Flt_Rounds FLT_ROUNDS 434#else 435#define Flt_Rounds 1 436#endif 437#endif /*Flt_Rounds*/ 438 439#ifdef Honor_FLT_ROUNDS 440#undef Check_FLT_ROUNDS 441#define Check_FLT_ROUNDS 442#else 443#define Rounding Flt_Rounds 444#endif 445 446#else /* ifndef IEEE_Arith */ 447#undef Check_FLT_ROUNDS 448#undef Honor_FLT_ROUNDS 449#undef SET_INEXACT 450#undef Sudden_Underflow 451#define Sudden_Underflow 452#ifdef IBM 453#undef Flt_Rounds 454#define Flt_Rounds 0 455#define Exp_shift 24 456#define Exp_shift1 24 457#define Exp_msk1 0x1000000 458#define Exp_msk11 0x1000000 459#define Exp_mask 0x7f000000 460#define P 14 461#define Nbits 56 462#define Bias 65 463#define Emax 248 464#define Emin (-260) 465#define Exp_1 0x41000000 466#define Exp_11 0x41000000 467#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 468#define Frac_mask 0xffffff 469#define Frac_mask1 0xffffff 470#define Bletch 4 471#define Ten_pmax 22 472#define Bndry_mask 0xefffff 473#define Bndry_mask1 0xffffff 474#define LSB 1 475#define Sign_bit 0x80000000 476#define Log2P 4 477#define Tiny0 0x100000 478#define Tiny1 0 479#define Quick_max 14 480#define Int_max 15 481#else /* VAX */ 482#undef Flt_Rounds 483#define Flt_Rounds 1 484#define Exp_shift 23 485#define Exp_shift1 7 486#define Exp_msk1 0x80 487#define Exp_msk11 0x800000 488#define Exp_mask 0x7f80 489#define P 56 490#define Nbits 56 491#define Bias 129 492#define Emax 126 493#define Emin (-129) 494#define Exp_1 0x40800000 495#define Exp_11 0x4080 496#define Ebits 8 497#define Frac_mask 0x7fffff 498#define Frac_mask1 0xffff007f 499#define Ten_pmax 24 500#define Bletch 2 501#define Bndry_mask 0xffff007f 502#define Bndry_mask1 0xffff007f 503#define LSB 0x10000 504#define Sign_bit 0x8000 505#define Log2P 1 506#define Tiny0 0x80 507#define Tiny1 0 508#define Quick_max 15 509#define Int_max 15 510#endif /* IBM, VAX */ 511#endif /* IEEE_Arith */ 512 513#ifndef IEEE_Arith 514#define ROUND_BIASED 515#else 516#ifdef ROUND_BIASED_without_Round_Up 517#undef ROUND_BIASED 518#define ROUND_BIASED 519#endif 520#endif 521 522#ifdef RND_PRODQUOT 523#define rounded_product(a,b) a = rnd_prod(a, b) 524#define rounded_quotient(a,b) a = rnd_quot(a, b) 525#ifdef KR_headers 526extern double rnd_prod(), rnd_quot(); 527#else 528extern double rnd_prod(double, double), rnd_quot(double, double); 529#endif 530#else 531#define rounded_product(a,b) a *= b 532#define rounded_quotient(a,b) a /= b 533#endif 534 535#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 536#define Big1 0xffffffff 537 538#ifndef Pack_32 539#define Pack_32 540#endif 541 542typedef struct BCinfo BCinfo; 543 struct 544BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflchk; }; 545 546#ifdef KR_headers 547#define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff) 548#else 549#define FFFFFFFF 0xffffffffUL 550#endif 551 552#ifdef NO_LONG_LONG 553#undef ULLong 554#ifdef Just_16 555#undef Pack_32 556/* When Pack_32 is not defined, we store 16 bits per 32-bit Long. 557 * This makes some inner loops simpler and sometimes saves work 558 * during multiplications, but it often seems to make things slightly 559 * slower. Hence the default is now to store 32 bits per Long. 560 */ 561#endif 562#else /* long long available */ 563#ifndef Llong 564#define Llong long long 565#endif 566#ifndef ULLong 567#define ULLong unsigned Llong 568#endif 569#endif /* NO_LONG_LONG */ 570 571#define MULTIPLE_THREADS 572static PLock dtoaLocks[2]; 573#define ACQUIRE_DTOA_LOCK(n) { dtoaLocks[n].Lock(); } 574#define FREE_DTOA_LOCK(n) { dtoaLocks[n].Unlock(); } 575 576#ifndef MULTIPLE_THREADS 577#define ACQUIRE_DTOA_LOCK(n) /*nothing*/ 578#define FREE_DTOA_LOCK(n) /*nothing*/ 579#endif 580 581#define Kmax 7 582 583#ifdef __cplusplus 584extern "C" double strtod(const char *s00, char **se); 585extern "C" char *dtoa(double d, int mode, int ndigits, 586 int *decpt, int *sign, char **rve); 587#endif 588 589 struct 590Bigint { 591 struct Bigint *next; 592 int k, maxwds, sign, wds; 593 ULong x[1]; 594 }; 595 596 typedef struct Bigint Bigint; 597 598 static Bigint *freelist[Kmax+1]; 599 600 static Bigint * 601Balloc 602#ifdef KR_headers 603 (k) int k; 604#else 605 (int k) 606#endif 607{ 608 int x; 609 Bigint *rv; 610#ifndef Omit_Private_Memory 611 unsigned int len; 612#endif 613 614 ACQUIRE_DTOA_LOCK(0); 615 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */ 616 /* but this case seems very unlikely. */ 617 if (k <= Kmax && (rv = freelist[k])) 618 freelist[k] = rv->next; 619 else { 620 x = 1 << k; 621#ifdef Omit_Private_Memory 622 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); 623#else 624 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 625 /sizeof(double); 626 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) { 627 rv = (Bigint*)pmem_next; 628 pmem_next += len; 629 } 630 else 631 rv = (Bigint*)MALLOC(len*sizeof(double)); 632#endif 633 rv->k = k; 634 rv->maxwds = x; 635 } 636 FREE_DTOA_LOCK(0); 637 rv->sign = rv->wds = 0; 638 return rv; 639 } 640 641 static void 642Bfree 643#ifdef KR_headers 644 (v) Bigint *v; 645#else 646 (Bigint *v) 647#endif 648{ 649 if (v) { 650 if (v->k > Kmax) 651#ifdef FREE 652 FREE((void*)v); 653#else 654 free((void*)v); 655#endif 656 else { 657 ACQUIRE_DTOA_LOCK(0); 658 v->next = freelist[v->k]; 659 freelist[v->k] = v; 660 FREE_DTOA_LOCK(0); 661 } 662 } 663 } 664 665#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 666y->wds*sizeof(Long) + 2*sizeof(int)) 667 668 static Bigint * 669multadd 670#ifdef KR_headers 671 (b, m, a) Bigint *b; int m, a; 672#else 673 (Bigint *b, int m, int a) /* multiply by m and add a */ 674#endif 675{ 676 int i, wds; 677#ifdef ULLong 678 ULong *x; 679 ULLong carry, y; 680#else 681 ULong carry, *x, y; 682#ifdef Pack_32 683 ULong xi, z; 684#endif 685#endif 686 Bigint *b1; 687 688 wds = b->wds; 689 x = b->x; 690 i = 0; 691 carry = a; 692 do { 693#ifdef ULLong 694 y = *x * (ULLong)m + carry; 695 carry = y >> 32; 696 *x++ = y & FFFFFFFF; 697#else 698#ifdef Pack_32 699 xi = *x; 700 y = (xi & 0xffff) * m + carry; 701 z = (xi >> 16) * m + (y >> 16); 702 carry = z >> 16; 703 *x++ = (z << 16) + (y & 0xffff); 704#else 705 y = *x * m + carry; 706 carry = y >> 16; 707 *x++ = y & 0xffff; 708#endif 709#endif 710 } 711 while(++i < wds); 712 if (carry) { 713 if (wds >= b->maxwds) { 714 b1 = Balloc(b->k+1); 715 Bcopy(b1, b); 716 Bfree(b); 717 b = b1; 718 } 719 b->x[wds++] = carry; 720 b->wds = wds; 721 } 722 return b; 723 } 724 725#ifndef HAVE_STRTOD 726 static Bigint * 727s2b 728#ifdef KR_headers 729 (s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9; 730#else 731 (const char *s, int nd0, int nd, ULong y9, int dplen) 732#endif 733{ 734 Bigint *b; 735 int i, k; 736 Long x, y; 737 738 x = (nd + 8) / 9; 739 for(k = 0, y = 1; x > y; y <<= 1, k++) ; 740#ifdef Pack_32 741 b = Balloc(k); 742 b->x[0] = y9; 743 b->wds = 1; 744#else 745 b = Balloc(k+1); 746 b->x[0] = y9 & 0xffff; 747 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; 748#endif 749 750 i = 9; 751 if (9 < nd0) { 752 s += 9; 753 do b = multadd(b, 10, *s++ - '0'); 754 while(++i < nd0); 755 s += dplen; 756 } 757 else 758 s += dplen + 9; 759 for(; i < nd; i++) 760 b = multadd(b, 10, *s++ - '0'); 761 return b; 762 } 763 764#endif // HAVE_STRTOD 765 766 static int 767hi0bits 768#ifdef KR_headers 769 (x) ULong x; 770#else 771 (ULong x) 772#endif 773{ 774 int k = 0; 775 776 if (!(x & 0xffff0000)) { 777 k = 16; 778 x <<= 16; 779 } 780 if (!(x & 0xff000000)) { 781 k += 8; 782 x <<= 8; 783 } 784 if (!(x & 0xf0000000)) { 785 k += 4; 786 x <<= 4; 787 } 788 if (!(x & 0xc0000000)) { 789 k += 2; 790 x <<= 2; 791 } 792 if (!(x & 0x80000000)) { 793 k++; 794 if (!(x & 0x40000000)) 795 return 32; 796 } 797 return k; 798 } 799 800 static int 801lo0bits 802#ifdef KR_headers 803 (y) ULong *y; 804#else 805 (ULong *y) 806#endif 807{ 808 int k; 809 ULong x = *y; 810 811 if (x & 7) { 812 if (x & 1) 813 return 0; 814 if (x & 2) { 815 *y = x >> 1; 816 return 1; 817 } 818 *y = x >> 2; 819 return 2; 820 } 821 k = 0; 822 if (!(x & 0xffff)) { 823 k = 16; 824 x >>= 16; 825 } 826 if (!(x & 0xff)) { 827 k += 8; 828 x >>= 8; 829 } 830 if (!(x & 0xf)) { 831 k += 4; 832 x >>= 4; 833 } 834 if (!(x & 0x3)) { 835 k += 2; 836 x >>= 2; 837 } 838 if (!(x & 1)) { 839 k++; 840 x >>= 1; 841 if (!x) 842 return 32; 843 } 844 *y = x; 845 return k; 846 } 847 848 static Bigint * 849i2b 850#ifdef KR_headers 851 (i) int i; 852#else 853 (int i) 854#endif 855{ 856 Bigint *b; 857 858 b = Balloc(1); 859 b->x[0] = i; 860 b->wds = 1; 861 return b; 862 } 863 864 static Bigint * 865mult 866#ifdef KR_headers 867 (a, b) Bigint *a, *b; 868#else 869 (Bigint *a, Bigint *b) 870#endif 871{ 872 Bigint *c; 873 int k, wa, wb, wc; 874 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 875 ULong y; 876#ifdef ULLong 877 ULLong carry, z; 878#else 879 ULong carry, z; 880#ifdef Pack_32 881 ULong z2; 882#endif 883#endif 884 885 if (a->wds < b->wds) { 886 c = a; 887 a = b; 888 b = c; 889 } 890 k = a->k; 891 wa = a->wds; 892 wb = b->wds; 893 wc = wa + wb; 894 if (wc > a->maxwds) 895 k++; 896 c = Balloc(k); 897 for(x = c->x, xa = x + wc; x < xa; x++) 898 *x = 0; 899 xa = a->x; 900 xae = xa + wa; 901 xb = b->x; 902 xbe = xb + wb; 903 xc0 = c->x; 904#ifdef ULLong 905 for(; xb < xbe; xc0++) { 906 if ((y = *xb++)) { 907 x = xa; 908 xc = xc0; 909 carry = 0; 910 do { 911 z = *x++ * (ULLong)y + *xc + carry; 912 carry = z >> 32; 913 *xc++ = z & FFFFFFFF; 914 } 915 while(x < xae); 916 *xc = carry; 917 } 918 } 919#else 920#ifdef Pack_32 921 for(; xb < xbe; xb++, xc0++) { 922 if (y = *xb & 0xffff) { 923 x = xa; 924 xc = xc0; 925 carry = 0; 926 do { 927 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 928 carry = z >> 16; 929 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 930 carry = z2 >> 16; 931 Storeinc(xc, z2, z); 932 } 933 while(x < xae); 934 *xc = carry; 935 } 936 if (y = *xb >> 16) { 937 x = xa; 938 xc = xc0; 939 carry = 0; 940 z2 = *xc; 941 do { 942 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 943 carry = z >> 16; 944 Storeinc(xc, z, z2); 945 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 946 carry = z2 >> 16; 947 } 948 while(x < xae); 949 *xc = z2; 950 } 951 } 952#else 953 for(; xb < xbe; xc0++) { 954 if (y = *xb++) { 955 x = xa; 956 xc = xc0; 957 carry = 0; 958 do { 959 z = *x++ * y + *xc + carry; 960 carry = z >> 16; 961 *xc++ = z & 0xffff; 962 } 963 while(x < xae); 964 *xc = carry; 965 } 966 } 967#endif 968#endif 969 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 970 c->wds = wc; 971 return c; 972 } 973 974 static Bigint *p5s; 975 976 static Bigint * 977pow5mult 978#ifdef KR_headers 979 (b, k) Bigint *b; int k; 980#else 981 (Bigint *b, int k) 982#endif 983{ 984 Bigint *b1, *p5, *p51; 985 int i; 986 static int p05[3] = { 5, 25, 125 }; 987 988 if ((i = k & 3)) 989 b = multadd(b, p05[i-1], 0); 990 991 if (!(k >>= 2)) 992 return b; 993 if (!(p5 = p5s)) { 994 /* first time */ 995#ifdef MULTIPLE_THREADS 996 ACQUIRE_DTOA_LOCK(1); 997 if (!(p5 = p5s)) { 998 p5 = p5s = i2b(625); 999 p5->next = 0; 1000 } 1001 FREE_DTOA_LOCK(1); 1002#else 1003 p5 = p5s = i2b(625); 1004 p5->next = 0; 1005#endif 1006 } 1007 for(;;) { 1008 if (k & 1) { 1009 b1 = mult(b, p5); 1010 Bfree(b); 1011 b = b1; 1012 } 1013 if (!(k >>= 1)) 1014 break; 1015 if (!(p51 = p5->next)) { 1016#ifdef MULTIPLE_THREADS 1017 ACQUIRE_DTOA_LOCK(1); 1018 if (!(p51 = p5->next)) { 1019 p51 = p5->next = mult(p5,p5); 1020 p51->next = 0; 1021 } 1022 FREE_DTOA_LOCK(1); 1023#else 1024 p51 = p5->next = mult(p5,p5); 1025 p51->next = 0; 1026#endif 1027 } 1028 p5 = p51; 1029 } 1030 return b; 1031 } 1032 1033 static Bigint * 1034lshift 1035#ifdef KR_headers 1036 (b, k) Bigint *b; int k; 1037#else 1038 (Bigint *b, int k) 1039#endif 1040{ 1041 int i, k1, n, n1; 1042 Bigint *b1; 1043 ULong *x, *x1, *xe, z; 1044 1045#ifdef Pack_32 1046 n = k >> 5; 1047#else 1048 n = k >> 4; 1049#endif 1050 k1 = b->k; 1051 n1 = n + b->wds + 1; 1052 for(i = b->maxwds; n1 > i; i <<= 1) 1053 k1++; 1054 b1 = Balloc(k1); 1055 x1 = b1->x; 1056 for(i = 0; i < n; i++) 1057 *x1++ = 0; 1058 x = b->x; 1059 xe = x + b->wds; 1060#ifdef Pack_32 1061 if (k &= 0x1f) { 1062 k1 = 32 - k; 1063 z = 0; 1064 do { 1065 *x1++ = *x << k | z; 1066 z = *x++ >> k1; 1067 } 1068 while(x < xe); 1069 if ((*x1 = z)) 1070 ++n1; 1071 } 1072#else 1073 if (k &= 0xf) { 1074 k1 = 16 - k; 1075 z = 0; 1076 do { 1077 *x1++ = *x << k & 0xffff | z; 1078 z = *x++ >> k1; 1079 } 1080 while(x < xe); 1081 if (*x1 = z) 1082 ++n1; 1083 } 1084#endif 1085 else do 1086 *x1++ = *x++; 1087 while(x < xe); 1088 b1->wds = n1 - 1; 1089 Bfree(b); 1090 return b1; 1091 } 1092 1093 static int 1094cmp 1095#ifdef KR_headers 1096 (a, b) Bigint *a, *b; 1097#else 1098 (Bigint *a, Bigint *b) 1099#endif 1100{ 1101 ULong *xa, *xa0, *xb, *xb0; 1102 int i, j; 1103 1104 i = a->wds; 1105 j = b->wds; 1106#ifdef DEBUG 1107 if (i > 1 && !a->x[i-1]) 1108 Bug("cmp called with a->x[a->wds-1] == 0"); 1109 if (j > 1 && !b->x[j-1]) 1110 Bug("cmp called with b->x[b->wds-1] == 0"); 1111#endif 1112 if (i -= j) 1113 return i; 1114 xa0 = a->x; 1115 xa = xa0 + j; 1116 xb0 = b->x; 1117 xb = xb0 + j; 1118 for(;;) { 1119 if (*--xa != *--xb) 1120 return *xa < *xb ? -1 : 1; 1121 if (xa <= xa0) 1122 break; 1123 } 1124 return 0; 1125 } 1126 1127 static Bigint * 1128diff 1129#ifdef KR_headers 1130 (a, b) Bigint *a, *b; 1131#else 1132 (Bigint *a, Bigint *b) 1133#endif 1134{ 1135 Bigint *c; 1136 int i, wa, wb; 1137 ULong *xa, *xae, *xb, *xbe, *xc; 1138#ifdef ULLong 1139 ULLong borrow, y; 1140#else 1141 ULong borrow, y; 1142#ifdef Pack_32 1143 ULong z; 1144#endif 1145#endif 1146 1147 i = cmp(a,b); 1148 if (!i) { 1149 c = Balloc(0); 1150 c->wds = 1; 1151 c->x[0] = 0; 1152 return c; 1153 } 1154 if (i < 0) { 1155 c = a; 1156 a = b; 1157 b = c; 1158 i = 1; 1159 } 1160 else 1161 i = 0; 1162 c = Balloc(a->k); 1163 c->sign = i; 1164 wa = a->wds; 1165 xa = a->x; 1166 xae = xa + wa; 1167 wb = b->wds; 1168 xb = b->x; 1169 xbe = xb + wb; 1170 xc = c->x; 1171 borrow = 0; 1172#ifdef ULLong 1173 do { 1174 y = (ULLong)*xa++ - *xb++ - borrow; 1175 borrow = y >> 32 & (ULong)1; 1176 *xc++ = y & FFFFFFFF; 1177 } 1178 while(xb < xbe); 1179 while(xa < xae) { 1180 y = *xa++ - borrow; 1181 borrow = y >> 32 & (ULong)1; 1182 *xc++ = y & FFFFFFFF; 1183 } 1184#else 1185#ifdef Pack_32 1186 do { 1187 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 1188 borrow = (y & 0x10000) >> 16; 1189 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 1190 borrow = (z & 0x10000) >> 16; 1191 Storeinc(xc, z, y); 1192 } 1193 while(xb < xbe); 1194 while(xa < xae) { 1195 y = (*xa & 0xffff) - borrow; 1196 borrow = (y & 0x10000) >> 16; 1197 z = (*xa++ >> 16) - borrow; 1198 borrow = (z & 0x10000) >> 16; 1199 Storeinc(xc, z, y); 1200 } 1201#else 1202 do { 1203 y = *xa++ - *xb++ - borrow; 1204 borrow = (y & 0x10000) >> 16; 1205 *xc++ = y & 0xffff; 1206 } 1207 while(xb < xbe); 1208 while(xa < xae) { 1209 y = *xa++ - borrow; 1210 borrow = (y & 0x10000) >> 16; 1211 *xc++ = y & 0xffff; 1212 } 1213#endif 1214#endif 1215 while(!*--xc) 1216 wa--; 1217 c->wds = wa; 1218 return c; 1219 } 1220 1221#ifndef HAVE_STRTOD 1222 1223 static double 1224ulp 1225#ifdef KR_headers 1226 (x) U *x; 1227#else 1228 (U *x) 1229#endif 1230{ 1231 Long L; 1232 U u; 1233 1234 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 1235#ifndef Avoid_Underflow 1236#ifndef Sudden_Underflow 1237 if (L > 0) { 1238#endif 1239#endif 1240#ifdef IBM 1241 L |= Exp_msk1 >> 4; 1242#endif 1243 word0(&u) = L; 1244 word1(&u) = 0; 1245#ifndef Avoid_Underflow 1246#ifndef Sudden_Underflow 1247 } 1248 else { 1249 L = -L >> Exp_shift; 1250 if (L < Exp_shift) { 1251 word0(&u) = 0x80000 >> L; 1252 word1(&u) = 0; 1253 } 1254 else { 1255 word0(&u) = 0; 1256 L -= Exp_shift; 1257 word1(&u) = L >= 31 ? 1 : 1 << 31 - L; 1258 } 1259 } 1260#endif 1261#endif 1262 return dval(&u); 1263 } 1264 1265 static double 1266b2d 1267#ifdef KR_headers 1268 (a, e) Bigint *a; int *e; 1269#else 1270 (Bigint *a, int *e) 1271#endif 1272{ 1273 ULong *xa, *xa0, w, y, z; 1274 int k; 1275 U d; 1276#ifdef VAX 1277 ULong d0, d1; 1278#else 1279#define d0 word0(&d) 1280#define d1 word1(&d) 1281#endif 1282 1283 xa0 = a->x; 1284 xa = xa0 + a->wds; 1285 y = *--xa; 1286#ifdef DEBUG 1287 if (!y) Bug("zero y in b2d"); 1288#endif 1289 k = hi0bits(y); 1290 *e = 32 - k; 1291#ifdef Pack_32 1292 if (k < Ebits) { 1293 d0 = Exp_1 | y >> (Ebits - k); 1294 w = xa > xa0 ? *--xa : 0; 1295 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k); 1296 goto ret_d; 1297 } 1298 z = xa > xa0 ? *--xa : 0; 1299 if (k -= Ebits) { 1300 d0 = Exp_1 | y << k | z >> (32 - k); 1301 y = xa > xa0 ? *--xa : 0; 1302 d1 = z << k | y >> (32 - k); 1303 } 1304 else { 1305 d0 = Exp_1 | y; 1306 d1 = z; 1307 } 1308#else 1309 if (k < Ebits + 16) { 1310 z = xa > xa0 ? *--xa : 0; 1311 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 1312 w = xa > xa0 ? *--xa : 0; 1313 y = xa > xa0 ? *--xa : 0; 1314 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 1315 goto ret_d; 1316 } 1317 z = xa > xa0 ? *--xa : 0; 1318 w = xa > xa0 ? *--xa : 0; 1319 k -= Ebits + 16; 1320 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 1321 y = xa > xa0 ? *--xa : 0; 1322 d1 = w << k + 16 | y << k; 1323#endif 1324 ret_d: 1325#ifdef VAX 1326 word0(&d) = d0 >> 16 | d0 << 16; 1327 word1(&d) = d1 >> 16 | d1 << 16; 1328#else 1329#undef d0 1330#undef d1 1331#endif 1332 return dval(&d); 1333 } 1334 1335#endif // HAVE_STRTOD 1336 1337 static Bigint * 1338d2b 1339#ifdef KR_headers 1340 (d, e, bits) U *d; int *e, *bits; 1341#else 1342 (U *d, int *e, int *bits) 1343#endif 1344{ 1345 Bigint *b; 1346 int de, k; 1347 ULong *x, y, z; 1348#ifndef Sudden_Underflow 1349 int i; 1350#endif 1351#ifdef VAX 1352 ULong d0, d1; 1353 d0 = word0(d) >> 16 | word0(d) << 16; 1354 d1 = word1(d) >> 16 | word1(d) << 16; 1355#else 1356#define d0 word0(d) 1357#define d1 word1(d) 1358#endif 1359 1360#ifdef Pack_32 1361 b = Balloc(1); 1362#else 1363 b = Balloc(2); 1364#endif 1365 x = b->x; 1366 1367 z = d0 & Frac_mask; 1368 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 1369#ifdef Sudden_Underflow 1370 de = (int)(d0 >> Exp_shift); 1371#ifndef IBM 1372 z |= Exp_msk11; 1373#endif 1374#else 1375 if ((de = (int)(d0 >> Exp_shift))) 1376 z |= Exp_msk1; 1377#endif 1378#ifdef Pack_32 1379 if ((y = d1)) { 1380 if ((k = lo0bits(&y))) { 1381 x[0] = y | z << (32 - k); 1382 z >>= k; 1383 } 1384 else 1385 x[0] = y; 1386#ifndef Sudden_Underflow 1387 i = 1388#endif 1389 b->wds = (x[1] = z) ? 2 : 1; 1390 } 1391 else { 1392 k = lo0bits(&z); 1393 x[0] = z; 1394#ifndef Sudden_Underflow 1395 i = 1396#endif 1397 b->wds = 1; 1398 k += 32; 1399 } 1400#else 1401 if (y = d1) { 1402 if (k = lo0bits(&y)) 1403 if (k >= 16) { 1404 x[0] = y | z << 32 - k & 0xffff; 1405 x[1] = z >> k - 16 & 0xffff; 1406 x[2] = z >> k; 1407 i = 2; 1408 } 1409 else { 1410 x[0] = y & 0xffff; 1411 x[1] = y >> 16 | z << 16 - k & 0xffff; 1412 x[2] = z >> k & 0xffff; 1413 x[3] = z >> k+16; 1414 i = 3; 1415 } 1416 else { 1417 x[0] = y & 0xffff; 1418 x[1] = y >> 16; 1419 x[2] = z & 0xffff; 1420 x[3] = z >> 16; 1421 i = 3; 1422 } 1423 } 1424 else { 1425#ifdef DEBUG 1426 if (!z) 1427 Bug("Zero passed to d2b"); 1428#endif 1429 k = lo0bits(&z); 1430 if (k >= 16) { 1431 x[0] = z; 1432 i = 0; 1433 } 1434 else { 1435 x[0] = z & 0xffff; 1436 x[1] = z >> 16; 1437 i = 1; 1438 } 1439 k += 32; 1440 } 1441 while(!x[i]) 1442 --i; 1443 b->wds = i + 1; 1444#endif 1445#ifndef Sudden_Underflow 1446 if (de) { 1447#endif 1448#ifdef IBM 1449 *e = (de - Bias - (P-1) << 2) + k; 1450 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 1451#else 1452 *e = de - Bias - (P-1) + k; 1453 *bits = P - k; 1454#endif 1455#ifndef Sudden_Underflow 1456 } 1457 else { 1458 *e = de - Bias - (P-1) + 1 + k; 1459#ifdef Pack_32 1460 *bits = 32*i - hi0bits(x[i-1]); 1461#else 1462 *bits = (i+2)*16 - hi0bits(x[i]); 1463#endif 1464 } 1465#endif 1466 return b; 1467 } 1468#undef d0 1469#undef d1 1470 1471#ifndef HAVE_STRTOD 1472 static double 1473ratio 1474#ifdef KR_headers 1475 (a, b) Bigint *a, *b; 1476#else 1477 (Bigint *a, Bigint *b) 1478#endif 1479{ 1480 U da, db; 1481 int k, ka, kb; 1482 1483 dval(&da) = b2d(a, &ka); 1484 dval(&db) = b2d(b, &kb); 1485#ifdef Pack_32 1486 k = ka - kb + 32*(a->wds - b->wds); 1487#else 1488 k = ka - kb + 16*(a->wds - b->wds); 1489#endif 1490#ifdef IBM 1491 if (k > 0) { 1492 word0(&da) += (k >> 2)*Exp_msk1; 1493 if (k &= 3) 1494 dval(&da) *= 1 << k; 1495 } 1496 else { 1497 k = -k; 1498 word0(&db) += (k >> 2)*Exp_msk1; 1499 if (k &= 3) 1500 dval(&db) *= 1 << k; 1501 } 1502#else 1503 if (k > 0) 1504 word0(&da) += k*Exp_msk1; 1505 else { 1506 k = -k; 1507 word0(&db) += k*Exp_msk1; 1508 } 1509#endif 1510 return dval(&da) / dval(&db); 1511 } 1512#endif // HAVE_STRTOD 1513 1514 static CONST double 1515tens[] = { 1516 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1517 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1518 1e20, 1e21, 1e22 1519#ifdef VAX 1520 , 1e23, 1e24 1521#endif 1522 }; 1523 1524 static CONST double 1525#ifdef IEEE_Arith 1526bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1527static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1528#ifdef Avoid_Underflow 1529 9007199254740992.*9007199254740992.e-256 1530 /* = 2^106 * 1e-256 */ 1531#else 1532 1e-256 1533#endif 1534 }; 1535/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 1536/* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 1537#define Scale_Bit 0x10 1538#define n_bigtens 5 1539#else 1540#ifdef IBM 1541bigtens[] = { 1e16, 1e32, 1e64 }; 1542static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 1543#define n_bigtens 3 1544#else 1545bigtens[] = { 1e16, 1e32 }; 1546static CONST double tinytens[] = { 1e-16, 1e-32 }; 1547#define n_bigtens 2 1548#endif 1549#endif 1550 1551#undef Need_Hexdig 1552#ifdef INFNAN_CHECK 1553#ifndef No_Hex_NaN 1554#define Need_Hexdig 1555#endif 1556#endif 1557 1558#ifndef Need_Hexdig 1559#ifndef NO_HEX_FP 1560#define Need_Hexdig 1561#endif 1562#endif 1563 1564#ifdef Need_Hexdig /*{*/ 1565#if 0 1566static unsigned char hexdig[256]; 1567 1568 static void 1569htinit(unsigned char *h, unsigned char *s, int inc) 1570{ 1571 int i, j; 1572 for(i = 0; (j = s[i]) !=0; i++) 1573 h[j] = i + inc; 1574 } 1575 1576 static void 1577hexdig_init(void) /* Use of hexdig_init omitted 20121220 to avoid a */ 1578 /* race condition when multiple threads are used. */ 1579{ 1580#define USC (unsigned char *) 1581 htinit(hexdig, USC "0123456789", 0x10); 1582 htinit(hexdig, USC "abcdef", 0x10 + 10); 1583 htinit(hexdig, USC "ABCDEF", 0x10 + 10); 1584 } 1585#else 1586static unsigned char hexdig[256] = { 1587 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1588 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1589 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1590 16,17,18,19,20,21,22,23,24,25,0,0,0,0,0,0, 1591 0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0, 1592 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1593 0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0, 1594 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1595 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1596 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1597 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1598 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1599 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1600 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1601 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1602 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 1603 }; 1604#endif 1605#endif /* } Need_Hexdig */ 1606 1607#ifdef INFNAN_CHECK 1608 1609#ifndef NAN_WORD0 1610#define NAN_WORD0 0x7ff80000 1611#endif 1612 1613#ifndef NAN_WORD1 1614#define NAN_WORD1 0 1615#endif 1616 1617#ifndef HAVE_STRTOD 1618 static int 1619match 1620#ifdef KR_headers 1621 (sp, t) char **sp, *t; 1622#else 1623 (const char **sp, const char *t) 1624#endif 1625{ 1626 int c, d; 1627 CONST char *s = *sp; 1628 1629 while((d = *t++)) { 1630 if ((c = *++s) >= 'A' && c <= 'Z') 1631 c += 'a' - 'A'; 1632 if (c != d) 1633 return 0; 1634 } 1635 *sp = s + 1; 1636 return 1; 1637 } 1638 1639#ifndef No_Hex_NaN 1640 static void 1641hexnan 1642#ifdef KR_headers 1643 (rvp, sp) U *rvp; CONST char **sp; 1644#else 1645 (U *rvp, const char **sp) 1646#endif 1647{ 1648 ULong c, x[2]; 1649 CONST char *s; 1650 int c1, havedig, udx0, xshift; 1651 1652 /**** if (!hexdig['0']) hexdig_init(); ****/ 1653 x[0] = x[1] = 0; 1654 havedig = xshift = 0; 1655 udx0 = 1; 1656 s = *sp; 1657 /* allow optional initial 0x or 0X */ 1658 while((c = *(CONST unsigned char*)(s+1)) && c <= ' ') 1659 ++s; 1660 if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')) 1661 s += 2; 1662 while((c = *(CONST unsigned char*)++s)) { 1663 if ((c1 = hexdig[c])) 1664 c = c1 & 0xf; 1665 else if (c <= ' ') { 1666 if (udx0 && havedig) { 1667 udx0 = 0; 1668 xshift = 1; 1669 } 1670 continue; 1671 } 1672#ifdef GDTOA_NON_PEDANTIC_NANCHECK 1673 else if (/*(*/ c == ')' && havedig) { 1674 *sp = s + 1; 1675 break; 1676 } 1677 else 1678 return; /* invalid form: don't change *sp */ 1679#else 1680 else { 1681 do { 1682 if (/*(*/ c == ')') { 1683 *sp = s + 1; 1684 break; 1685 } 1686 } while((c = *++s)); 1687 break; 1688 } 1689#endif 1690 havedig = 1; 1691 if (xshift) { 1692 xshift = 0; 1693 x[0] = x[1]; 1694 x[1] = 0; 1695 } 1696 if (udx0) 1697 x[0] = (x[0] << 4) | (x[1] >> 28); 1698 x[1] = (x[1] << 4) | c; 1699 } 1700 if ((x[0] &= 0xfffff) || x[1]) { 1701 word0(rvp) = Exp_mask | x[0]; 1702 word1(rvp) = x[1]; 1703 } 1704 } 1705#endif /*No_Hex_NaN*/ 1706#endif /* INFNAN_CHECK */ 1707 1708#endif // HAVE_STRTOD 1709 1710#ifdef Pack_32 1711#define ULbits 32 1712#define kshift 5 1713#define kmask 31 1714#else 1715#define ULbits 16 1716#define kshift 4 1717#define kmask 15 1718#endif 1719 1720#if !defined(NO_HEX_FP) || defined(Honor_FLT_ROUNDS) /*{*/ 1721 static Bigint * 1722#ifdef KR_headers 1723increment(b) Bigint *b; 1724#else 1725increment(Bigint *b) 1726#endif 1727{ 1728 ULong *x, *xe; 1729 Bigint *b1; 1730 1731 x = b->x; 1732 xe = x + b->wds; 1733 do { 1734 if (*x < (ULong)0xffffffffL) { 1735 ++*x; 1736 return b; 1737 } 1738 *x++ = 0; 1739 } while(x < xe); 1740 { 1741 if (b->wds >= b->maxwds) { 1742 b1 = Balloc(b->k+1); 1743 Bcopy(b1,b); 1744 Bfree(b); 1745 b = b1; 1746 } 1747 b->x[b->wds++] = 1; 1748 } 1749 return b; 1750 } 1751 1752#endif /*}*/ 1753 1754#ifndef NO_HEX_FP /*{*/ 1755 1756 static void 1757#ifdef KR_headers 1758rshift(b, k) Bigint *b; int k; 1759#else 1760rshift(Bigint *b, int k) 1761#endif 1762{ 1763 ULong *x, *x1, *xe, y; 1764 int n; 1765 1766 x = x1 = b->x; 1767 n = k >> kshift; 1768 if (n < b->wds) { 1769 xe = x + b->wds; 1770 x += n; 1771 if (k &= kmask) { 1772 n = 32 - k; 1773 y = *x++ >> k; 1774 while(x < xe) { 1775 *x1++ = (y | (*x << n)) & 0xffffffff; 1776 y = *x++ >> k; 1777 } 1778 if ((*x1 = y) !=0) 1779 x1++; 1780 } 1781 else 1782 while(x < xe) 1783 *x1++ = *x++; 1784 } 1785 if ((b->wds = x1 - b->x) == 0) 1786 b->x[0] = 0; 1787 } 1788 1789 static ULong 1790#ifdef KR_headers 1791any_on(b, k) Bigint *b; int k; 1792#else 1793any_on(Bigint *b, int k) 1794#endif 1795{ 1796 int n, nwds; 1797 ULong *x, *x0, x1, x2; 1798 1799 x = b->x; 1800 nwds = b->wds; 1801 n = k >> kshift; 1802 if (n > nwds) 1803 n = nwds; 1804 else if (n < nwds && (k &= kmask)) { 1805 x1 = x2 = x[n]; 1806 x1 >>= k; 1807 x1 <<= k; 1808 if (x1 != x2) 1809 return 1; 1810 } 1811 x0 = x; 1812 x += n; 1813 while(x > x0) 1814 if (*--x) 1815 return 1; 1816 return 0; 1817 } 1818 1819enum { /* rounding values: same as FLT_ROUNDS */ 1820 Round_zero = 0, 1821 Round_near = 1, 1822 Round_up = 2, 1823 Round_down = 3 1824 }; 1825 1826 void 1827#ifdef KR_headers 1828gethex(sp, rvp, rounding, sign) 1829 CONST char **sp; U *rvp; int rounding, sign; 1830#else 1831gethex( CONST char **sp, U *rvp, int rounding, int sign) 1832#endif 1833{ 1834 Bigint *b; 1835 CONST unsigned char *decpt, *s0, *s, *s1; 1836 Long e, e1; 1837 ULong L, lostbits, *x; 1838 int big, denorm, esign, havedig, k, n, nbits, up, zret; 1839#ifdef IBM 1840 int j; 1841#endif 1842 enum { 1843#ifdef IEEE_Arith /*{{*/ 1844 emax = 0x7fe - Bias - P + 1, 1845 emin = Emin - P + 1 1846#else /*}{*/ 1847 emin = Emin - P, 1848#ifdef VAX 1849 emax = 0x7ff - Bias - P + 1 1850#endif 1851#ifdef IBM 1852 emax = 0x7f - Bias - P 1853#endif 1854#endif /*}}*/ 1855 }; 1856#ifdef USE_LOCALE 1857 int i; 1858#ifdef NO_LOCALE_CACHE 1859 const unsigned char *decimalpoint = (unsigned char*) 1860 localeconv()->decimal_point; 1861#else 1862 const unsigned char *decimalpoint; 1863 static unsigned char *decimalpoint_cache; 1864 if (!(s0 = decimalpoint_cache)) { 1865 s0 = (unsigned char*)localeconv()->decimal_point; 1866 if ((decimalpoint_cache = (unsigned char*) 1867 MALLOC(strlen((CONST char*)s0) + 1))) { 1868 strcpy((char*)decimalpoint_cache, (CONST char*)s0); 1869 s0 = decimalpoint_cache; 1870 } 1871 } 1872 decimalpoint = s0; 1873#endif 1874#endif 1875 1876 /**** if (!hexdig['0']) hexdig_init(); ****/ 1877 havedig = 0; 1878 s0 = *(CONST unsigned char **)sp + 2; 1879 while(s0[havedig] == '0') 1880 havedig++; 1881 s0 += havedig; 1882 s = s0; 1883 decpt = 0; 1884 zret = 0; 1885 e = 0; 1886 if (hexdig[*s]) 1887 havedig++; 1888 else { 1889 zret = 1; 1890#ifdef USE_LOCALE 1891 for(i = 0; decimalpoint[i]; ++i) { 1892 if (s[i] != decimalpoint[i]) 1893 goto pcheck; 1894 } 1895 decpt = s += i; 1896#else 1897 if (*s != '.') 1898 goto pcheck; 1899 decpt = ++s; 1900#endif 1901 if (!hexdig[*s]) 1902 goto pcheck; 1903 while(*s == '0') 1904 s++; 1905 if (hexdig[*s]) 1906 zret = 0; 1907 havedig = 1; 1908 s0 = s; 1909 } 1910 while(hexdig[*s]) 1911 s++; 1912#ifdef USE_LOCALE 1913 if (*s == *decimalpoint && !decpt) { 1914 for(i = 1; decimalpoint[i]; ++i) { 1915 if (s[i] != decimalpoint[i]) 1916 goto pcheck; 1917 } 1918 decpt = s += i; 1919#else 1920 if (*s == '.' && !decpt) { 1921 decpt = ++s; 1922#endif 1923 while(hexdig[*s]) 1924 s++; 1925 }/*}*/ 1926 if (decpt) 1927 e = -(((Long)(s-decpt)) << 2); 1928 pcheck: 1929 s1 = s; 1930 big = esign = 0; 1931 switch(*s) { 1932 case 'p': 1933 case 'P': 1934 switch(*++s) { 1935 case '-': 1936 esign = 1; 1937 /* no break */ 1938 case '+': 1939 s++; 1940 } 1941 if ((n = hexdig[*s]) == 0 || n > 0x19) { 1942 s = s1; 1943 break; 1944 } 1945 e1 = n - 0x10; 1946 while((n = hexdig[*++s]) !=0 && n <= 0x19) { 1947 if (e1 & 0xf8000000) 1948 big = 1; 1949 e1 = 10*e1 + n - 0x10; 1950 } 1951 if (esign) 1952 e1 = -e1; 1953 e += e1; 1954 } 1955 *sp = (char*)s; 1956 if (!havedig) 1957 *sp = (char*)s0 - 1; 1958 if (zret) 1959 goto retz1; 1960 if (big) { 1961 if (esign) { 1962#ifdef IEEE_Arith 1963 switch(rounding) { 1964 case Round_up: 1965 if (sign) 1966 break; 1967 goto ret_tiny; 1968 case Round_down: 1969 if (!sign) 1970 break; 1971 goto ret_tiny; 1972 } 1973#endif 1974 goto retz; 1975#ifdef IEEE_Arith 1976 ret_tinyf: 1977 Bfree(b); 1978 ret_tiny: 1979#ifndef NO_ERRNO 1980 errno = ERANGE; 1981#endif 1982 word0(rvp) = 0; 1983 word1(rvp) = 1; 1984 return; 1985#endif /* IEEE_Arith */ 1986 } 1987 switch(rounding) { 1988 case Round_near: 1989 goto ovfl1; 1990 case Round_up: 1991 if (!sign) 1992 goto ovfl1; 1993 goto ret_big; 1994 case Round_down: 1995 if (sign) 1996 goto ovfl1; 1997 goto ret_big; 1998 } 1999 ret_big: 2000 word0(rvp) = Big0; 2001 word1(rvp) = Big1; 2002 return; 2003 } 2004 n = s1 - s0 - 1; 2005 for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1) 2006 k++; 2007 b = Balloc(k); 2008 x = b->x; 2009 n = 0; 2010 L = 0; 2011#ifdef USE_LOCALE 2012 for(i = 0; decimalpoint[i+1]; ++i); 2013#endif 2014 while(s1 > s0) { 2015#ifdef USE_LOCALE 2016 if (*--s1 == decimalpoint[i]) { 2017 s1 -= i; 2018 continue; 2019 } 2020#else 2021 if (*--s1 == '.') 2022 continue; 2023#endif 2024 if (n == ULbits) { 2025 *x++ = L; 2026 L = 0; 2027 n = 0; 2028 } 2029 L |= (hexdig[*s1] & 0x0f) << n; 2030 n += 4; 2031 } 2032 *x++ = L; 2033 b->wds = n = x - b->x; 2034 n = ULbits*n - hi0bits(L); 2035 nbits = Nbits; 2036 lostbits = 0; 2037 x = b->x; 2038 if (n > nbits) { 2039 n -= nbits; 2040 if (any_on(b,n)) { 2041 lostbits = 1; 2042 k = n - 1; 2043 if (x[k>>kshift] & 1 << (k & kmask)) { 2044 lostbits = 2; 2045 if (k > 0 && any_on(b,k)) 2046 lostbits = 3; 2047 } 2048 } 2049 rshift(b, n); 2050 e += n; 2051 } 2052 else if (n < nbits) { 2053 n = nbits - n; 2054 b = lshift(b, n); 2055 e -= n; 2056 x = b->x; 2057 } 2058 if (e > Emax) { 2059 ovfl: 2060 Bfree(b); 2061 ovfl1: 2062#ifndef NO_ERRNO 2063 errno = ERANGE; 2064#endif 2065 word0(rvp) = Exp_mask; 2066 word1(rvp) = 0; 2067 return; 2068 } 2069 denorm = 0; 2070 if (e < emin) { 2071 denorm = 1; 2072 n = emin - e; 2073 if (n >= nbits) { 2074#ifdef IEEE_Arith /*{*/ 2075 switch (rounding) { 2076 case Round_near: 2077 if (n == nbits && (n < 2 || any_on(b,n-1))) 2078 goto ret_tinyf; 2079 break; 2080 case Round_up: 2081 if (!sign) 2082 goto ret_tinyf; 2083 break; 2084 case Round_down: 2085 if (sign) 2086 goto ret_tinyf; 2087 } 2088#endif /* } IEEE_Arith */ 2089 Bfree(b); 2090 retz: 2091#ifndef NO_ERRNO 2092 errno = ERANGE; 2093#endif 2094 retz1: 2095 rvp->d = 0.; 2096 return; 2097 } 2098 k = n - 1; 2099 if (lostbits) 2100 lostbits = 1; 2101 else if (k > 0) 2102 lostbits = any_on(b,k); 2103 if (x[k>>kshift] & 1 << (k & kmask)) 2104 lostbits |= 2; 2105 nbits -= n; 2106 rshift(b,n); 2107 e = emin; 2108 } 2109 if (lostbits) { 2110 up = 0; 2111 switch(rounding) { 2112 case Round_zero: 2113 break; 2114 case Round_near: 2115 if (lostbits & 2 2116 && (lostbits & 1) | (x[0] & 1)) 2117 up = 1; 2118 break; 2119 case Round_up: 2120 up = 1 - sign; 2121 break; 2122 case Round_down: 2123 up = sign; 2124 } 2125 if (up) { 2126 k = b->wds; 2127 b = increment(b); 2128 x = b->x; 2129 if (denorm) { 2130#if 0 2131 if (nbits == Nbits - 1 2132 && x[nbits >> kshift] & 1 << (nbits & kmask)) 2133 denorm = 0; /* not currently used */ 2134#endif 2135 } 2136 else if (b->wds > k 2137 || ((n = nbits & kmask) !=0 2138 && hi0bits(x[k-1]) < 32-n)) { 2139 rshift(b,1); 2140 if (++e > Emax) 2141 goto ovfl; 2142 } 2143 } 2144 } 2145#ifdef IEEE_Arith 2146 if (denorm) 2147 word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0; 2148 else 2149 word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20); 2150 word1(rvp) = b->x[0]; 2151#endif 2152#ifdef IBM 2153 if ((j = e & 3)) { 2154 k = b->x[0] & ((1 << j) - 1); 2155 rshift(b,j); 2156 if (k) { 2157 switch(rounding) { 2158 case Round_up: 2159 if (!sign) 2160 increment(b); 2161 break; 2162 case Round_down: 2163 if (sign) 2164 increment(b); 2165 break; 2166 case Round_near: 2167 j = 1 << (j-1); 2168 if (k & j && ((k & (j-1)) | lostbits)) 2169 increment(b); 2170 } 2171 } 2172 } 2173 e >>= 2; 2174 word0(rvp) = b->x[1] | ((e + 65 + 13) << 24); 2175 word1(rvp) = b->x[0]; 2176#endif 2177#ifdef VAX 2178 /* The next two lines ignore swap of low- and high-order 2 bytes. */ 2179 /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */ 2180 /* word1(rvp) = b->x[0]; */ 2181 word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16); 2182 word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16); 2183#endif 2184 Bfree(b); 2185 } 2186#endif /*!NO_HEX_FP}*/ 2187 2188 static int 2189#ifdef KR_headers 2190dshift(b, p2) Bigint *b; int p2; 2191#else 2192dshift(Bigint *b, int p2) 2193#endif 2194{ 2195 int rv = hi0bits(b->x[b->wds-1]) - 4; 2196 if (p2 > 0) 2197 rv -= p2; 2198 return rv & kmask; 2199 } 2200 2201 static int 2202quorem 2203#ifdef KR_headers 2204 (b, S) Bigint *b, *S; 2205#else 2206 (Bigint *b, Bigint *S) 2207#endif 2208{ 2209 int n; 2210 ULong *bx, *bxe, q, *sx, *sxe; 2211#ifdef ULLong 2212 ULLong borrow, carry, y, ys; 2213#else 2214 ULong borrow, carry, y, ys; 2215#ifdef Pack_32 2216 ULong si, z, zs; 2217#endif 2218#endif 2219 2220 n = S->wds; 2221#ifdef DEBUG 2222 /*debug*/ if (b->wds > n) 2223 /*debug*/ Bug("oversize b in quorem"); 2224#endif 2225 if (b->wds < n) 2226 return 0; 2227 sx = S->x; 2228 sxe = sx + --n; 2229 bx = b->x; 2230 bxe = bx + n; 2231 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 2232#ifdef DEBUG 2233#ifdef NO_STRTOD_BIGCOMP 2234 /*debug*/ if (q > 9) 2235#else 2236 /* An oversized q is possible when quorem is called from bigcomp and */ 2237 /* the input is near, e.g., twice the smallest denormalized number. */ 2238 /*debug*/ if (q > 15) 2239#endif 2240 /*debug*/ Bug("oversized quotient in quorem"); 2241#endif 2242 if (q) { 2243 borrow = 0; 2244 carry = 0; 2245 do { 2246#ifdef ULLong 2247 ys = *sx++ * (ULLong)q + carry; 2248 carry = ys >> 32; 2249 y = *bx - (ys & FFFFFFFF) - borrow; 2250 borrow = y >> 32 & (ULong)1; 2251 *bx++ = y & FFFFFFFF; 2252#else 2253#ifdef Pack_32 2254 si = *sx++; 2255 ys = (si & 0xffff) * q + carry; 2256 zs = (si >> 16) * q + (ys >> 16); 2257 carry = zs >> 16; 2258 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 2259 borrow = (y & 0x10000) >> 16; 2260 z = (*bx >> 16) - (zs & 0xffff) - borrow; 2261 borrow = (z & 0x10000) >> 16; 2262 Storeinc(bx, z, y); 2263#else 2264 ys = *sx++ * q + carry; 2265 carry = ys >> 16; 2266 y = *bx - (ys & 0xffff) - borrow; 2267 borrow = (y & 0x10000) >> 16; 2268 *bx++ = y & 0xffff; 2269#endif 2270#endif 2271 } 2272 while(sx <= sxe); 2273 if (!*bxe) { 2274 bx = b->x; 2275 while(--bxe > bx && !*bxe) 2276 --n; 2277 b->wds = n; 2278 } 2279 } 2280 if (cmp(b, S) >= 0) { 2281 q++; 2282 borrow = 0; 2283 carry = 0; 2284 bx = b->x; 2285 sx = S->x; 2286 do { 2287#ifdef ULLong 2288 ys = *sx++ + carry; 2289 carry = ys >> 32; 2290 y = *bx - (ys & FFFFFFFF) - borrow; 2291 borrow = y >> 32 & (ULong)1; 2292 *bx++ = y & FFFFFFFF; 2293#else 2294#ifdef Pack_32 2295 si = *sx++; 2296 ys = (si & 0xffff) + carry; 2297 zs = (si >> 16) + (ys >> 16); 2298 carry = zs >> 16; 2299 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 2300 borrow = (y & 0x10000) >> 16; 2301 z = (*bx >> 16) - (zs & 0xffff) - borrow; 2302 borrow = (z & 0x10000) >> 16; 2303 Storeinc(bx, z, y); 2304#else 2305 ys = *sx++ + carry; 2306 carry = ys >> 16; 2307 y = *bx - (ys & 0xffff) - borrow; 2308 borrow = (y & 0x10000) >> 16; 2309 *bx++ = y & 0xffff; 2310#endif 2311#endif 2312 } 2313 while(sx <= sxe); 2314 bx = b->x; 2315 bxe = bx + n; 2316 if (!*bxe) { 2317 while(--bxe > bx && !*bxe) 2318 --n; 2319 b->wds = n; 2320 } 2321 } 2322 return q; 2323 } 2324 2325#ifndef HAVE_STRTOD 2326 2327#if defined(Avoid_Underflow) || !defined(NO_STRTOD_BIGCOMP) /*{*/ 2328 static double 2329sulp 2330#ifdef KR_headers 2331 (x, bc) U *x; BCinfo *bc; 2332#else 2333 (U *x, BCinfo *bc) 2334#endif 2335{ 2336 U u; 2337 double rv; 2338 int i; 2339 2340 rv = ulp(x); 2341 if (!bc->scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0) 2342 return rv; /* Is there an example where i <= 0 ? */ 2343 word0(&u) = Exp_1 + (i << Exp_shift); 2344 word1(&u) = 0; 2345 return rv * u.d; 2346 } 2347#endif /*}*/ 2348 2349#ifndef NO_STRTOD_BIGCOMP 2350 static void 2351bigcomp 2352#ifdef KR_headers 2353 (rv, s0, bc) 2354 U *rv; CONST char *s0; BCinfo *bc; 2355#else 2356 (U *rv, const char *s0, BCinfo *bc) 2357#endif 2358{ 2359 Bigint *b, *d; 2360 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase; 2361 2362 dsign = bc->dsign; 2363 nd = bc->nd; 2364 nd0 = bc->nd0; 2365 p5 = nd + bc->e0 - 1; 2366 speccase = 0; 2367#ifndef Sudden_Underflow 2368 if (rv->d == 0.) { /* special case: value near underflow-to-zero */ 2369 /* threshold was rounded to zero */ 2370 b = i2b(1); 2371 p2 = Emin - P + 1; 2372 bbits = 1; 2373#ifdef Avoid_Underflow 2374 word0(rv) = (P+2) << Exp_shift; 2375#else 2376 word1(rv) = 1; 2377#endif 2378 i = 0; 2379#ifdef Honor_FLT_ROUNDS 2380 if (bc->rounding == 1) 2381#endif 2382 { 2383 speccase = 1; 2384 --p2; 2385 dsign = 0; 2386 goto have_i; 2387 } 2388 } 2389 else 2390#endif 2391 b = d2b(rv, &p2, &bbits); 2392#ifdef Avoid_Underflow 2393 p2 -= bc->scale; 2394#endif 2395 /* floor(log2(rv)) == bbits - 1 + p2 */ 2396 /* Check for denormal case. */ 2397 i = P - bbits; 2398 if (i > (j = P - Emin - 1 + p2)) { 2399#ifdef Sudden_Underflow 2400 Bfree(b); 2401 b = i2b(1); 2402 p2 = Emin; 2403 i = P - 1; 2404#ifdef Avoid_Underflow 2405 word0(rv) = (1 + bc->scale) << Exp_shift; 2406#else 2407 word0(rv) = Exp_msk1; 2408#endif 2409 word1(rv) = 0; 2410#else 2411 i = j; 2412#endif 2413 } 2414#ifdef Honor_FLT_ROUNDS 2415 if (bc->rounding != 1) { 2416 if (i > 0) 2417 b = lshift(b, i); 2418 if (dsign) 2419 b = increment(b); 2420 } 2421 else 2422#endif 2423 { 2424 b = lshift(b, ++i); 2425 b->x[0] |= 1; 2426 } 2427#ifndef Sudden_Underflow 2428 have_i: 2429#endif 2430 p2 -= p5 + i; 2431 d = i2b(1); 2432 /* Arrange for convenient computation of quotients: 2433 * shift left if necessary so divisor has 4 leading 0 bits. 2434 */ 2435 if (p5 > 0) 2436 d = pow5mult(d, p5); 2437 else if (p5 < 0) 2438 b = pow5mult(b, -p5); 2439 if (p2 > 0) { 2440 b2 = p2; 2441 d2 = 0; 2442 } 2443 else { 2444 b2 = 0; 2445 d2 = -p2; 2446 } 2447 i = dshift(d, d2); 2448 if ((b2 += i) > 0) 2449 b = lshift(b, b2); 2450 if ((d2 += i) > 0) 2451 d = lshift(d, d2); 2452 2453 /* Now b/d = exactly half-way between the two floating-point values */ 2454 /* on either side of the input string. Compute first digit of b/d. */ 2455 2456 if (!(dig = quorem(b,d))) { 2457 b = multadd(b, 10, 0); /* very unlikely */ 2458 dig = quorem(b,d); 2459 } 2460 2461 /* Compare b/d with s0 */ 2462 2463 for(i = 0; i < nd0; ) { 2464 if ((dd = s0[i++] - '0' - dig)) 2465 goto ret; 2466 if (!b->x[0] && b->wds == 1) { 2467 if (i < nd) 2468 dd = 1; 2469 goto ret; 2470 } 2471 b = multadd(b, 10, 0); 2472 dig = quorem(b,d); 2473 } 2474 for(j = bc->dp1; i++ < nd;) { 2475 if ((dd = s0[j++] - '0' - dig)) 2476 goto ret; 2477 if (!b->x[0] && b->wds == 1) { 2478 if (i < nd) 2479 dd = 1; 2480 goto ret; 2481 } 2482 b = multadd(b, 10, 0); 2483 dig = quorem(b,d); 2484 } 2485 if (dig > 0 || b->x[0] || b->wds > 1) 2486 dd = -1; 2487 ret: 2488 Bfree(b); 2489 Bfree(d); 2490#ifdef Honor_FLT_ROUNDS 2491 if (bc->rounding != 1) { 2492 if (dd < 0) { 2493 if (bc->rounding == 0) { 2494 if (!dsign) 2495 goto retlow1; 2496 } 2497 else if (dsign) 2498 goto rethi1; 2499 } 2500 else if (dd > 0) { 2501 if (bc->rounding == 0) { 2502 if (dsign) 2503 goto rethi1; 2504 goto ret1; 2505 } 2506 if (!dsign) 2507 goto rethi1; 2508 dval(rv) += 2.*sulp(rv,bc); 2509 } 2510 else { 2511 bc->inexact = 0; 2512 if (dsign) 2513 goto rethi1; 2514 } 2515 } 2516 else 2517#endif 2518 if (speccase) { 2519 if (dd <= 0) 2520 rv->d = 0.; 2521 } 2522 else if (dd < 0) { 2523 if (!dsign) /* does not happen for round-near */ 2524retlow1: 2525 dval(rv) -= sulp(rv,bc); 2526 } 2527 else if (dd > 0) { 2528 if (dsign) { 2529 rethi1: 2530 dval(rv) += sulp(rv,bc); 2531 } 2532 } 2533 else { 2534 /* Exact half-way case: apply round-even rule. */ 2535 if ((j = ((word0(rv) & Exp_mask) >> Exp_shift) - bc->scale) <= 0) { 2536 i = 1 - j; 2537 if (i <= 31) { 2538 if (word1(rv) & (0x1 << i)) 2539 goto odd; 2540 } 2541 else if (word0(rv) & (0x1 << (i-32))) 2542 goto odd; 2543 } 2544 else if (word1(rv) & 1) { 2545 odd: 2546 if (dsign) 2547 goto rethi1; 2548 goto retlow1; 2549 } 2550 } 2551 2552#ifdef Honor_FLT_ROUNDS 2553 ret1: 2554#endif 2555 return; 2556 } 2557#endif /* NO_STRTOD_BIGCOMP */ 2558 2559 double 2560poly_strtod 2561#ifdef KR_headers 2562 (s00, se) CONST char *s00; char **se; 2563#else 2564 (const char *s00, char **se) 2565#endif 2566{ 2567 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1; 2568 int esign, i, j, k, nd, nd0, nf, nz, nz0, nz1, sign; 2569 CONST char *s, *s0, *s1; 2570 double aadj, aadj1; 2571 Long L; 2572 U aadj2, adj, rv, rv0; 2573 ULong y, z; 2574 BCinfo bc; 2575 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 2576#ifdef Avoid_Underflow 2577 ULong Lsb, Lsb1; 2578#endif 2579#ifdef SET_INEXACT 2580 int oldinexact; 2581#endif 2582#ifndef NO_STRTOD_BIGCOMP 2583 int req_bigcomp = 0; 2584#endif 2585#ifdef Honor_FLT_ROUNDS /*{*/ 2586#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ 2587 bc.rounding = Flt_Rounds; 2588#else /*}{*/ 2589 bc.rounding = 1; 2590 switch(fegetround()) { 2591 case FE_TOWARDZERO: bc.rounding = 0; break; 2592 case FE_UPWARD: bc.rounding = 2; break; 2593 case FE_DOWNWARD: bc.rounding = 3; 2594 } 2595#endif /*}}*/ 2596#endif /*}*/ 2597#ifdef USE_LOCALE 2598 CONST char *s2; 2599#endif 2600 2601 sign = nz0 = nz1 = nz = bc.dplen = bc.uflchk = 0; 2602 dval(&rv) = 0.; 2603 for(s = s00;;s++) switch(*s) { 2604 case '-': 2605 sign = 1; 2606 /* no break */ 2607 case '+': 2608 if (*++s) 2609 goto break2; 2610 /* no break */ 2611 case 0: 2612 goto ret0; 2613 case '\t': 2614 case '\n': 2615 case '\v': 2616 case '\f': 2617 case '\r': 2618 case ' ': 2619 continue; 2620 default: 2621 goto break2; 2622 } 2623 break2: 2624 if (*s == '0') { 2625#ifndef NO_HEX_FP /*{*/ 2626 switch(s[1]) { 2627 case 'x': 2628 case 'X': 2629#ifdef Honor_FLT_ROUNDS 2630 gethex(&s, &rv, bc.rounding, sign); 2631#else 2632 gethex(&s, &rv, 1, sign); 2633#endif 2634 goto ret; 2635 } 2636#endif /*}*/ 2637 nz0 = 1; 2638 while(*++s == '0') ; 2639 if (!*s) 2640 goto ret; 2641 } 2642 s0 = s; 2643 y = z = 0; 2644 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 2645 if (nd < 9) 2646 y = 10*y + c - '0'; 2647 else if (nd < DBL_DIG + 2) 2648 z = 10*z + c - '0'; 2649 nd0 = nd; 2650 bc.dp0 = bc.dp1 = s - s0; 2651 for(s1 = s; s1 > s0 && *--s1 == '0'; ) 2652 ++nz1; 2653#ifdef USE_LOCALE 2654 s1 = localeconv()->decimal_point; 2655 if (c == *s1) { 2656 c = '.'; 2657 if (*++s1) { 2658 s2 = s; 2659 for(;;) { 2660 if (*++s2 != *s1) { 2661 c = 0; 2662 break; 2663 } 2664 if (!*++s1) { 2665 s = s2; 2666 break; 2667 } 2668 } 2669 } 2670 } 2671#endif 2672 if (c == '.') { 2673 c = *++s; 2674 bc.dp1 = s - s0; 2675 bc.dplen = bc.dp1 - bc.dp0; 2676 if (!nd) { 2677 for(; c == '0'; c = *++s) 2678 nz++; 2679 if (c > '0' && c <= '9') { 2680 bc.dp0 = s0 - s; 2681 bc.dp1 = bc.dp0 + bc.dplen; 2682 s0 = s; 2683 nf += nz; 2684 nz = 0; 2685 goto have_dig; 2686 } 2687 goto dig_done; 2688 } 2689 for(; c >= '0' && c <= '9'; c = *++s) { 2690 have_dig: 2691 nz++; 2692 if (c -= '0') { 2693 nf += nz; 2694 for(i = 1; i < nz; i++) 2695 if (nd++ < 9) 2696 y *= 10; 2697 else if (nd <= DBL_DIG + 2) 2698 z *= 10; 2699 if (nd++ < 9) 2700 y = 10*y + c; 2701 else if (nd <= DBL_DIG + 2) 2702 z = 10*z + c; 2703 nz = nz1 = 0; 2704 } 2705 } 2706 } 2707 dig_done: 2708 e = 0; 2709 if (c == 'e' || c == 'E') { 2710 if (!nd && !nz && !nz0) { 2711 goto ret0; 2712 } 2713 s00 = s; 2714 esign = 0; 2715 switch(c = *++s) { 2716 case '-': 2717 esign = 1; 2718 case '+': 2719 c = *++s; 2720 } 2721 if (c >= '0' && c <= '9') { 2722 while(c == '0') 2723 c = *++s; 2724 if (c > '0' && c <= '9') { 2725 L = c - '0'; 2726 s1 = s; 2727 while((c = *++s) >= '0' && c <= '9') 2728 L = 10*L + c - '0'; 2729 if (s - s1 > 8 || L > 19999) 2730 /* Avoid confusion from exponents 2731 * so large that e might overflow. 2732 */ 2733 e = 19999; /* safe for 16 bit ints */ 2734 else 2735 e = (int)L; 2736 if (esign) 2737 e = -e; 2738 } 2739 else 2740 e = 0; 2741 } 2742 else 2743 s = s00; 2744 } 2745 if (!nd) { 2746 if (!nz && !nz0) { 2747#ifdef INFNAN_CHECK 2748 /* Check for Nan and Infinity */ 2749 if (!bc.dplen) 2750 switch(c) { 2751 case 'i': 2752 case 'I': 2753 if (match(&s,"nf")) { 2754 --s; 2755 if (!match(&s,"inity")) 2756 ++s; 2757 word0(&rv) = 0x7ff00000; 2758 word1(&rv) = 0; 2759 goto ret; 2760 } 2761 break; 2762 case 'n': 2763 case 'N': 2764 if (match(&s, "an")) { 2765 word0(&rv) = NAN_WORD0; 2766 word1(&rv) = NAN_WORD1; 2767#ifndef No_Hex_NaN 2768 if (*s == '(') /*)*/ 2769 hexnan(&rv, &s); 2770#endif 2771 goto ret; 2772 } 2773 } 2774#endif /* INFNAN_CHECK */ 2775 ret0: 2776 s = s00; 2777 sign = 0; 2778 } 2779 goto ret; 2780 } 2781 bc.e0 = e1 = e -= nf; 2782 2783 /* Now we have nd0 digits, starting at s0, followed by a 2784 * decimal point, followed by nd-nd0 digits. The number we're 2785 * after is the integer represented by those digits times 2786 * 10**e */ 2787 2788 if (!nd0) 2789 nd0 = nd; 2790 k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2; 2791 dval(&rv) = y; 2792 if (k > 9) { 2793#ifdef SET_INEXACT 2794 if (k > DBL_DIG) 2795 oldinexact = get_inexact(); 2796#endif 2797 dval(&rv) = tens[k - 9] * dval(&rv) + z; 2798 } 2799 bd0 = 0; 2800 if (nd <= DBL_DIG 2801#ifndef RND_PRODQUOT 2802#ifndef Honor_FLT_ROUNDS 2803 && Flt_Rounds == 1 2804#endif 2805#endif 2806 ) { 2807 if (!e) 2808 goto ret; 2809#ifndef ROUND_BIASED_without_Round_Up 2810 if (e > 0) { 2811 if (e <= Ten_pmax) { 2812#ifdef VAX 2813 goto vax_ovfl_check; 2814#else 2815#ifdef Honor_FLT_ROUNDS 2816 /* round correctly FLT_ROUNDS = 2 or 3 */ 2817 if (sign) { 2818 rv.d = -rv.d; 2819 sign = 0; 2820 } 2821#endif 2822 /* rv = */ rounded_product(dval(&rv), tens[e]); 2823 goto ret; 2824#endif 2825 } 2826 i = DBL_DIG - nd; 2827 if (e <= Ten_pmax + i) { 2828 /* A fancier test would sometimes let us do 2829 * this for larger i values. 2830 */ 2831#ifdef Honor_FLT_ROUNDS 2832 /* round correctly FLT_ROUNDS = 2 or 3 */ 2833 if (sign) { 2834 rv.d = -rv.d; 2835 sign = 0; 2836 } 2837#endif 2838 e -= i; 2839 dval(&rv) *= tens[i]; 2840#ifdef VAX 2841 /* VAX exponent range is so narrow we must 2842 * worry about overflow here... 2843 */ 2844 vax_ovfl_check: 2845 word0(&rv) -= P*Exp_msk1; 2846 /* rv = */ rounded_product(dval(&rv), tens[e]); 2847 if ((word0(&rv) & Exp_mask) 2848 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 2849 goto ovfl; 2850 word0(&rv) += P*Exp_msk1; 2851#else 2852 /* rv = */ rounded_product(dval(&rv), tens[e]); 2853#endif 2854 goto ret; 2855 } 2856 } 2857#ifndef Inaccurate_Divide 2858 else if (e >= -Ten_pmax) { 2859#ifdef Honor_FLT_ROUNDS 2860 /* round correctly FLT_ROUNDS = 2 or 3 */ 2861 if (sign) { 2862 rv.d = -rv.d; 2863 sign = 0; 2864 } 2865#endif 2866 /* rv = */ rounded_quotient(dval(&rv), tens[-e]); 2867 goto ret; 2868 } 2869#endif 2870#endif /* ROUND_BIASED_without_Round_Up */ 2871 } 2872 e1 += nd - k; 2873 2874#ifdef IEEE_Arith 2875#ifdef SET_INEXACT 2876 bc.inexact = 1; 2877 if (k <= DBL_DIG) 2878 oldinexact = get_inexact(); 2879#endif 2880#ifdef Avoid_Underflow 2881 bc.scale = 0; 2882#endif 2883#ifdef Honor_FLT_ROUNDS 2884 if (bc.rounding >= 2) { 2885 if (sign) 2886 bc.rounding = bc.rounding == 2 ? 0 : 2; 2887 else 2888 if (bc.rounding != 2) 2889 bc.rounding = 0; 2890 } 2891#endif 2892#endif /*IEEE_Arith*/ 2893 2894 /* Get starting approximation = rv * 10**e1 */ 2895 2896 if (e1 > 0) { 2897 if ((i = e1 & 15)) 2898 dval(&rv) *= tens[i]; 2899 if (e1 &= ~15) { 2900 if (e1 > DBL_MAX_10_EXP) { 2901 ovfl: 2902 /* Can't trust HUGE_VAL */ 2903#ifdef IEEE_Arith 2904#ifdef Honor_FLT_ROUNDS 2905 switch(bc.rounding) { 2906 case 0: /* toward 0 */ 2907 case 3: /* toward -infinity */ 2908 word0(&rv) = Big0; 2909 word1(&rv) = Big1; 2910 break; 2911 default: 2912 word0(&rv) = Exp_mask; 2913 word1(&rv) = 0; 2914 } 2915#else /*Honor_FLT_ROUNDS*/ 2916 word0(&rv) = Exp_mask; 2917 word1(&rv) = 0; 2918#endif /*Honor_FLT_ROUNDS*/ 2919#ifdef SET_INEXACT 2920 /* set overflow bit */ 2921 dval(&rv0) = 1e300; 2922 dval(&rv0) *= dval(&rv0); 2923#endif 2924#else /*IEEE_Arith*/ 2925 word0(&rv) = Big0; 2926 word1(&rv) = Big1; 2927#endif /*IEEE_Arith*/ 2928 range_err: 2929 if (bd0) { 2930 Bfree(bb); 2931 Bfree(bd); 2932 Bfree(bs); 2933 Bfree(bd0); 2934 Bfree(delta); 2935 } 2936#ifndef NO_ERRNO 2937 errno = ERANGE; 2938#endif 2939 goto ret; 2940 } 2941 e1 >>= 4; 2942 for(j = 0; e1 > 1; j++, e1 >>= 1) 2943 if (e1 & 1) 2944 dval(&rv) *= bigtens[j]; 2945 /* The last multiplication could overflow. */ 2946 word0(&rv) -= P*Exp_msk1; 2947 dval(&rv) *= bigtens[j]; 2948 if ((z = word0(&rv) & Exp_mask) 2949 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 2950 goto ovfl; 2951 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 2952 /* set to largest number */ 2953 /* (Can't trust DBL_MAX) */ 2954 word0(&rv) = Big0; 2955 word1(&rv) = Big1; 2956 } 2957 else 2958 word0(&rv) += P*Exp_msk1; 2959 } 2960 } 2961 else if (e1 < 0) { 2962 e1 = -e1; 2963 if ((i = e1 & 15)) 2964 dval(&rv) /= tens[i]; 2965 if (e1 >>= 4) { 2966 if (e1 >= 1 << n_bigtens) 2967 goto undfl; 2968#ifdef Avoid_Underflow 2969 if (e1 & Scale_Bit) 2970 bc.scale = 2*P; 2971 for(j = 0; e1 > 0; j++, e1 >>= 1) 2972 if (e1 & 1) 2973 dval(&rv) *= tinytens[j]; 2974 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) 2975 >> Exp_shift)) > 0) { 2976 /* scaled rv is denormal; clear j low bits */ 2977 if (j >= 32) { 2978 if (j > 54) 2979 goto undfl; 2980 word1(&rv) = 0; 2981 if (j >= 53) 2982 word0(&rv) = (P+2)*Exp_msk1; 2983 else 2984 word0(&rv) &= 0xffffffff << (j-32); 2985 } 2986 else 2987 word1(&rv) &= 0xffffffff << j; 2988 } 2989#else 2990 for(j = 0; e1 > 1; j++, e1 >>= 1) 2991 if (e1 & 1) 2992 dval(&rv) *= tinytens[j]; 2993 /* The last multiplication could underflow. */ 2994 dval(&rv0) = dval(&rv); 2995 dval(&rv) *= tinytens[j]; 2996 if (!dval(&rv)) { 2997 dval(&rv) = 2.*dval(&rv0); 2998 dval(&rv) *= tinytens[j]; 2999#endif 3000 if (!dval(&rv)) { 3001 undfl: 3002 dval(&rv) = 0.; 3003 goto range_err; 3004 } 3005#ifndef Avoid_Underflow 3006 word0(&rv) = Tiny0; 3007 word1(&rv) = Tiny1; 3008 /* The refinement below will clean 3009 * this approximation up. 3010 */ 3011 } 3012#endif 3013 } 3014 } 3015 3016 /* Now the hard part -- adjusting rv to the correct value.*/ 3017 3018 /* Put digits into bd: true value = bd * 10^e */ 3019 3020 bc.nd = nd - nz1; 3021#ifndef NO_STRTOD_BIGCOMP 3022 bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */ 3023 /* to silence an erroneous warning about bc.nd0 */ 3024 /* possibly not being initialized. */ 3025 if (nd > strtod_diglim) { 3026 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */ 3027 /* minimum number of decimal digits to distinguish double values */ 3028 /* in IEEE arithmetic. */ 3029 i = j = 18; 3030 if (i > nd0) 3031 j += bc.dplen; 3032 for(;;) { 3033 if (--j < bc.dp1 && j >= bc.dp0) 3034 j = bc.dp0 - 1; 3035 if (s0[j] != '0') 3036 break; 3037 --i; 3038 } 3039 e += nd - i; 3040 nd = i; 3041 if (nd0 > nd) 3042 nd0 = nd; 3043 if (nd < 9) { /* must recompute y */ 3044 y = 0; 3045 for(i = 0; i < nd0; ++i) 3046 y = 10*y + s0[i] - '0'; 3047 for(j = bc.dp1; i < nd; ++i) 3048 y = 10*y + s0[j++] - '0'; 3049 } 3050 } 3051#endif 3052 bd0 = s2b(s0, nd0, nd, y, bc.dplen); 3053 3054 for(;;) { 3055 bd = Balloc(bd0->k); 3056 Bcopy(bd, bd0); 3057 bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */ 3058 bs = i2b(1); 3059 3060 if (e >= 0) { 3061 bb2 = bb5 = 0; 3062 bd2 = bd5 = e; 3063 } 3064 else { 3065 bb2 = bb5 = -e; 3066 bd2 = bd5 = 0; 3067 } 3068 if (bbe >= 0) 3069 bb2 += bbe; 3070 else 3071 bd2 -= bbe; 3072 bs2 = bb2; 3073#ifdef Honor_FLT_ROUNDS 3074 if (bc.rounding != 1) 3075 bs2++; 3076#endif 3077#ifdef Avoid_Underflow 3078 Lsb = LSB; 3079 Lsb1 = 0; 3080 j = bbe - bc.scale; 3081 i = j + bbbits - 1; /* logb(rv) */ 3082 j = P + 1 - bbbits; 3083 if (i < Emin) { /* denormal */ 3084 i = Emin - i; 3085 j -= i; 3086 if (i < 32) 3087 Lsb <<= i; 3088 else if (i < 52) 3089 Lsb1 = Lsb << (i-32); 3090 else 3091 Lsb1 = Exp_mask; 3092 } 3093#else /*Avoid_Underflow*/ 3094#ifdef Sudden_Underflow 3095#ifdef IBM 3096 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 3097#else 3098 j = P + 1 - bbbits; 3099#endif 3100#else /*Sudden_Underflow*/ 3101 j = bbe; 3102 i = j + bbbits - 1; /* logb(rv) */ 3103 if (i < Emin) /* denormal */ 3104 j += P - Emin; 3105 else 3106 j = P + 1 - bbbits; 3107#endif /*Sudden_Underflow*/ 3108#endif /*Avoid_Underflow*/ 3109 bb2 += j; 3110 bd2 += j; 3111#ifdef Avoid_Underflow 3112 bd2 += bc.scale; 3113#endif 3114 i = bb2 < bd2 ? bb2 : bd2; 3115 if (i > bs2) 3116 i = bs2; 3117 if (i > 0) { 3118 bb2 -= i; 3119 bd2 -= i; 3120 bs2 -= i; 3121 } 3122 if (bb5 > 0) { 3123 bs = pow5mult(bs, bb5); 3124 bb1 = mult(bs, bb); 3125 Bfree(bb); 3126 bb = bb1; 3127 } 3128 if (bb2 > 0) 3129 bb = lshift(bb, bb2); 3130 if (bd5 > 0) 3131 bd = pow5mult(bd, bd5); 3132 if (bd2 > 0) 3133 bd = lshift(bd, bd2); 3134 if (bs2 > 0) 3135 bs = lshift(bs, bs2); 3136 delta = diff(bb, bd); 3137 bc.dsign = delta->sign; 3138 delta->sign = 0; 3139 i = cmp(delta, bs); 3140#ifndef NO_STRTOD_BIGCOMP /*{*/ 3141 if (bc.nd > nd && i <= 0) { 3142 if (bc.dsign) { 3143 /* Must use bigcomp(). */ 3144 req_bigcomp = 1; 3145 break; 3146 } 3147#ifdef Honor_FLT_ROUNDS 3148 if (bc.rounding != 1) { 3149 if (i < 0) { 3150 req_bigcomp = 1; 3151 break; 3152 } 3153 } 3154 else 3155#endif 3156 i = -1; /* Discarded digits make delta smaller. */ 3157 } 3158#endif /*}*/ 3159#ifdef Honor_FLT_ROUNDS /*{*/ 3160 if (bc.rounding != 1) { 3161 if (i < 0) { 3162 /* Error is less than an ulp */ 3163 if (!delta->x[0] && delta->wds <= 1) { 3164 /* exact */ 3165#ifdef SET_INEXACT 3166 bc.inexact = 0; 3167#endif 3168 break; 3169 } 3170 if (bc.rounding) { 3171 if (bc.dsign) { 3172 adj.d = 1.; 3173 goto apply_adj; 3174 } 3175 } 3176 else if (!bc.dsign) { 3177 adj.d = -1.; 3178 if (!word1(&rv) 3179 && !(word0(&rv) & Frac_mask)) { 3180 y = word0(&rv) & Exp_mask; 3181#ifdef Avoid_Underflow 3182 if (!bc.scale || y > 2*P*Exp_msk1) 3183#else 3184 if (y) 3185#endif 3186 { 3187 delta = lshift(delta,Log2P); 3188 if (cmp(delta, bs) <= 0) 3189 adj.d = -0.5; 3190 } 3191 } 3192 apply_adj: 3193#ifdef Avoid_Underflow /*{*/ 3194 if (bc.scale && (y = word0(&rv) & Exp_mask) 3195 <= 2*P*Exp_msk1) 3196 word0(&adj) += (2*P+1)*Exp_msk1 - y; 3197#else 3198#ifdef Sudden_Underflow 3199 if ((word0(&rv) & Exp_mask) <= 3200 P*Exp_msk1) { 3201 word0(&rv) += P*Exp_msk1; 3202 dval(&rv) += adj.d*ulp(dval(&rv)); 3203 word0(&rv) -= P*Exp_msk1; 3204 } 3205 else 3206#endif /*Sudden_Underflow*/ 3207#endif /*Avoid_Underflow}*/ 3208 dval(&rv) += adj.d*ulp(&rv); 3209 } 3210 break; 3211 } 3212 adj.d = ratio(delta, bs); 3213 if (adj.d < 1.) 3214 adj.d = 1.; 3215 if (adj.d <= 0x7ffffffe) { 3216 /* adj = rounding ? ceil(adj) : floor(adj); */ 3217 y = adj.d; 3218 if (y != adj.d) { 3219 if (!((bc.rounding>>1) ^ bc.dsign)) 3220 y++; 3221 adj.d = y; 3222 } 3223 } 3224#ifdef Avoid_Underflow /*{*/ 3225 if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) 3226 word0(&adj) += (2*P+1)*Exp_msk1 - y; 3227#else 3228#ifdef Sudden_Underflow 3229 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { 3230 word0(&rv) += P*Exp_msk1; 3231 adj.d *= ulp(dval(&rv)); 3232 if (bc.dsign) 3233 dval(&rv) += adj.d; 3234 else 3235 dval(&rv) -= adj.d; 3236 word0(&rv) -= P*Exp_msk1; 3237 goto cont; 3238 } 3239#endif /*Sudden_Underflow*/ 3240#endif /*Avoid_Underflow}*/ 3241 adj.d *= ulp(&rv); 3242 if (bc.dsign) { 3243 if (word0(&rv) == Big0 && word1(&rv) == Big1) 3244 goto ovfl; 3245 dval(&rv) += adj.d; 3246 } 3247 else 3248 dval(&rv) -= adj.d; 3249 goto cont; 3250 } 3251#endif /*}Honor_FLT_ROUNDS*/ 3252 3253 if (i < 0) { 3254 /* Error is less than half an ulp -- check for 3255 * special case of mantissa a power of two. 3256 */ 3257 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask 3258#ifdef IEEE_Arith /*{*/ 3259#ifdef Avoid_Underflow 3260 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 3261#else 3262 || (word0(&rv) & Exp_mask) <= Exp_msk1 3263#endif 3264#endif /*}*/ 3265 ) { 3266#ifdef SET_INEXACT 3267 if (!delta->x[0] && delta->wds <= 1) 3268 bc.inexact = 0; 3269#endif 3270 break; 3271 } 3272 if (!delta->x[0] && delta->wds <= 1) { 3273 /* exact result */ 3274#ifdef SET_INEXACT 3275 bc.inexact = 0; 3276#endif 3277 break; 3278 } 3279 delta = lshift(delta,Log2P); 3280 if (cmp(delta, bs) > 0) 3281 goto drop_down; 3282 break; 3283 } 3284 if (i == 0) { 3285 /* exactly half-way between */ 3286 if (bc.dsign) { 3287 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 3288 && word1(&rv) == ( 3289#ifdef Avoid_Underflow 3290 (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) 3291 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 3292#endif 3293 0xffffffff)) { 3294 /*boundary case -- increment exponent*/ 3295 if (word0(&rv) == Big0 && word1(&rv) == Big1) 3296 goto ovfl; 3297 word0(&rv) = (word0(&rv) & Exp_mask) 3298 + Exp_msk1 3299#ifdef IBM 3300 | Exp_msk1 >> 4 3301#endif 3302 ; 3303 word1(&rv) = 0; 3304#ifdef Avoid_Underflow 3305 bc.dsign = 0; 3306#endif 3307 break; 3308 } 3309 } 3310 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) { 3311 drop_down: 3312 /* boundary case -- decrement exponent */ 3313#ifdef Sudden_Underflow /*{{*/ 3314 L = word0(&rv) & Exp_mask; 3315#ifdef IBM 3316 if (L < Exp_msk1) 3317#else 3318#ifdef Avoid_Underflow 3319 if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 3320#else 3321 if (L <= Exp_msk1) 3322#endif /*Avoid_Underflow*/ 3323#endif /*IBM*/ 3324 { 3325 if (bc.nd >nd) { 3326 bc.uflchk = 1; 3327 break; 3328 } 3329 goto undfl; 3330 } 3331 L -= Exp_msk1; 3332#else /*Sudden_Underflow}{*/ 3333#ifdef Avoid_Underflow 3334 if (bc.scale) { 3335 L = word0(&rv) & Exp_mask; 3336 if (L <= (2*P+1)*Exp_msk1) { 3337 if (L > (P+2)*Exp_msk1) 3338 /* round even ==> */ 3339 /* accept rv */ 3340 break; 3341 /* rv = smallest denormal */ 3342 if (bc.nd >nd) { 3343 bc.uflchk = 1; 3344 break; 3345 } 3346 goto undfl; 3347 } 3348 } 3349#endif /*Avoid_Underflow*/ 3350 L = (word0(&rv) & Exp_mask) - Exp_msk1; 3351#endif /*Sudden_Underflow}}*/ 3352 word0(&rv) = L | Bndry_mask1; 3353 word1(&rv) = 0xffffffff; 3354#ifdef IBM 3355 goto cont; 3356#else 3357#ifndef NO_STRTOD_BIGCOMP 3358 if (bc.nd > nd) 3359 goto cont; 3360#endif 3361 break; 3362#endif 3363 } 3364#ifndef ROUND_BIASED 3365#ifdef Avoid_Underflow 3366 if (Lsb1) { 3367 if (!(word0(&rv) & Lsb1)) 3368 break; 3369 } 3370 else if (!(word1(&rv) & Lsb)) 3371 break; 3372#else 3373 if (!(word1(&rv) & LSB)) 3374 break; 3375#endif 3376#endif 3377 if (bc.dsign) 3378#ifdef Avoid_Underflow 3379 dval(&rv) += sulp(&rv, &bc); 3380#else 3381 dval(&rv) += ulp(&rv); 3382#endif 3383#ifndef ROUND_BIASED 3384 else { 3385#ifdef Avoid_Underflow 3386 dval(&rv) -= sulp(&rv, &bc); 3387#else 3388 dval(&rv) -= ulp(&rv); 3389#endif 3390#ifndef Sudden_Underflow 3391 if (!dval(&rv)) { 3392 if (bc.nd >nd) { 3393 bc.uflchk = 1; 3394 break; 3395 } 3396 goto undfl; 3397 } 3398#endif 3399 } 3400#ifdef Avoid_Underflow 3401 bc.dsign = 1 - bc.dsign; 3402#endif 3403#endif 3404 break; 3405 } 3406 if ((aadj = ratio(delta, bs)) <= 2.) { 3407 if (bc.dsign) 3408 aadj = aadj1 = 1.; 3409 else if (word1(&rv) || word0(&rv) & Bndry_mask) { 3410#ifndef Sudden_Underflow 3411 if (word1(&rv) == Tiny1 && !word0(&rv)) { 3412 if (bc.nd >nd) { 3413 bc.uflchk = 1; 3414 break; 3415 } 3416 goto undfl; 3417 } 3418#endif 3419 aadj = 1.; 3420 aadj1 = -1.; 3421 } 3422 else { 3423 /* special case -- power of FLT_RADIX to be */ 3424 /* rounded down... */ 3425 3426 if (aadj < 2./FLT_RADIX) 3427 aadj = 1./FLT_RADIX; 3428 else 3429 aadj *= 0.5; 3430 aadj1 = -aadj; 3431 } 3432 } 3433 else { 3434 aadj *= 0.5; 3435 aadj1 = bc.dsign ? aadj : -aadj; 3436#ifdef Check_FLT_ROUNDS 3437 switch(bc.rounding) { 3438 case 2: /* towards +infinity */ 3439 aadj1 -= 0.5; 3440 break; 3441 case 0: /* towards 0 */ 3442 case 3: /* towards -infinity */ 3443 aadj1 += 0.5; 3444 } 3445#else 3446 if (Flt_Rounds == 0) 3447 aadj1 += 0.5; 3448#endif /*Check_FLT_ROUNDS*/ 3449 } 3450 y = word0(&rv) & Exp_mask; 3451 3452 /* Check for overflow */ 3453 3454 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 3455 dval(&rv0) = dval(&rv); 3456 word0(&rv) -= P*Exp_msk1; 3457 adj.d = aadj1 * ulp(&rv); 3458 dval(&rv) += adj.d; 3459 if ((word0(&rv) & Exp_mask) >= 3460 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 3461 if (word0(&rv0) == Big0 && word1(&rv0) == Big1) 3462 goto ovfl; 3463 word0(&rv) = Big0; 3464 word1(&rv) = Big1; 3465 goto cont; 3466 } 3467 else 3468 word0(&rv) += P*Exp_msk1; 3469 } 3470 else { 3471#ifdef Avoid_Underflow 3472 if (bc.scale && y <= 2*P*Exp_msk1) { 3473 if (aadj <= 0x7fffffff) { 3474 if ((z = aadj) <= 0) 3475 z = 1; 3476 aadj = z; 3477 aadj1 = bc.dsign ? aadj : -aadj; 3478 } 3479 dval(&aadj2) = aadj1; 3480 word0(&aadj2) += (2*P+1)*Exp_msk1 - y; 3481 aadj1 = dval(&aadj2); 3482 adj.d = aadj1 * ulp(&rv); 3483 dval(&rv) += adj.d; 3484 if (rv.d == 0.) 3485#ifdef NO_STRTOD_BIGCOMP 3486 goto undfl; 3487#else 3488 { 3489 req_bigcomp = 1; 3490 break; 3491 } 3492#endif 3493 } 3494 else { 3495 adj.d = aadj1 * ulp(&rv); 3496 dval(&rv) += adj.d; 3497 } 3498#else 3499#ifdef Sudden_Underflow 3500 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { 3501 dval(&rv0) = dval(&rv); 3502 word0(&rv) += P*Exp_msk1; 3503 adj.d = aadj1 * ulp(&rv); 3504 dval(&rv) += adj.d; 3505#ifdef IBM 3506 if ((word0(&rv) & Exp_mask) < P*Exp_msk1) 3507#else 3508 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) 3509#endif 3510 { 3511 if (word0(&rv0) == Tiny0 3512 && word1(&rv0) == Tiny1) { 3513 if (bc.nd >nd) { 3514 bc.uflchk = 1; 3515 break; 3516 } 3517 goto undfl; 3518 } 3519 word0(&rv) = Tiny0; 3520 word1(&rv) = Tiny1; 3521 goto cont; 3522 } 3523 else 3524 word0(&rv) -= P*Exp_msk1; 3525 } 3526 else { 3527 adj.d = aadj1 * ulp(&rv); 3528 dval(&rv) += adj.d; 3529 } 3530#else /*Sudden_Underflow*/ 3531 /* Compute adj so that the IEEE rounding rules will 3532 * correctly round rv + adj in some half-way cases. 3533 * If rv * ulp(rv) is denormalized (i.e., 3534 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 3535 * trouble from bits lost to denormalization; 3536 * example: 1.2e-307 . 3537 */ 3538 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 3539 aadj1 = (double)(int)(aadj + 0.5); 3540 if (!bc.dsign) 3541 aadj1 = -aadj1; 3542 } 3543 adj.d = aadj1 * ulp(&rv); 3544 dval(&rv) += adj.d; 3545#endif /*Sudden_Underflow*/ 3546#endif /*Avoid_Underflow*/ 3547 } 3548 z = word0(&rv) & Exp_mask; 3549#ifndef SET_INEXACT 3550 if (bc.nd == nd) { 3551#ifdef Avoid_Underflow 3552 if (!bc.scale) 3553#endif 3554 if (y == z) { 3555 /* Can we stop now? */ 3556 L = (Long)aadj; 3557 aadj -= L; 3558 /* The tolerances below are conservative. */ 3559 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) { 3560 if (aadj < .4999999 || aadj > .5000001) 3561 break; 3562 } 3563 else if (aadj < .4999999/FLT_RADIX) 3564 break; 3565 } 3566 } 3567#endif 3568 cont: 3569 Bfree(bb); 3570 Bfree(bd); 3571 Bfree(bs); 3572 Bfree(delta); 3573 } 3574 Bfree(bb); 3575 Bfree(bd); 3576 Bfree(bs); 3577 Bfree(bd0); 3578 Bfree(delta); 3579#ifndef NO_STRTOD_BIGCOMP 3580 if (req_bigcomp) { 3581 bd0 = 0; 3582 bc.e0 += nz1; 3583 bigcomp(&rv, s0, &bc); 3584 y = word0(&rv) & Exp_mask; 3585 if (y == Exp_mask) 3586 goto ovfl; 3587 if (y == 0 && rv.d == 0.) 3588 goto undfl; 3589 } 3590#endif 3591#ifdef SET_INEXACT 3592 if (bc.inexact) { 3593 if (!oldinexact) { 3594 word0(&rv0) = Exp_1 + (70 << Exp_shift); 3595 word1(&rv0) = 0; 3596 dval(&rv0) += 1.; 3597 } 3598 } 3599 else if (!oldinexact) 3600 clear_inexact(); 3601#endif 3602#ifdef Avoid_Underflow 3603 if (bc.scale) { 3604 word0(&rv0) = Exp_1 - 2*P*Exp_msk1; 3605 word1(&rv0) = 0; 3606 dval(&rv) *= dval(&rv0); 3607#ifndef NO_ERRNO 3608 /* try to avoid the bug of testing an 8087 register value */ 3609#ifdef IEEE_Arith 3610 if (!(word0(&rv) & Exp_mask)) 3611#else 3612 if (word0(&rv) == 0 && word1(&rv) == 0) 3613#endif 3614 errno = ERANGE; 3615#endif 3616 } 3617#endif /* Avoid_Underflow */ 3618#ifdef SET_INEXACT 3619 if (bc.inexact && !(word0(&rv) & Exp_mask)) { 3620 /* set underflow bit */ 3621 dval(&rv0) = 1e-300; 3622 dval(&rv0) *= dval(&rv0); 3623 } 3624#endif 3625 ret: 3626 if (se) 3627 *se = (char *)s; 3628 return sign ? -dval(&rv) : dval(&rv); 3629 } 3630 3631#endif // HAVE_STRTOD 3632 3633#ifndef MULTIPLE_THREADS 3634 static char *dtoa_result; 3635#endif 3636 3637 static char * 3638#ifdef KR_headers 3639rv_alloc(i) int i; 3640#else 3641rv_alloc(int i) 3642#endif 3643{ 3644 int j, k, *r; 3645 3646 j = sizeof(ULong); 3647 for(k = 0; 3648 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i; 3649 j <<= 1) 3650 k++; 3651 r = (int*)Balloc(k); 3652 *r = k; 3653 return 3654#ifndef MULTIPLE_THREADS 3655 dtoa_result = 3656#endif 3657 (char *)(r+1); 3658 } 3659 3660 static char * 3661#ifdef KR_headers 3662nrv_alloc(s, rve, n) char *s, **rve; int n; 3663#else 3664nrv_alloc(const char *s, char **rve, int n) 3665#endif 3666{ 3667 char *rv, *t; 3668 3669 t = rv = rv_alloc(n); 3670 while((*t = *s++)) t++; 3671 if (rve) 3672 *rve = t; 3673 return rv; 3674 } 3675 3676/* freedtoa(s) must be used to free values s returned by dtoa 3677 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 3678 * but for consistency with earlier versions of dtoa, it is optional 3679 * when MULTIPLE_THREADS is not defined. 3680 */ 3681 3682 void 3683#ifdef KR_headers 3684poly_freedtoa(s) char *s; 3685#else 3686poly_freedtoa(char *s) 3687#endif 3688{ 3689 Bigint *b = (Bigint *)((int *)s - 1); 3690 b->maxwds = 1 << (b->k = *(int*)b); 3691 Bfree(b); 3692#ifndef MULTIPLE_THREADS 3693 if (s == dtoa_result) 3694 dtoa_result = 0; 3695#endif 3696 } 3697 3698/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 3699 * 3700 * Inspired by "How to Print Floating-Point Numbers Accurately" by 3701 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. 3702 * 3703 * Modifications: 3704 * 1. Rather than iterating, we use a simple numeric overestimate 3705 * to determine k = floor(log10(d)). We scale relevant 3706 * quantities using O(log2(k)) rather than O(k) multiplications. 3707 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 3708 * try to generate digits strictly left to right. Instead, we 3709 * compute with fewer bits and propagate the carry if necessary 3710 * when rounding the final digit up. This is often faster. 3711 * 3. Under the assumption that input will be rounded nearest, 3712 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 3713 * That is, we allow equality in stopping tests when the 3714 * round-nearest rule will give the same floating-point value 3715 * as would satisfaction of the stopping test with strict 3716 * inequality. 3717 * 4. We remove common factors of powers of 2 from relevant 3718 * quantities. 3719 * 5. When converting floating-point integers less than 1e16, 3720 * we use floating-point arithmetic rather than resorting 3721 * to multiple-precision integers. 3722 * 6. When asked to produce fewer than 15 digits, we first try 3723 * to get by with floating-point arithmetic; we resort to 3724 * multiple-precision integer arithmetic only if we cannot 3725 * guarantee that the floating-point calculation has given 3726 * the correctly rounded result. For k requested digits and 3727 * "uniformly" distributed input, the probability is 3728 * something like 10^(k-15) that we must resort to the Long 3729 * calculation. 3730 */ 3731 3732 char * 3733poly_dtoa 3734#ifdef KR_headers 3735 (dd, mode, ndigits, decpt, sign, rve) 3736 double dd; int mode, ndigits, *decpt, *sign; char **rve; 3737#else 3738 (double dd, int mode, int ndigits, int *decpt, int *sign, char **rve) 3739#endif 3740{ 3741 /* Arguments ndigits, decpt, sign are similar to those 3742 of ecvt and fcvt; trailing zeros are suppressed from 3743 the returned string. If not null, *rve is set to point 3744 to the end of the return value. If d is +-Infinity or NaN, 3745 then *decpt is set to 9999. 3746 3747 mode: 3748 0 ==> shortest string that yields d when read in 3749 and rounded to nearest. 3750 1 ==> like 0, but with Steele & White stopping rule; 3751 e.g. with IEEE P754 arithmetic , mode 0 gives 3752 1e23 whereas mode 1 gives 9.999999999999999e22. 3753 2 ==> max(1,ndigits) significant digits. This gives a 3754 return value similar to that of ecvt, except 3755 that trailing zeros are suppressed. 3756 3 ==> through ndigits past the decimal point. This 3757 gives a return value similar to that from fcvt, 3758 except that trailing zeros are suppressed, and 3759 ndigits can be negative. 3760 4,5 ==> similar to 2 and 3, respectively, but (in 3761 round-nearest mode) with the tests of mode 0 to 3762 possibly return a shorter string that rounds to d. 3763 With IEEE arithmetic and compilation with 3764 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same 3765 as modes 2 and 3 when FLT_ROUNDS != 1. 3766 6-9 ==> Debugging modes similar to mode - 4: don't try 3767 fast floating-point estimate (if applicable). 3768 3769 Values of mode other than 0-9 are treated as mode 0. 3770 3771 Sufficient space is allocated to the return value 3772 to hold the suppressed trailing zeros. 3773 */ 3774 3775 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, 3776 j, j1=0, k, k0, k_check, leftright, m2, m5, s2, s5, 3777 spec_case, try_quick; 3778 Long L; 3779#ifndef Sudden_Underflow 3780 int denorm; 3781 ULong x; 3782#endif 3783 Bigint *b, *b1, *delta, *mlo=0, *mhi, *S; 3784 U d2, eps, u; 3785 double ds; 3786 char *s, *s0; 3787#ifndef No_leftright 3788#ifdef IEEE_Arith 3789 U eps1; 3790#endif 3791#endif 3792#ifdef SET_INEXACT 3793 int inexact, oldinexact; 3794#endif 3795#ifdef Honor_FLT_ROUNDS /*{*/ 3796 int Rounding; 3797#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ 3798 Rounding = Flt_Rounds; 3799#else /*}{*/ 3800 Rounding = 1; 3801 switch(fegetround()) { 3802 case FE_TOWARDZERO: Rounding = 0; break; 3803 case FE_UPWARD: Rounding = 2; break; 3804 case FE_DOWNWARD: Rounding = 3; 3805 } 3806#endif /*}}*/ 3807#endif /*}*/ 3808 3809#ifndef MULTIPLE_THREADS 3810 if (dtoa_result) { 3811 freedtoa(dtoa_result); 3812 dtoa_result = 0; 3813 } 3814#endif 3815 3816 u.d = dd; 3817 if (word0(&u) & Sign_bit) { 3818 /* set sign for everything, including 0's and NaNs */ 3819 *sign = 1; 3820 word0(&u) &= ~Sign_bit; /* clear sign bit */ 3821 } 3822 else 3823 *sign = 0; 3824 3825#if defined(IEEE_Arith) + defined(VAX) 3826#ifdef IEEE_Arith 3827 if ((word0(&u) & Exp_mask) == Exp_mask) 3828#else 3829 if (word0(&u) == 0x8000) 3830#endif 3831 { 3832 /* Infinity or NaN */ 3833 *decpt = 9999; 3834#ifdef IEEE_Arith 3835 if (!word1(&u) && !(word0(&u) & 0xfffff)) 3836 return nrv_alloc("Infinity", rve, 8); 3837#endif 3838 return nrv_alloc("NaN", rve, 3); 3839 } 3840#endif 3841#ifdef IBM 3842 dval(&u) += 0; /* normalize */ 3843#endif 3844 if (!dval(&u)) { 3845 *decpt = 1; 3846 return nrv_alloc("0", rve, 1); 3847 } 3848 3849#ifdef SET_INEXACT 3850 try_quick = oldinexact = get_inexact(); 3851 inexact = 1; 3852#endif 3853#ifdef Honor_FLT_ROUNDS 3854 if (Rounding >= 2) { 3855 if (*sign) 3856 Rounding = Rounding == 2 ? 0 : 2; 3857 else 3858 if (Rounding != 2) 3859 Rounding = 0; 3860 } 3861#endif 3862 3863 b = d2b(&u, &be, &bbits); 3864#ifdef Sudden_Underflow 3865 i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); 3866#else 3867 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) { 3868#endif 3869 dval(&d2) = dval(&u); 3870 word0(&d2) &= Frac_mask1; 3871 word0(&d2) |= Exp_11; 3872#ifdef IBM 3873 if (j = 11 - hi0bits(word0(&d2) & Frac_mask)) 3874 dval(&d2) /= 1 << j; 3875#endif 3876 3877 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 3878 * log10(x) = log(x) / log(10) 3879 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 3880 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 3881 * 3882 * This suggests computing an approximation k to log10(d) by 3883 * 3884 * k = (i - Bias)*0.301029995663981 3885 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 3886 * 3887 * We want k to be too large rather than too small. 3888 * The error in the first-order Taylor series approximation 3889 * is in our favor, so we just round up the constant enough 3890 * to compensate for any error in the multiplication of 3891 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 3892 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 3893 * adding 1e-13 to the constant term more than suffices. 3894 * Hence we adjust the constant term to 0.1760912590558. 3895 * (We could get a more accurate k by invoking log10, 3896 * but this is probably not worthwhile.) 3897 */ 3898 3899 i -= Bias; 3900#ifdef IBM 3901 i <<= 2; 3902 i += j; 3903#endif 3904#ifndef Sudden_Underflow 3905 denorm = 0; 3906 } 3907 else { 3908 /* d is denormalized */ 3909 3910 i = bbits + be + (Bias + (P-1) - 1); 3911 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32) 3912 : word1(&u) << (32 - i); 3913 dval(&d2) = x; 3914 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */ 3915 i -= (Bias + (P-1) - 1) + 1; 3916 denorm = 1; 3917 } 3918#endif 3919 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; 3920 k = (int)ds; 3921 if (ds < 0. && ds != k) 3922 k--; /* want k = floor(ds) */ 3923 k_check = 1; 3924 if (k >= 0 && k <= Ten_pmax) { 3925 if (dval(&u) < tens[k]) 3926 k--; 3927 k_check = 0; 3928 } 3929 j = bbits - i - 1; 3930 if (j >= 0) { 3931 b2 = 0; 3932 s2 = j; 3933 } 3934 else { 3935 b2 = -j; 3936 s2 = 0; 3937 } 3938 if (k >= 0) { 3939 b5 = 0; 3940 s5 = k; 3941 s2 += k; 3942 } 3943 else { 3944 b2 -= k; 3945 b5 = -k; 3946 s5 = 0; 3947 } 3948 if (mode < 0 || mode > 9) 3949 mode = 0; 3950 3951#ifndef SET_INEXACT 3952#ifdef Check_FLT_ROUNDS 3953 try_quick = Rounding == 1; 3954#else 3955 try_quick = 1; 3956#endif 3957#endif /*SET_INEXACT*/ 3958 3959 if (mode > 5) { 3960 mode -= 4; 3961 try_quick = 0; 3962 } 3963 leftright = 1; 3964 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */ 3965 /* silence erroneous "gcc -Wall" warning. */ 3966 switch(mode) { 3967 case 0: 3968 case 1: 3969 i = 18; 3970 ndigits = 0; 3971 break; 3972 case 2: 3973 leftright = 0; 3974 /* no break */ 3975 case 4: 3976 if (ndigits <= 0) 3977 ndigits = 1; 3978 ilim = ilim1 = i = ndigits; 3979 break; 3980 case 3: 3981 leftright = 0; 3982 /* no break */ 3983 case 5: 3984 i = ndigits + k + 1; 3985 ilim = i; 3986 ilim1 = i - 1; 3987 if (i <= 0) 3988 i = 1; 3989 } 3990 s = s0 = rv_alloc(i); 3991 3992#ifdef Honor_FLT_ROUNDS 3993 if (mode > 1 && Rounding != 1) 3994 leftright = 0; 3995#endif 3996 3997 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 3998 3999 /* Try to get by with floating-point arithmetic. */ 4000 4001 i = 0; 4002 dval(&d2) = dval(&u); 4003 k0 = k; 4004 ilim0 = ilim; 4005 ieps = 2; /* conservative */ 4006 if (k > 0) { 4007 ds = tens[k&0xf]; 4008 j = k >> 4; 4009 if (j & Bletch) { 4010 /* prevent overflows */ 4011 j &= Bletch - 1; 4012 dval(&u) /= bigtens[n_bigtens-1]; 4013 ieps++; 4014 } 4015 for(; j; j >>= 1, i++) 4016 if (j & 1) { 4017 ieps++; 4018 ds *= bigtens[i]; 4019 } 4020 dval(&u) /= ds; 4021 } 4022 else if ((j1 = -k)) { 4023 dval(&u) *= tens[j1 & 0xf]; 4024 for(j = j1 >> 4; j; j >>= 1, i++) 4025 if (j & 1) { 4026 ieps++; 4027 dval(&u) *= bigtens[i]; 4028 } 4029 } 4030 if (k_check && dval(&u) < 1. && ilim > 0) { 4031 if (ilim1 <= 0) 4032 goto fast_failed; 4033 ilim = ilim1; 4034 k--; 4035 dval(&u) *= 10.; 4036 ieps++; 4037 } 4038 dval(&eps) = ieps*dval(&u) + 7.; 4039 word0(&eps) -= (P-1)*Exp_msk1; 4040 if (ilim == 0) { 4041 S = mhi = 0; 4042 dval(&u) -= 5.; 4043 if (dval(&u) > dval(&eps)) 4044 goto one_digit; 4045 if (dval(&u) < -dval(&eps)) 4046 goto no_digits; 4047 goto fast_failed; 4048 } 4049#ifndef No_leftright 4050 if (leftright) { 4051 /* Use Steele & White method of only 4052 * generating digits needed. 4053 */ 4054 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps); 4055#ifdef IEEE_Arith 4056 if (k0 < 0 && j1 >= 307) { 4057 eps1.d = 1.01e256; /* 1.01 allows roundoff in the next few lines */ 4058 word0(&eps1) -= Exp_msk1 * (Bias+P-1); 4059 dval(&eps1) *= tens[j1 & 0xf]; 4060 for(i = 0, j = (j1-256) >> 4; j; j >>= 1, i++) 4061 if (j & 1) 4062 dval(&eps1) *= bigtens[i]; 4063 if (eps.d < eps1.d) 4064 eps.d = eps1.d; 4065 } 4066#endif 4067 for(i = 0;;) { 4068 L = dval(&u); 4069 dval(&u) -= L; 4070 *s++ = '0' + (int)L; 4071 if (1. - dval(&u) < dval(&eps)) 4072 goto bump_up; 4073 if (dval(&u) < dval(&eps)) 4074 goto ret1; 4075 if (++i >= ilim) 4076 break; 4077 dval(&eps) *= 10.; 4078 dval(&u) *= 10.; 4079 } 4080 } 4081 else { 4082#endif 4083 /* Generate ilim digits, then fix them up. */ 4084 dval(&eps) *= tens[ilim-1]; 4085 for(i = 1;; i++, dval(&u) *= 10.) { 4086 L = (Long)(dval(&u)); 4087 if (!(dval(&u) -= L)) 4088 ilim = i; 4089 *s++ = '0' + (int)L; 4090 if (i == ilim) { 4091 if (dval(&u) > 0.5 + dval(&eps)) 4092 goto bump_up; 4093 else if (dval(&u) < 0.5 - dval(&eps)) { 4094 while(*--s == '0'); 4095 s++; 4096 goto ret1; 4097 } 4098 break; 4099 } 4100 } 4101#ifndef No_leftright 4102 } 4103#endif 4104 fast_failed: 4105 s = s0; 4106 dval(&u) = dval(&d2); 4107 k = k0; 4108 ilim = ilim0; 4109 } 4110 4111 /* Do we have a "small" integer? */ 4112 4113 if (be >= 0 && k <= Int_max) { 4114 /* Yes. */ 4115 ds = tens[k]; 4116 if (ndigits < 0 && ilim <= 0) { 4117 S = mhi = 0; 4118 if (ilim < 0 || dval(&u) <= 5*ds) 4119 goto no_digits; 4120 goto one_digit; 4121 } 4122 for(i = 1;; i++, dval(&u) *= 10.) { 4123 L = (Long)(dval(&u) / ds); 4124 dval(&u) -= L*ds; 4125#ifdef Check_FLT_ROUNDS 4126 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 4127 if (dval(&u) < 0) { 4128 L--; 4129 dval(&u) += ds; 4130 } 4131#endif 4132 *s++ = '0' + (int)L; 4133 if (!dval(&u)) { 4134#ifdef SET_INEXACT 4135 inexact = 0; 4136#endif 4137 break; 4138 } 4139 if (i == ilim) { 4140#ifdef Honor_FLT_ROUNDS 4141 if (mode > 1) 4142 switch(Rounding) { 4143 case 0: goto ret1; 4144 case 2: goto bump_up; 4145 } 4146#endif 4147 dval(&u) += dval(&u); 4148#ifdef ROUND_BIASED 4149 if (dval(&u) >= ds) 4150#else 4151 if (dval(&u) > ds || (dval(&u) == ds && L & 1)) 4152#endif 4153 { 4154 bump_up: 4155 while(*--s == '9') 4156 if (s == s0) { 4157 k++; 4158 *s = '0'; 4159 break; 4160 } 4161 ++*s++; 4162 } 4163 break; 4164 } 4165 } 4166 goto ret1; 4167 } 4168 4169 m2 = b2; 4170 m5 = b5; 4171 mhi = mlo = 0; 4172 if (leftright) { 4173 i = 4174#ifndef Sudden_Underflow 4175 denorm ? be + (Bias + (P-1) - 1 + 1) : 4176#endif 4177#ifdef IBM 4178 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 4179#else 4180 1 + P - bbits; 4181#endif 4182 b2 += i; 4183 s2 += i; 4184 mhi = i2b(1); 4185 } 4186 if (m2 > 0 && s2 > 0) { 4187 i = m2 < s2 ? m2 : s2; 4188 b2 -= i; 4189 m2 -= i; 4190 s2 -= i; 4191 } 4192 if (b5 > 0) { 4193 if (leftright) { 4194 if (m5 > 0) { 4195 mhi = pow5mult(mhi, m5); 4196 b1 = mult(mhi, b); 4197 Bfree(b); 4198 b = b1; 4199 } 4200 if ((j = b5 - m5)) 4201 b = pow5mult(b, j); 4202 } 4203 else 4204 b = pow5mult(b, b5); 4205 } 4206 S = i2b(1); 4207 if (s5 > 0) 4208 S = pow5mult(S, s5); 4209 4210 /* Check for special case that d is a normalized power of 2. */ 4211 4212 spec_case = 0; 4213 if ((mode < 2 || leftright) 4214#ifdef Honor_FLT_ROUNDS 4215 && Rounding == 1 4216#endif 4217 ) { 4218 if (!word1(&u) && !(word0(&u) & Bndry_mask) 4219#ifndef Sudden_Underflow 4220 && word0(&u) & (Exp_mask & ~Exp_msk1) 4221#endif 4222 ) { 4223 /* The special case */ 4224 b2 += Log2P; 4225 s2 += Log2P; 4226 spec_case = 1; 4227 } 4228 } 4229 4230 /* Arrange for convenient computation of quotients: 4231 * shift left if necessary so divisor has 4 leading 0 bits. 4232 * 4233 * Perhaps we should just compute leading 28 bits of S once 4234 * and for all and pass them and a shift to quorem, so it 4235 * can do shifts and ors to compute the numerator for q. 4236 */ 4237 i = dshift(S, s2); 4238 b2 += i; 4239 m2 += i; 4240 s2 += i; 4241 if (b2 > 0) 4242 b = lshift(b, b2); 4243 if (s2 > 0) 4244 S = lshift(S, s2); 4245 if (k_check) { 4246 if (cmp(b,S) < 0) { 4247 k--; 4248 b = multadd(b, 10, 0); /* we botched the k estimate */ 4249 if (leftright) 4250 mhi = multadd(mhi, 10, 0); 4251 ilim = ilim1; 4252 } 4253 } 4254 if (ilim <= 0 && (mode == 3 || mode == 5)) { 4255 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 4256 /* no digits, fcvt style */ 4257 no_digits: 4258 k = -1 - ndigits; 4259 goto ret; 4260 } 4261 one_digit: 4262 *s++ = '1'; 4263 k++; 4264 goto ret; 4265 } 4266 if (leftright) { 4267 if (m2 > 0) 4268 mhi = lshift(mhi, m2); 4269 4270 /* Compute mlo -- check for special case 4271 * that d is a normalized power of 2. 4272 */ 4273 4274 mlo = mhi; 4275 if (spec_case) { 4276 mhi = Balloc(mhi->k); 4277 Bcopy(mhi, mlo); 4278 mhi = lshift(mhi, Log2P); 4279 } 4280 4281 for(i = 1;;i++) { 4282 dig = quorem(b,S) + '0'; 4283 /* Do we yet have the shortest decimal string 4284 * that will round to d? 4285 */ 4286 j = cmp(b, mlo); 4287 delta = diff(S, mhi); 4288 j1 = delta->sign ? 1 : cmp(b, delta); 4289 Bfree(delta); 4290#ifndef ROUND_BIASED 4291 if (j1 == 0 && mode != 1 && !(word1(&u) & 1) 4292#ifdef Honor_FLT_ROUNDS 4293 && Rounding >= 1 4294#endif 4295 ) { 4296 if (dig == '9') 4297 goto round_9_up; 4298 if (j > 0) 4299 dig++; 4300#ifdef SET_INEXACT 4301 else if (!b->x[0] && b->wds <= 1) 4302 inexact = 0; 4303#endif 4304 *s++ = dig; 4305 goto ret; 4306 } 4307#endif 4308 if (j < 0 || (j == 0 && mode != 1 4309#ifndef ROUND_BIASED 4310 && !(word1(&u) & 1) 4311#endif 4312 )) { 4313 if (!b->x[0] && b->wds <= 1) { 4314#ifdef SET_INEXACT 4315 inexact = 0; 4316#endif 4317 goto accept_dig; 4318 } 4319#ifdef Honor_FLT_ROUNDS 4320 if (mode > 1) 4321 switch(Rounding) { 4322 case 0: goto accept_dig; 4323 case 2: goto keep_dig; 4324 } 4325#endif /*Honor_FLT_ROUNDS*/ 4326 if (j1 > 0) { 4327 b = lshift(b, 1); 4328 j1 = cmp(b, S); 4329#ifdef ROUND_BIASED 4330 if (j1 >= 0 /*)*/ 4331#else 4332 if ((j1 > 0 || (j1 == 0 && dig & 1)) 4333#endif 4334 && dig++ == '9') 4335 goto round_9_up; 4336 } 4337 accept_dig: 4338 *s++ = dig; 4339 goto ret; 4340 } 4341 if (j1 > 0) { 4342#ifdef Honor_FLT_ROUNDS 4343 if (!Rounding) 4344 goto accept_dig; 4345#endif 4346 if (dig == '9') { /* possible if i == 1 */ 4347 round_9_up: 4348 *s++ = '9'; 4349 goto roundoff; 4350 } 4351 *s++ = dig + 1; 4352 goto ret; 4353 } 4354#ifdef Honor_FLT_ROUNDS 4355 keep_dig: 4356#endif 4357 *s++ = dig; 4358 if (i == ilim) 4359 break; 4360 b = multadd(b, 10, 0); 4361 if (mlo == mhi) 4362 mlo = mhi = multadd(mhi, 10, 0); 4363 else { 4364 mlo = multadd(mlo, 10, 0); 4365 mhi = multadd(mhi, 10, 0); 4366 } 4367 } 4368 } 4369 else 4370 for(i = 1;; i++) { 4371 *s++ = dig = quorem(b,S) + '0'; 4372 if (!b->x[0] && b->wds <= 1) { 4373#ifdef SET_INEXACT 4374 inexact = 0; 4375#endif 4376 goto ret; 4377 } 4378 if (i >= ilim) 4379 break; 4380 b = multadd(b, 10, 0); 4381 } 4382 4383 /* Round off last digit */ 4384 4385#ifdef Honor_FLT_ROUNDS 4386 switch(Rounding) { 4387 case 0: goto trimzeros; 4388 case 2: goto roundoff; 4389 } 4390#endif 4391 b = lshift(b, 1); 4392 j = cmp(b, S); 4393#ifdef ROUND_BIASED 4394 if (j >= 0) 4395#else 4396 if (j > 0 || (j == 0 && dig & 1)) 4397#endif 4398 { 4399 roundoff: 4400 while(*--s == '9') 4401 if (s == s0) { 4402 k++; 4403 *s++ = '1'; 4404 goto ret; 4405 } 4406 ++*s++; 4407 } 4408 else { 4409#ifdef Honor_FLT_ROUNDS 4410 trimzeros: 4411#endif 4412 while(*--s == '0'); 4413 s++; 4414 } 4415 ret: 4416 Bfree(S); 4417 if (mhi) { 4418 if (mlo && mlo != mhi) 4419 Bfree(mlo); 4420 Bfree(mhi); 4421 } 4422 ret1: 4423#ifdef SET_INEXACT 4424 if (inexact) { 4425 if (!oldinexact) { 4426 word0(&u) = Exp_1 + (70 << Exp_shift); 4427 word1(&u) = 0; 4428 dval(&u) += 1.; 4429 } 4430 } 4431 else if (!oldinexact) 4432 clear_inexact(); 4433#endif 4434 Bfree(b); 4435 *s = 0; 4436 *decpt = k + 1; 4437 if (rve) 4438 *rve = s; 4439 return s0; 4440 } 4441#ifdef __cplusplus 4442} 4443#endif 4444