1 /* $NetBSD: strtod.c,v 1.45.2.1 2005/04/19 13:35:54 tron Exp $ */ 2 3 /**************************************************************** 4 * 5 * The author of this software is David M. Gay. 6 * 7 * Copyright (c) 1991 by AT&T. 8 * 9 * Permission to use, copy, modify, and distribute this software for any 10 * purpose without fee is hereby granted, provided that this entire notice 11 * is included in all copies of any software which is or includes a copy 12 * or modification of this software and in all copies of the supporting 13 * documentation for such software. 14 * 15 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 16 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY 17 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 18 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 19 * 20 ***************************************************************/ 21 22 /* Please send bug reports to 23 David M. Gay 24 AT&T Bell Laboratories, Room 2C-463 25 600 Mountain Avenue 26 Murray Hill, NJ 07974-2070 27 U.S.A. 28 dmg (at) research.att.com or research!dmg 29 */ 30 31 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 32 * 33 * This strtod returns a nearest machine number to the input decimal 34 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 35 * broken by the IEEE round-even rule. Otherwise ties are broken by 36 * biased rounding (add half and chop). 37 * 38 * Inspired loosely by William D. Clinger's paper "How to Read Floating 39 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 40 * 41 * Modifications: 42 * 43 * 1. We only require IEEE, IBM, or VAX double-precision 44 * arithmetic (not IEEE double-extended). 45 * 2. We get by with floating-point arithmetic in a case that 46 * Clinger missed -- when we're computing d * 10^n 47 * for a small integer d and the integer n is not too 48 * much larger than 22 (the maximum integer k for which 49 * we can represent 10^k exactly), we may be able to 50 * compute (d*10^k) * 10^(e-k) with just one roundoff. 51 * 3. Rather than a bit-at-a-time adjustment of the binary 52 * result in the hard case, we use floating-point 53 * arithmetic to determine the adjustment to within 54 * one bit; only in really hard cases do we need to 55 * compute a second residual. 56 * 4. Because of 3., we don't need a large table of powers of 10 57 * for ten-to-e (just some small tables, e.g. of 10^k 58 * for 0 <= k <= 22). 59 */ 60 61 /* 62 * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least 63 * significant byte has the lowest address. 64 * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most 65 * significant byte has the lowest address. 66 * #define Long int on machines with 32-bit ints and 64-bit longs. 67 * #define Sudden_Underflow for IEEE-format machines without gradual 68 * underflow (i.e., that flush to zero on underflow). 69 * #define IBM for IBM mainframe-style floating-point arithmetic. 70 * #define VAX for VAX-style floating-point arithmetic. 71 * #define Unsigned_Shifts if >> does treats its left operand as unsigned. 72 * #define No_leftright to omit left-right logic in fast floating-point 73 * computation of dtoa. 74 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3. 75 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 76 * that use extended-precision instructions to compute rounded 77 * products and quotients) with IBM. 78 * #define ROUND_BIASED for IEEE-format with biased rounding. 79 * #define Inaccurate_Divide for IEEE-format with correctly rounded 80 * products but inaccurate quotients, e.g., for Intel i860. 81 * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision 82 * integer arithmetic. Whether this speeds things up or slows things 83 * down depends on the machine and the number being converted. 84 * #define KR_headers for old-style C function headers. 85 * #define Bad_float_h if your system lacks a float.h or if it does not 86 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, 87 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. 88 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) 89 * if memory is available and otherwise does something you deem 90 * appropriate. If MALLOC is undefined, malloc will be invoked 91 * directly -- and assumed always to succeed. 92 */ 93 94 #define ANDROID_CHANGES 95 96 #ifdef ANDROID_CHANGES 97 /* Needs to be above math.h include below */ 98 #include "fpmath.h" 99 100 #include <pthread.h> 101 #define mutex_lock(x) pthread_mutex_lock(x) 102 #define mutex_unlock(x) pthread_mutex_unlock(x) 103 #endif 104 105 #include <sys/cdefs.h> 106 #if defined(LIBC_SCCS) && !defined(lint) 107 __RCSID("$NetBSD: strtod.c,v 1.45.2.1 2005/04/19 13:35:54 tron Exp $"); 108 #endif /* LIBC_SCCS and not lint */ 109 110 #define Unsigned_Shifts 111 #if defined(__m68k__) || defined(__sparc__) || defined(__i386__) || \ 112 defined(__mips__) || defined(__ns32k__) || defined(__alpha__) || \ 113 defined(__powerpc__) || defined(__sh__) || defined(__x86_64__) || \ 114 defined(__hppa__) || \ 115 (defined(__arm__) && defined(__VFP_FP__)) || defined(__aarch64__) || \ 116 defined(__le32__) || defined(__le64__) 117 #include <endian.h> 118 #if BYTE_ORDER == BIG_ENDIAN 119 #define IEEE_BIG_ENDIAN 120 #else 121 #define IEEE_LITTLE_ENDIAN 122 #endif 123 #endif 124 125 #if defined(__arm__) && !defined(__VFP_FP__) 126 /* 127 * Although the CPU is little endian the FP has different 128 * byte and word endianness. The byte order is still little endian 129 * but the word order is big endian. 130 */ 131 #define IEEE_BIG_ENDIAN 132 #endif 133 134 #ifdef __vax__ 135 #define VAX 136 #endif 137 138 #if defined(__hppa__) || defined(__mips__) || defined(__sh__) 139 #define NAN_WORD0 0x7ff40000 140 #else 141 #define NAN_WORD0 0x7ff80000 142 #endif 143 #define NAN_WORD1 0 144 145 #define Long int32_t 146 #define ULong u_int32_t 147 148 #ifdef DEBUG 149 #include "stdio.h" 150 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 151 #define BugPrintf(x, v) {fprintf(stderr, x, v); exit(1);} 152 #endif 153 154 #ifdef __cplusplus 155 #include "malloc.h" 156 #include "memory.h" 157 #else 158 #ifndef KR_headers 159 #include "stdlib.h" 160 #include "string.h" 161 #ifndef ANDROID_CHANGES 162 #include "locale.h" 163 #endif /* ANDROID_CHANGES */ 164 #else 165 #include "malloc.h" 166 #include "memory.h" 167 #endif 168 #endif 169 #ifndef ANDROID_CHANGES 170 #include "extern.h" 171 #include "reentrant.h" 172 #endif /* ANDROID_CHANGES */ 173 174 #ifdef MALLOC 175 #ifdef KR_headers 176 extern char *MALLOC(); 177 #else 178 extern void *MALLOC(size_t); 179 #endif 180 #else 181 #define MALLOC malloc 182 #endif 183 184 #include "ctype.h" 185 #include "errno.h" 186 #include "float.h" 187 188 #ifndef __MATH_H__ 189 #include "math.h" 190 #endif 191 192 #ifdef __cplusplus 193 extern "C" { 194 #endif 195 196 #ifndef CONST 197 #ifdef KR_headers 198 #define CONST /* blank */ 199 #else 200 #define CONST const 201 #endif 202 #endif 203 204 #ifdef Unsigned_Shifts 205 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000; 206 #else 207 #define Sign_Extend(a,b) /*no-op*/ 208 #endif 209 210 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + \ 211 defined(IBM) != 1 212 Exactly one of IEEE_LITTLE_ENDIAN IEEE_BIG_ENDIAN, VAX, or 213 IBM should be defined. 214 #endif 215 216 typedef union { 217 double d; 218 ULong ul[2]; 219 } _double; 220 #define value(x) ((x).d) 221 #ifdef IEEE_LITTLE_ENDIAN 222 #define word0(x) ((x).ul[1]) 223 #define word1(x) ((x).ul[0]) 224 #else 225 #define word0(x) ((x).ul[0]) 226 #define word1(x) ((x).ul[1]) 227 #endif 228 229 /* The following definition of Storeinc is appropriate for MIPS processors. 230 * An alternative that might be better on some machines is 231 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 232 */ 233 #if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm__) 234 #define Storeinc(a,b,c) \ 235 (((u_short *)(void *)a)[1] = \ 236 (u_short)b, ((u_short *)(void *)a)[0] = (u_short)c, a++) 237 #else 238 #define Storeinc(a,b,c) \ 239 (((u_short *)(void *)a)[0] = \ 240 (u_short)b, ((u_short *)(void *)a)[1] = (u_short)c, a++) 241 #endif 242 243 /* #define P DBL_MANT_DIG */ 244 /* Ten_pmax = floor(P*log(2)/log(5)) */ 245 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 246 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 247 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 248 249 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) 250 #define Exp_shift 20 251 #define Exp_shift1 20 252 #define Exp_msk1 0x100000 253 #define Exp_msk11 0x100000 254 #define Exp_mask 0x7ff00000 255 #define P 53 256 #define Bias 1023 257 #define IEEE_Arith 258 #define Emin (-1022) 259 #define Exp_1 0x3ff00000 260 #define Exp_11 0x3ff00000 261 #define Ebits 11 262 #define Frac_mask 0xfffff 263 #define Frac_mask1 0xfffff 264 #define Ten_pmax 22 265 #define Bletch 0x10 266 #define Bndry_mask 0xfffff 267 #define Bndry_mask1 0xfffff 268 #define LSB 1 269 #define Sign_bit 0x80000000 270 #define Log2P 1 271 #define Tiny0 0 272 #define Tiny1 1 273 #define Quick_max 14 274 #define Int_max 14 275 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */ 276 #else 277 #undef Sudden_Underflow 278 #define Sudden_Underflow 279 #ifdef IBM 280 #define Exp_shift 24 281 #define Exp_shift1 24 282 #define Exp_msk1 0x1000000 283 #define Exp_msk11 0x1000000 284 #define Exp_mask 0x7f000000 285 #define P 14 286 #define Bias 65 287 #define Exp_1 0x41000000 288 #define Exp_11 0x41000000 289 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 290 #define Frac_mask 0xffffff 291 #define Frac_mask1 0xffffff 292 #define Bletch 4 293 #define Ten_pmax 22 294 #define Bndry_mask 0xefffff 295 #define Bndry_mask1 0xffffff 296 #define LSB 1 297 #define Sign_bit 0x80000000 298 #define Log2P 4 299 #define Tiny0 0x100000 300 #define Tiny1 0 301 #define Quick_max 14 302 #define Int_max 15 303 #else /* VAX */ 304 #define Exp_shift 23 305 #define Exp_shift1 7 306 #define Exp_msk1 0x80 307 #define Exp_msk11 0x800000 308 #define Exp_mask 0x7f80 309 #define P 56 310 #define Bias 129 311 #define Exp_1 0x40800000 312 #define Exp_11 0x4080 313 #define Ebits 8 314 #define Frac_mask 0x7fffff 315 #define Frac_mask1 0xffff007f 316 #define Ten_pmax 24 317 #define Bletch 2 318 #define Bndry_mask 0xffff007f 319 #define Bndry_mask1 0xffff007f 320 #define LSB 0x10000 321 #define Sign_bit 0x8000 322 #define Log2P 1 323 #define Tiny0 0x80 324 #define Tiny1 0 325 #define Quick_max 15 326 #define Int_max 15 327 #endif 328 #endif 329 330 #ifndef IEEE_Arith 331 #define ROUND_BIASED 332 #endif 333 334 #ifdef RND_PRODQUOT 335 #define rounded_product(a,b) a = rnd_prod(a, b) 336 #define rounded_quotient(a,b) a = rnd_quot(a, b) 337 #ifdef KR_headers 338 extern double rnd_prod(), rnd_quot(); 339 #else 340 extern double rnd_prod(double, double), rnd_quot(double, double); 341 #endif 342 #else 343 #define rounded_product(a,b) a *= b 344 #define rounded_quotient(a,b) a /= b 345 #endif 346 347 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 348 #define Big1 0xffffffff 349 350 #ifndef Just_16 351 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long. 352 * This makes some inner loops simpler and sometimes saves work 353 * during multiplications, but it often seems to make things slightly 354 * slower. Hence the default is now to store 32 bits per Long. 355 */ 356 #ifndef Pack_32 357 #define Pack_32 358 #endif 359 #endif 360 361 #define Kmax 15 362 363 #ifdef Pack_32 364 #define ULbits 32 365 #define kshift 5 366 #define kmask 31 367 #define ALL_ON 0xffffffff 368 #else 369 #define ULbits 16 370 #define kshift 4 371 #define kmask 15 372 #define ALL_ON 0xffff 373 #endif 374 375 #define Kmax 15 376 377 enum { /* return values from strtodg */ 378 STRTOG_Zero = 0, 379 STRTOG_Normal = 1, 380 STRTOG_Denormal = 2, 381 STRTOG_Infinite = 3, 382 STRTOG_NaN = 4, 383 STRTOG_NaNbits = 5, 384 STRTOG_NoNumber = 6, 385 STRTOG_Retmask = 7, 386 387 /* The following may be or-ed into one of the above values. */ 388 389 STRTOG_Neg = 0x08, /* does not affect STRTOG_Inexlo or STRTOG_Inexhi */ 390 STRTOG_Inexlo = 0x10, /* returned result rounded toward zero */ 391 STRTOG_Inexhi = 0x20, /* returned result rounded away from zero */ 392 STRTOG_Inexact = 0x30, 393 STRTOG_Underflow= 0x40, 394 STRTOG_Overflow = 0x80 395 }; 396 397 typedef struct 398 FPI { 399 int nbits; 400 int emin; 401 int emax; 402 int rounding; 403 int sudden_underflow; 404 } FPI; 405 406 enum { /* FPI.rounding values: same as FLT_ROUNDS */ 407 FPI_Round_zero = 0, 408 FPI_Round_near = 1, 409 FPI_Round_up = 2, 410 FPI_Round_down = 3 411 }; 412 413 #undef SI 414 #ifdef Sudden_Underflow 415 #define SI 1 416 #else 417 #define SI 0 418 #endif 419 420 #ifdef __cplusplus 421 extern "C" double strtod(const char *s00, char **se); 422 extern "C" char *__dtoa(double d, int mode, int ndigits, 423 int *decpt, int *sign, char **rve); 424 #endif 425 426 struct 427 Bigint { 428 struct Bigint *next; 429 int k, maxwds, sign, wds; 430 ULong x[1]; 431 }; 432 433 typedef struct Bigint Bigint; 434 435 CONST unsigned char hexdig[256] = { 436 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 437 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 438 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 439 0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0x18, 0x19, 0, 0, 0, 0, 0, 0, 440 0, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f, 0, 0, 0, 0, 0, 0, 0, 0, 0, 441 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 442 0, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f, 0, 0, 0, 0, 0, 0, 0, 0, 0, 443 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 444 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 445 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 446 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 447 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 448 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 449 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 450 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 451 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 452 }; 453 454 static int 455 gethex(CONST char **, CONST FPI *, Long *, Bigint **, int, locale_t); 456 457 458 static Bigint *freelist[Kmax+1]; 459 460 #ifdef ANDROID_CHANGES 461 static pthread_mutex_t freelist_mutex = PTHREAD_MUTEX_INITIALIZER; 462 #else 463 #ifdef _REENTRANT 464 static mutex_t freelist_mutex = MUTEX_INITIALIZER; 465 #endif 466 #endif 467 468 /* Special value used to indicate an invalid Bigint value, 469 * e.g. when a memory allocation fails. The idea is that we 470 * want to avoid introducing NULL checks everytime a bigint 471 * computation is performed. Also the NULL value can also be 472 * already used to indicate "value not initialized yet" and 473 * returning NULL might alter the execution code path in 474 * case of OOM. 475 */ 476 #define BIGINT_INVALID ((Bigint *)&bigint_invalid_value) 477 478 static const Bigint bigint_invalid_value; 479 480 481 static void 482 copybits(ULong *c, int n, Bigint *b) 483 { 484 ULong *ce, *x, *xe; 485 #ifdef Pack_16 486 int nw, nw1; 487 #endif 488 489 ce = c + ((n-1) >> kshift) + 1; 490 x = b->x; 491 #ifdef Pack_32 492 xe = x + b->wds; 493 while(x < xe) 494 *c++ = *x++; 495 #else 496 nw = b->wds; 497 nw1 = nw & 1; 498 for(xe = x + (nw - nw1); x < xe; x += 2) 499 Storeinc(c, x[1], x[0]); 500 if (nw1) 501 *c++ = *x; 502 #endif 503 while(c < ce) 504 *c++ = 0; 505 } 506 507 ULong 508 any_on(Bigint *b, int k) 509 { 510 int n, nwds; 511 ULong *x, *x0, x1, x2; 512 513 x = b->x; 514 nwds = b->wds; 515 n = k >> kshift; 516 if (n > nwds) 517 n = nwds; 518 else if (n < nwds && (k &= kmask)) { 519 x1 = x2 = x[n]; 520 x1 >>= k; 521 x1 <<= k; 522 if (x1 != x2) 523 return 1; 524 } 525 x0 = x; 526 x += n; 527 while(x > x0) 528 if (*--x) 529 return 1; 530 return 0; 531 } 532 533 void 534 rshift(Bigint *b, int k) 535 { 536 ULong *x, *x1, *xe, y; 537 int n; 538 539 x = x1 = b->x; 540 n = k >> kshift; 541 if (n < b->wds) { 542 xe = x + b->wds; 543 x += n; 544 if (k &= kmask) { 545 n = ULbits - k; 546 y = *x++ >> k; 547 while(x < xe) { 548 *x1++ = (y | (*x << n)) & ALL_ON; 549 y = *x++ >> k; 550 } 551 if ((*x1 = y) !=0) 552 x1++; 553 } 554 else 555 while(x < xe) 556 *x1++ = *x++; 557 } 558 if ((b->wds = x1 - b->x) == 0) 559 b->x[0] = 0; 560 } 561 562 563 typedef union { double d; ULong L[2]; } U; 564 565 static void 566 ULtod(ULong *L, ULong *bits, Long exp, int k) 567 { 568 # define _0 1 569 # define _1 0 570 571 switch(k & STRTOG_Retmask) { 572 case STRTOG_NoNumber: 573 case STRTOG_Zero: 574 L[0] = L[1] = 0; 575 break; 576 577 case STRTOG_Denormal: 578 L[_1] = bits[0]; 579 L[_0] = bits[1]; 580 break; 581 582 case STRTOG_Normal: 583 case STRTOG_NaNbits: 584 L[_1] = bits[0]; 585 L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20); 586 break; 587 588 case STRTOG_Infinite: 589 L[_0] = 0x7ff00000; 590 L[_1] = 0; 591 break; 592 593 #define d_QNAN0 0x7ff80000 594 #define d_QNAN1 0x0 595 case STRTOG_NaN: 596 L[0] = d_QNAN0; 597 L[1] = d_QNAN1; 598 } 599 if (k & STRTOG_Neg) 600 L[_0] |= 0x80000000L; 601 } 602 603 604 605 /* Return BIGINT_INVALID on allocation failure. 606 * 607 * Most of the code here depends on the fact that this function 608 * never returns NULL. 609 */ 610 static Bigint * 611 Balloc 612 #ifdef KR_headers 613 (k) int k; 614 #else 615 (int k) 616 #endif 617 { 618 int x; 619 Bigint *rv; 620 621 mutex_lock(&freelist_mutex); 622 623 if ((rv = freelist[k]) != NULL) { 624 freelist[k] = rv->next; 625 } 626 else { 627 x = 1 << k; 628 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(Long)); 629 if (rv == NULL) { 630 rv = BIGINT_INVALID; 631 goto EXIT; 632 } 633 rv->k = k; 634 rv->maxwds = x; 635 } 636 rv->sign = rv->wds = 0; 637 EXIT: 638 mutex_unlock(&freelist_mutex); 639 640 return rv; 641 } 642 643 static void 644 Bfree 645 #ifdef KR_headers 646 (v) Bigint *v; 647 #else 648 (Bigint *v) 649 #endif 650 { 651 if (v && v != BIGINT_INVALID) { 652 mutex_lock(&freelist_mutex); 653 654 v->next = freelist[v->k]; 655 freelist[v->k] = v; 656 657 mutex_unlock(&freelist_mutex); 658 } 659 } 660 661 #define Bcopy_valid(x,y) memcpy(&(x)->sign, &(y)->sign, \ 662 (y)->wds*sizeof(Long) + 2*sizeof(int)) 663 664 #define Bcopy(x,y) Bcopy_ptr(&(x),(y)) 665 666 static void 667 Bcopy_ptr(Bigint **px, Bigint *y) 668 { 669 if (*px == BIGINT_INVALID) 670 return; /* no space to store copy */ 671 if (y == BIGINT_INVALID) { 672 Bfree(*px); /* invalid input */ 673 *px = BIGINT_INVALID; 674 } else { 675 Bcopy_valid(*px,y); 676 } 677 } 678 679 static Bigint * 680 multadd 681 #ifdef KR_headers 682 (b, m, a) Bigint *b; int m, a; 683 #else 684 (Bigint *b, int m, int a) /* multiply by m and add a */ 685 #endif 686 { 687 int i, wds; 688 ULong *x, y; 689 #ifdef Pack_32 690 ULong xi, z; 691 #endif 692 Bigint *b1; 693 694 if (b == BIGINT_INVALID) 695 return b; 696 697 wds = b->wds; 698 x = b->x; 699 i = 0; 700 do { 701 #ifdef Pack_32 702 xi = *x; 703 y = (xi & 0xffff) * m + a; 704 z = (xi >> 16) * m + (y >> 16); 705 a = (int)(z >> 16); 706 *x++ = (z << 16) + (y & 0xffff); 707 #else 708 y = *x * m + a; 709 a = (int)(y >> 16); 710 *x++ = y & 0xffff; 711 #endif 712 } 713 while(++i < wds); 714 if (a) { 715 if (wds >= b->maxwds) { 716 b1 = Balloc(b->k+1); 717 if (b1 == BIGINT_INVALID) { 718 Bfree(b); 719 return b1; 720 } 721 Bcopy_valid(b1, b); 722 Bfree(b); 723 b = b1; 724 } 725 b->x[wds++] = a; 726 b->wds = wds; 727 } 728 return b; 729 } 730 731 Bigint * 732 increment(Bigint *b) 733 { 734 ULong *x, *xe; 735 Bigint *b1; 736 #ifdef Pack_16 737 ULong carry = 1, y; 738 #endif 739 740 x = b->x; 741 xe = x + b->wds; 742 #ifdef Pack_32 743 do { 744 if (*x < (ULong)0xffffffffL) { 745 ++*x; 746 return b; 747 } 748 *x++ = 0; 749 } while(x < xe); 750 #else 751 do { 752 y = *x + carry; 753 carry = y >> 16; 754 *x++ = y & 0xffff; 755 if (!carry) 756 return b; 757 } while(x < xe); 758 if (carry) 759 #endif 760 { 761 if (b->wds >= b->maxwds) { 762 b1 = Balloc(b->k+1); 763 Bcopy(b1,b); 764 Bfree(b); 765 b = b1; 766 } 767 b->x[b->wds++] = 1; 768 } 769 return b; 770 } 771 772 773 static Bigint * 774 s2b 775 #ifdef KR_headers 776 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9; 777 #else 778 (CONST char *s, int nd0, int nd, ULong y9) 779 #endif 780 { 781 Bigint *b; 782 int i, k; 783 Long x, y; 784 785 x = (nd + 8) / 9; 786 for(k = 0, y = 1; x > y; y <<= 1, k++) ; 787 #ifdef Pack_32 788 b = Balloc(k); 789 if (b == BIGINT_INVALID) 790 return b; 791 b->x[0] = y9; 792 b->wds = 1; 793 #else 794 b = Balloc(k+1); 795 if (b == BIGINT_INVALID) 796 return b; 797 798 b->x[0] = y9 & 0xffff; 799 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; 800 #endif 801 802 i = 9; 803 if (9 < nd0) { 804 s += 9; 805 do b = multadd(b, 10, *s++ - '0'); 806 while(++i < nd0); 807 s++; 808 } 809 else 810 s += 10; 811 for(; i < nd; i++) 812 b = multadd(b, 10, *s++ - '0'); 813 return b; 814 } 815 816 static int 817 hi0bits 818 #ifdef KR_headers 819 (x) ULong x; 820 #else 821 (ULong x) 822 #endif 823 { 824 int k = 0; 825 826 if (!(x & 0xffff0000)) { 827 k = 16; 828 x <<= 16; 829 } 830 if (!(x & 0xff000000)) { 831 k += 8; 832 x <<= 8; 833 } 834 if (!(x & 0xf0000000)) { 835 k += 4; 836 x <<= 4; 837 } 838 if (!(x & 0xc0000000)) { 839 k += 2; 840 x <<= 2; 841 } 842 if (!(x & 0x80000000)) { 843 k++; 844 if (!(x & 0x40000000)) 845 return 32; 846 } 847 return k; 848 } 849 850 static int 851 lo0bits 852 #ifdef KR_headers 853 (y) ULong *y; 854 #else 855 (ULong *y) 856 #endif 857 { 858 int k; 859 ULong x = *y; 860 861 if (x & 7) { 862 if (x & 1) 863 return 0; 864 if (x & 2) { 865 *y = x >> 1; 866 return 1; 867 } 868 *y = x >> 2; 869 return 2; 870 } 871 k = 0; 872 if (!(x & 0xffff)) { 873 k = 16; 874 x >>= 16; 875 } 876 if (!(x & 0xff)) { 877 k += 8; 878 x >>= 8; 879 } 880 if (!(x & 0xf)) { 881 k += 4; 882 x >>= 4; 883 } 884 if (!(x & 0x3)) { 885 k += 2; 886 x >>= 2; 887 } 888 if (!(x & 1)) { 889 k++; 890 x >>= 1; 891 if (!x & 1) 892 return 32; 893 } 894 *y = x; 895 return k; 896 } 897 898 static Bigint * 899 i2b 900 #ifdef KR_headers 901 (i) int i; 902 #else 903 (int i) 904 #endif 905 { 906 Bigint *b; 907 908 b = Balloc(1); 909 if (b != BIGINT_INVALID) { 910 b->x[0] = i; 911 b->wds = 1; 912 } 913 return b; 914 } 915 916 static Bigint * 917 mult 918 #ifdef KR_headers 919 (a, b) Bigint *a, *b; 920 #else 921 (Bigint *a, Bigint *b) 922 #endif 923 { 924 Bigint *c; 925 int k, wa, wb, wc; 926 ULong carry, y, z; 927 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 928 #ifdef Pack_32 929 ULong z2; 930 #endif 931 932 if (a == BIGINT_INVALID || b == BIGINT_INVALID) 933 return BIGINT_INVALID; 934 935 if (a->wds < b->wds) { 936 c = a; 937 a = b; 938 b = c; 939 } 940 k = a->k; 941 wa = a->wds; 942 wb = b->wds; 943 wc = wa + wb; 944 if (wc > a->maxwds) 945 k++; 946 c = Balloc(k); 947 if (c == BIGINT_INVALID) 948 return c; 949 for(x = c->x, xa = x + wc; x < xa; x++) 950 *x = 0; 951 xa = a->x; 952 xae = xa + wa; 953 xb = b->x; 954 xbe = xb + wb; 955 xc0 = c->x; 956 #ifdef Pack_32 957 for(; xb < xbe; xb++, xc0++) { 958 if ((y = *xb & 0xffff) != 0) { 959 x = xa; 960 xc = xc0; 961 carry = 0; 962 do { 963 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 964 carry = z >> 16; 965 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 966 carry = z2 >> 16; 967 Storeinc(xc, z2, z); 968 } 969 while(x < xae); 970 *xc = carry; 971 } 972 if ((y = *xb >> 16) != 0) { 973 x = xa; 974 xc = xc0; 975 carry = 0; 976 z2 = *xc; 977 do { 978 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 979 carry = z >> 16; 980 Storeinc(xc, z, z2); 981 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 982 carry = z2 >> 16; 983 } 984 while(x < xae); 985 *xc = z2; 986 } 987 } 988 #else 989 for(; xb < xbe; xc0++) { 990 if (y = *xb++) { 991 x = xa; 992 xc = xc0; 993 carry = 0; 994 do { 995 z = *x++ * y + *xc + carry; 996 carry = z >> 16; 997 *xc++ = z & 0xffff; 998 } 999 while(x < xae); 1000 *xc = carry; 1001 } 1002 } 1003 #endif 1004 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 1005 c->wds = wc; 1006 return c; 1007 } 1008 1009 static Bigint *p5s; 1010 static pthread_mutex_t p5s_mutex = PTHREAD_MUTEX_INITIALIZER; 1011 1012 static Bigint * 1013 pow5mult 1014 #ifdef KR_headers 1015 (b, k) Bigint *b; int k; 1016 #else 1017 (Bigint *b, int k) 1018 #endif 1019 { 1020 Bigint *b1, *p5, *p51; 1021 int i; 1022 static const int p05[3] = { 5, 25, 125 }; 1023 1024 if (b == BIGINT_INVALID) 1025 return b; 1026 1027 if ((i = k & 3) != 0) 1028 b = multadd(b, p05[i-1], 0); 1029 1030 if (!(k = (unsigned int) k >> 2)) 1031 return b; 1032 mutex_lock(&p5s_mutex); 1033 if (!(p5 = p5s)) { 1034 /* first time */ 1035 p5 = i2b(625); 1036 if (p5 == BIGINT_INVALID) { 1037 Bfree(b); 1038 mutex_unlock(&p5s_mutex); 1039 return p5; 1040 } 1041 p5s = p5; 1042 p5->next = 0; 1043 } 1044 for(;;) { 1045 if (k & 1) { 1046 b1 = mult(b, p5); 1047 Bfree(b); 1048 b = b1; 1049 } 1050 if (!(k = (unsigned int) k >> 1)) 1051 break; 1052 if (!(p51 = p5->next)) { 1053 p51 = mult(p5,p5); 1054 if (p51 == BIGINT_INVALID) { 1055 Bfree(b); 1056 mutex_unlock(&p5s_mutex); 1057 return p51; 1058 } 1059 p5->next = p51; 1060 p51->next = 0; 1061 } 1062 p5 = p51; 1063 } 1064 mutex_unlock(&p5s_mutex); 1065 return b; 1066 } 1067 1068 static Bigint * 1069 lshift 1070 #ifdef KR_headers 1071 (b, k) Bigint *b; int k; 1072 #else 1073 (Bigint *b, int k) 1074 #endif 1075 { 1076 int i, k1, n, n1; 1077 Bigint *b1; 1078 ULong *x, *x1, *xe, z; 1079 1080 if (b == BIGINT_INVALID) 1081 return b; 1082 1083 #ifdef Pack_32 1084 n = (unsigned int)k >> 5; 1085 #else 1086 n = (unsigned int)k >> 4; 1087 #endif 1088 k1 = b->k; 1089 n1 = n + b->wds + 1; 1090 for(i = b->maxwds; n1 > i; i <<= 1) 1091 k1++; 1092 b1 = Balloc(k1); 1093 if (b1 == BIGINT_INVALID) { 1094 Bfree(b); 1095 return b1; 1096 } 1097 x1 = b1->x; 1098 for(i = 0; i < n; i++) 1099 *x1++ = 0; 1100 x = b->x; 1101 xe = x + b->wds; 1102 #ifdef Pack_32 1103 if (k &= 0x1f) { 1104 k1 = 32 - k; 1105 z = 0; 1106 do { 1107 *x1++ = *x << k | z; 1108 z = *x++ >> k1; 1109 } 1110 while(x < xe); 1111 if ((*x1 = z) != 0) 1112 ++n1; 1113 } 1114 #else 1115 if (k &= 0xf) { 1116 k1 = 16 - k; 1117 z = 0; 1118 do { 1119 *x1++ = *x << k & 0xffff | z; 1120 z = *x++ >> k1; 1121 } 1122 while(x < xe); 1123 if (*x1 = z) 1124 ++n1; 1125 } 1126 #endif 1127 else do 1128 *x1++ = *x++; 1129 while(x < xe); 1130 b1->wds = n1 - 1; 1131 Bfree(b); 1132 return b1; 1133 } 1134 1135 static int 1136 cmp 1137 #ifdef KR_headers 1138 (a, b) Bigint *a, *b; 1139 #else 1140 (Bigint *a, Bigint *b) 1141 #endif 1142 { 1143 ULong *xa, *xa0, *xb, *xb0; 1144 int i, j; 1145 1146 if (a == BIGINT_INVALID || b == BIGINT_INVALID) 1147 #ifdef DEBUG 1148 Bug("cmp called with a or b invalid"); 1149 #else 1150 return 0; /* equal - the best we can do right now */ 1151 #endif 1152 1153 i = a->wds; 1154 j = b->wds; 1155 #ifdef DEBUG 1156 if (i > 1 && !a->x[i-1]) 1157 Bug("cmp called with a->x[a->wds-1] == 0"); 1158 if (j > 1 && !b->x[j-1]) 1159 Bug("cmp called with b->x[b->wds-1] == 0"); 1160 #endif 1161 if (i -= j) 1162 return i; 1163 xa0 = a->x; 1164 xa = xa0 + j; 1165 xb0 = b->x; 1166 xb = xb0 + j; 1167 for(;;) { 1168 if (*--xa != *--xb) 1169 return *xa < *xb ? -1 : 1; 1170 if (xa <= xa0) 1171 break; 1172 } 1173 return 0; 1174 } 1175 1176 static Bigint * 1177 diff 1178 #ifdef KR_headers 1179 (a, b) Bigint *a, *b; 1180 #else 1181 (Bigint *a, Bigint *b) 1182 #endif 1183 { 1184 Bigint *c; 1185 int i, wa, wb; 1186 Long borrow, y; /* We need signed shifts here. */ 1187 ULong *xa, *xae, *xb, *xbe, *xc; 1188 #ifdef Pack_32 1189 Long z; 1190 #endif 1191 1192 if (a == BIGINT_INVALID || b == BIGINT_INVALID) 1193 return BIGINT_INVALID; 1194 1195 i = cmp(a,b); 1196 if (!i) { 1197 c = Balloc(0); 1198 if (c != BIGINT_INVALID) { 1199 c->wds = 1; 1200 c->x[0] = 0; 1201 } 1202 return c; 1203 } 1204 if (i < 0) { 1205 c = a; 1206 a = b; 1207 b = c; 1208 i = 1; 1209 } 1210 else 1211 i = 0; 1212 c = Balloc(a->k); 1213 if (c == BIGINT_INVALID) 1214 return c; 1215 c->sign = i; 1216 wa = a->wds; 1217 xa = a->x; 1218 xae = xa + wa; 1219 wb = b->wds; 1220 xb = b->x; 1221 xbe = xb + wb; 1222 xc = c->x; 1223 borrow = 0; 1224 #ifdef Pack_32 1225 do { 1226 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow; 1227 borrow = (ULong)y >> 16; 1228 Sign_Extend(borrow, y); 1229 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow; 1230 borrow = (ULong)z >> 16; 1231 Sign_Extend(borrow, z); 1232 Storeinc(xc, z, y); 1233 } 1234 while(xb < xbe); 1235 while(xa < xae) { 1236 y = (*xa & 0xffff) + borrow; 1237 borrow = (ULong)y >> 16; 1238 Sign_Extend(borrow, y); 1239 z = (*xa++ >> 16) + borrow; 1240 borrow = (ULong)z >> 16; 1241 Sign_Extend(borrow, z); 1242 Storeinc(xc, z, y); 1243 } 1244 #else 1245 do { 1246 y = *xa++ - *xb++ + borrow; 1247 borrow = y >> 16; 1248 Sign_Extend(borrow, y); 1249 *xc++ = y & 0xffff; 1250 } 1251 while(xb < xbe); 1252 while(xa < xae) { 1253 y = *xa++ + borrow; 1254 borrow = y >> 16; 1255 Sign_Extend(borrow, y); 1256 *xc++ = y & 0xffff; 1257 } 1258 #endif 1259 while(!*--xc) 1260 wa--; 1261 c->wds = wa; 1262 return c; 1263 } 1264 1265 static double 1266 ulp 1267 #ifdef KR_headers 1268 (_x) double _x; 1269 #else 1270 (double _x) 1271 #endif 1272 { 1273 _double x; 1274 Long L; 1275 _double a; 1276 1277 value(x) = _x; 1278 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 1279 #ifndef Sudden_Underflow 1280 if (L > 0) { 1281 #endif 1282 #ifdef IBM 1283 L |= Exp_msk1 >> 4; 1284 #endif 1285 word0(a) = L; 1286 word1(a) = 0; 1287 #ifndef Sudden_Underflow 1288 } 1289 else { 1290 L = (ULong)-L >> Exp_shift; 1291 if (L < Exp_shift) { 1292 word0(a) = 0x80000 >> L; 1293 word1(a) = 0; 1294 } 1295 else { 1296 word0(a) = 0; 1297 L -= Exp_shift; 1298 word1(a) = L >= 31 ? 1 : 1 << (31 - L); 1299 } 1300 } 1301 #endif 1302 return value(a); 1303 } 1304 1305 static double 1306 b2d 1307 #ifdef KR_headers 1308 (a, e) Bigint *a; int *e; 1309 #else 1310 (Bigint *a, int *e) 1311 #endif 1312 { 1313 ULong *xa, *xa0, w, y, z; 1314 int k; 1315 _double d; 1316 #ifdef VAX 1317 ULong d0, d1; 1318 #else 1319 #define d0 word0(d) 1320 #define d1 word1(d) 1321 #endif 1322 1323 if (a == BIGINT_INVALID) 1324 return NAN; 1325 1326 xa0 = a->x; 1327 xa = xa0 + a->wds; 1328 y = *--xa; 1329 #ifdef DEBUG 1330 if (!y) Bug("zero y in b2d"); 1331 #endif 1332 k = hi0bits(y); 1333 *e = 32 - k; 1334 #ifdef Pack_32 1335 if (k < Ebits) { 1336 d0 = Exp_1 | y >> (Ebits - k); 1337 w = xa > xa0 ? *--xa : 0; 1338 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k); 1339 goto ret_d; 1340 } 1341 z = xa > xa0 ? *--xa : 0; 1342 if (k -= Ebits) { 1343 d0 = Exp_1 | y << k | z >> (32 - k); 1344 y = xa > xa0 ? *--xa : 0; 1345 d1 = z << k | y >> (32 - k); 1346 } 1347 else { 1348 d0 = Exp_1 | y; 1349 d1 = z; 1350 } 1351 #else 1352 if (k < Ebits + 16) { 1353 z = xa > xa0 ? *--xa : 0; 1354 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 1355 w = xa > xa0 ? *--xa : 0; 1356 y = xa > xa0 ? *--xa : 0; 1357 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 1358 goto ret_d; 1359 } 1360 z = xa > xa0 ? *--xa : 0; 1361 w = xa > xa0 ? *--xa : 0; 1362 k -= Ebits + 16; 1363 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 1364 y = xa > xa0 ? *--xa : 0; 1365 d1 = w << k + 16 | y << k; 1366 #endif 1367 ret_d: 1368 #ifdef VAX 1369 word0(d) = d0 >> 16 | d0 << 16; 1370 word1(d) = d1 >> 16 | d1 << 16; 1371 #else 1372 #undef d0 1373 #undef d1 1374 #endif 1375 return value(d); 1376 } 1377 1378 static Bigint * 1379 d2b 1380 #ifdef KR_headers 1381 (_d, e, bits) double d; int *e, *bits; 1382 #else 1383 (double _d, int *e, int *bits) 1384 #endif 1385 { 1386 Bigint *b; 1387 int de, i, k; 1388 ULong *x, y, z; 1389 _double d; 1390 #ifdef VAX 1391 ULong d0, d1; 1392 #endif 1393 1394 value(d) = _d; 1395 #ifdef VAX 1396 d0 = word0(d) >> 16 | word0(d) << 16; 1397 d1 = word1(d) >> 16 | word1(d) << 16; 1398 #else 1399 #define d0 word0(d) 1400 #define d1 word1(d) 1401 #endif 1402 1403 #ifdef Pack_32 1404 b = Balloc(1); 1405 #else 1406 b = Balloc(2); 1407 #endif 1408 if (b == BIGINT_INVALID) 1409 return b; 1410 x = b->x; 1411 1412 z = d0 & Frac_mask; 1413 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 1414 #ifdef Sudden_Underflow 1415 de = (int)(d0 >> Exp_shift); 1416 #ifndef IBM 1417 z |= Exp_msk11; 1418 #endif 1419 #else 1420 if ((de = (int)(d0 >> Exp_shift)) != 0) 1421 z |= Exp_msk1; 1422 #endif 1423 #ifdef Pack_32 1424 if ((y = d1) != 0) { 1425 if ((k = lo0bits(&y)) != 0) { 1426 x[0] = y | z << (32 - k); 1427 z >>= k; 1428 } 1429 else 1430 x[0] = y; 1431 i = b->wds = (x[1] = z) ? 2 : 1; 1432 } 1433 else { 1434 #ifdef DEBUG 1435 if (!z) 1436 Bug("Zero passed to d2b"); 1437 #endif 1438 k = lo0bits(&z); 1439 x[0] = z; 1440 i = b->wds = 1; 1441 k += 32; 1442 } 1443 #else 1444 if (y = d1) { 1445 if (k = lo0bits(&y)) 1446 if (k >= 16) { 1447 x[0] = y | z << 32 - k & 0xffff; 1448 x[1] = z >> k - 16 & 0xffff; 1449 x[2] = z >> k; 1450 i = 2; 1451 } 1452 else { 1453 x[0] = y & 0xffff; 1454 x[1] = y >> 16 | z << 16 - k & 0xffff; 1455 x[2] = z >> k & 0xffff; 1456 x[3] = z >> k+16; 1457 i = 3; 1458 } 1459 else { 1460 x[0] = y & 0xffff; 1461 x[1] = y >> 16; 1462 x[2] = z & 0xffff; 1463 x[3] = z >> 16; 1464 i = 3; 1465 } 1466 } 1467 else { 1468 #ifdef DEBUG 1469 if (!z) 1470 Bug("Zero passed to d2b"); 1471 #endif 1472 k = lo0bits(&z); 1473 if (k >= 16) { 1474 x[0] = z; 1475 i = 0; 1476 } 1477 else { 1478 x[0] = z & 0xffff; 1479 x[1] = z >> 16; 1480 i = 1; 1481 } 1482 k += 32; 1483 } 1484 while(!x[i]) 1485 --i; 1486 b->wds = i + 1; 1487 #endif 1488 #ifndef Sudden_Underflow 1489 if (de) { 1490 #endif 1491 #ifdef IBM 1492 *e = (de - Bias - (P-1) << 2) + k; 1493 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 1494 #else 1495 *e = de - Bias - (P-1) + k; 1496 *bits = P - k; 1497 #endif 1498 #ifndef Sudden_Underflow 1499 } 1500 else { 1501 *e = de - Bias - (P-1) + 1 + k; 1502 #ifdef Pack_32 1503 *bits = 32*i - hi0bits(x[i-1]); 1504 #else 1505 *bits = (i+2)*16 - hi0bits(x[i]); 1506 #endif 1507 } 1508 #endif 1509 return b; 1510 } 1511 #undef d0 1512 #undef d1 1513 1514 static double 1515 ratio 1516 #ifdef KR_headers 1517 (a, b) Bigint *a, *b; 1518 #else 1519 (Bigint *a, Bigint *b) 1520 #endif 1521 { 1522 _double da, db; 1523 int k, ka, kb; 1524 1525 if (a == BIGINT_INVALID || b == BIGINT_INVALID) 1526 return NAN; /* for lack of better value ? */ 1527 1528 value(da) = b2d(a, &ka); 1529 value(db) = b2d(b, &kb); 1530 #ifdef Pack_32 1531 k = ka - kb + 32*(a->wds - b->wds); 1532 #else 1533 k = ka - kb + 16*(a->wds - b->wds); 1534 #endif 1535 #ifdef IBM 1536 if (k > 0) { 1537 word0(da) += (k >> 2)*Exp_msk1; 1538 if (k &= 3) 1539 da *= 1 << k; 1540 } 1541 else { 1542 k = -k; 1543 word0(db) += (k >> 2)*Exp_msk1; 1544 if (k &= 3) 1545 db *= 1 << k; 1546 } 1547 #else 1548 if (k > 0) 1549 word0(da) += k*Exp_msk1; 1550 else { 1551 k = -k; 1552 word0(db) += k*Exp_msk1; 1553 } 1554 #endif 1555 return value(da) / value(db); 1556 } 1557 1558 static CONST double 1559 tens[] = { 1560 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1561 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1562 1e20, 1e21, 1e22 1563 #ifdef VAX 1564 , 1e23, 1e24 1565 #endif 1566 }; 1567 1568 #ifdef IEEE_Arith 1569 static CONST double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1570 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 }; 1571 #define n_bigtens 5 1572 #else 1573 #ifdef IBM 1574 static CONST double bigtens[] = { 1e16, 1e32, 1e64 }; 1575 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 1576 #define n_bigtens 3 1577 #else 1578 static CONST double bigtens[] = { 1e16, 1e32 }; 1579 static CONST double tinytens[] = { 1e-16, 1e-32 }; 1580 #define n_bigtens 2 1581 #endif 1582 #endif 1583 1584 double 1585 strtod 1586 #ifdef KR_headers 1587 (s00, se) CONST char *s00; char **se; 1588 #else 1589 (CONST char *s00, char **se) 1590 #endif 1591 { 1592 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 1593 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 1594 CONST char *s, *s0, *s1; 1595 double aadj, aadj1, adj; 1596 _double rv, rv0; 1597 Long L; 1598 ULong y, z; 1599 Bigint *bb1, *bd0; 1600 Bigint *bb = NULL, *bd = NULL, *bs = NULL, *delta = NULL;/* pacify gcc */ 1601 1602 #ifdef ANDROID_CHANGES 1603 CONST char decimal_point = '.'; 1604 #else /* ANDROID_CHANGES */ 1605 #ifndef KR_headers 1606 CONST char decimal_point = localeconv()->decimal_point[0]; 1607 #else 1608 CONST char decimal_point = '.'; 1609 #endif 1610 1611 #endif /* ANDROID_CHANGES */ 1612 1613 sign = nz0 = nz = 0; 1614 value(rv) = 0.; 1615 1616 1617 for(s = s00; isspace((unsigned char) *s); s++) 1618 ; 1619 1620 if (*s == '-') { 1621 sign = 1; 1622 s++; 1623 } else if (*s == '+') { 1624 s++; 1625 } 1626 1627 if (*s == '\0') { 1628 s = s00; 1629 goto ret; 1630 } 1631 1632 /* "INF" or "INFINITY" */ 1633 if (tolower((unsigned char)*s) == 'i' && strncasecmp(s, "inf", 3) == 0) { 1634 if (strncasecmp(s + 3, "inity", 5) == 0) 1635 s += 8; 1636 else 1637 s += 3; 1638 1639 value(rv) = HUGE_VAL; 1640 goto ret; 1641 } 1642 1643 #ifdef IEEE_Arith 1644 /* "NAN" or "NAN(n-char-sequence-opt)" */ 1645 if (tolower((unsigned char)*s) == 'n' && strncasecmp(s, "nan", 3) == 0) { 1646 /* Build a quiet NaN. */ 1647 word0(rv) = NAN_WORD0; 1648 word1(rv) = NAN_WORD1; 1649 s+= 3; 1650 1651 /* Don't interpret (n-char-sequence-opt), for now. */ 1652 if (*s == '(') { 1653 s0 = s; 1654 for (s++; *s != ')' && *s != '\0'; s++) 1655 ; 1656 if (*s == ')') 1657 s++; /* Skip over closing paren ... */ 1658 else 1659 s = s0; /* ... otherwise go back. */ 1660 } 1661 1662 goto ret; 1663 } 1664 #endif 1665 1666 if (*s == '0') { 1667 #ifndef NO_HEX_FP /*{*/ 1668 { 1669 static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 1670 Long exp; 1671 ULong bits[2]; 1672 switch(s[1]) { 1673 case 'x': 1674 case 'X': 1675 { 1676 #ifdef Honor_FLT_ROUNDS 1677 FPI fpi1 = fpi; 1678 fpi1.rounding = Rounding; 1679 #else 1680 #define fpi1 fpi 1681 #endif 1682 switch((i = gethex(&s, &fpi1, &exp, &bb, sign, 0)) & STRTOG_Retmask) { 1683 case STRTOG_NoNumber: 1684 s = s00; 1685 sign = 0; 1686 case STRTOG_Zero: 1687 break; 1688 default: 1689 if (bb) { 1690 copybits(bits, fpi.nbits, bb); 1691 Bfree(bb); 1692 } 1693 ULtod(((U*)&rv)->L, bits, exp, i); 1694 }} 1695 goto ret; 1696 } 1697 } 1698 #endif /*}*/ 1699 nz0 = 1; 1700 while(*++s == '0') ; 1701 if (!*s) 1702 goto ret; 1703 } 1704 s0 = s; 1705 y = z = 0; 1706 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 1707 if (nd < 9) 1708 y = 10*y + c - '0'; 1709 else if (nd < 16) 1710 z = 10*z + c - '0'; 1711 nd0 = nd; 1712 if (c == decimal_point) { 1713 c = *++s; 1714 if (!nd) { 1715 for(; c == '0'; c = *++s) 1716 nz++; 1717 if (c > '0' && c <= '9') { 1718 s0 = s; 1719 nf += nz; 1720 nz = 0; 1721 goto have_dig; 1722 } 1723 goto dig_done; 1724 } 1725 for(; c >= '0' && c <= '9'; c = *++s) { 1726 have_dig: 1727 nz++; 1728 if (c -= '0') { 1729 nf += nz; 1730 for(i = 1; i < nz; i++) 1731 if (nd++ < 9) 1732 y *= 10; 1733 else if (nd <= DBL_DIG + 1) 1734 z *= 10; 1735 if (nd++ < 9) 1736 y = 10*y + c; 1737 else if (nd <= DBL_DIG + 1) 1738 z = 10*z + c; 1739 nz = 0; 1740 } 1741 } 1742 } 1743 dig_done: 1744 e = 0; 1745 if (c == 'e' || c == 'E') { 1746 if (!nd && !nz && !nz0) { 1747 s = s00; 1748 goto ret; 1749 } 1750 s00 = s; 1751 esign = 0; 1752 switch(c = *++s) { 1753 case '-': 1754 esign = 1; 1755 /* FALLTHROUGH */ 1756 case '+': 1757 c = *++s; 1758 } 1759 if (c >= '0' && c <= '9') { 1760 while(c == '0') 1761 c = *++s; 1762 if (c > '0' && c <= '9') { 1763 L = c - '0'; 1764 s1 = s; 1765 while((c = *++s) >= '0' && c <= '9') 1766 L = 10*L + c - '0'; 1767 if (s - s1 > 8 || L > 19999) 1768 /* Avoid confusion from exponents 1769 * so large that e might overflow. 1770 */ 1771 e = 19999; /* safe for 16 bit ints */ 1772 else 1773 e = (int)L; 1774 if (esign) 1775 e = -e; 1776 } 1777 else 1778 e = 0; 1779 } 1780 else 1781 s = s00; 1782 } 1783 if (!nd) { 1784 if (!nz && !nz0) 1785 s = s00; 1786 goto ret; 1787 } 1788 e1 = e -= nf; 1789 1790 /* Now we have nd0 digits, starting at s0, followed by a 1791 * decimal point, followed by nd-nd0 digits. The number we're 1792 * after is the integer represented by those digits times 1793 * 10**e */ 1794 1795 if (!nd0) 1796 nd0 = nd; 1797 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 1798 value(rv) = y; 1799 if (k > 9) 1800 value(rv) = tens[k - 9] * value(rv) + z; 1801 bd0 = 0; 1802 if (nd <= DBL_DIG 1803 #ifndef RND_PRODQUOT 1804 && FLT_ROUNDS == 1 1805 #endif 1806 ) { 1807 if (!e) 1808 goto ret; 1809 if (e > 0) { 1810 if (e <= Ten_pmax) { 1811 #ifdef VAX 1812 goto vax_ovfl_check; 1813 #else 1814 /* value(rv) = */ rounded_product(value(rv), 1815 tens[e]); 1816 goto ret; 1817 #endif 1818 } 1819 i = DBL_DIG - nd; 1820 if (e <= Ten_pmax + i) { 1821 /* A fancier test would sometimes let us do 1822 * this for larger i values. 1823 */ 1824 e -= i; 1825 value(rv) *= tens[i]; 1826 #ifdef VAX 1827 /* VAX exponent range is so narrow we must 1828 * worry about overflow here... 1829 */ 1830 vax_ovfl_check: 1831 word0(rv) -= P*Exp_msk1; 1832 /* value(rv) = */ rounded_product(value(rv), 1833 tens[e]); 1834 if ((word0(rv) & Exp_mask) 1835 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 1836 goto ovfl; 1837 word0(rv) += P*Exp_msk1; 1838 #else 1839 /* value(rv) = */ rounded_product(value(rv), 1840 tens[e]); 1841 #endif 1842 goto ret; 1843 } 1844 } 1845 #ifndef Inaccurate_Divide 1846 else if (e >= -Ten_pmax) { 1847 /* value(rv) = */ rounded_quotient(value(rv), 1848 tens[-e]); 1849 goto ret; 1850 } 1851 #endif 1852 } 1853 e1 += nd - k; 1854 1855 /* Get starting approximation = rv * 10**e1 */ 1856 1857 if (e1 > 0) { 1858 if ((i = e1 & 15) != 0) 1859 value(rv) *= tens[i]; 1860 if (e1 &= ~15) { 1861 if (e1 > DBL_MAX_10_EXP) { 1862 ovfl: 1863 errno = ERANGE; 1864 value(rv) = HUGE_VAL; 1865 if (bd0) 1866 goto retfree; 1867 goto ret; 1868 } 1869 if ((e1 = (unsigned int)e1 >> 4) != 0) { 1870 for(j = 0; e1 > 1; j++, 1871 e1 = (unsigned int)e1 >> 1) 1872 if (e1 & 1) 1873 value(rv) *= bigtens[j]; 1874 /* The last multiplication could overflow. */ 1875 word0(rv) -= P*Exp_msk1; 1876 value(rv) *= bigtens[j]; 1877 if ((z = word0(rv) & Exp_mask) 1878 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 1879 goto ovfl; 1880 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 1881 /* set to largest number */ 1882 /* (Can't trust DBL_MAX) */ 1883 word0(rv) = Big0; 1884 word1(rv) = Big1; 1885 } 1886 else 1887 word0(rv) += P*Exp_msk1; 1888 } 1889 } 1890 } 1891 else if (e1 < 0) { 1892 e1 = -e1; 1893 if ((i = e1 & 15) != 0) 1894 value(rv) /= tens[i]; 1895 if (e1 &= ~15) { 1896 e1 = (unsigned int)e1 >> 4; 1897 if (e1 >= 1 << n_bigtens) 1898 goto undfl; 1899 for(j = 0; e1 > 1; j++, 1900 e1 = (unsigned int)e1 >> 1) 1901 if (e1 & 1) 1902 value(rv) *= tinytens[j]; 1903 /* The last multiplication could underflow. */ 1904 value(rv0) = value(rv); 1905 value(rv) *= tinytens[j]; 1906 if (!value(rv)) { 1907 value(rv) = 2.*value(rv0); 1908 value(rv) *= tinytens[j]; 1909 if (!value(rv)) { 1910 undfl: 1911 value(rv) = 0.; 1912 errno = ERANGE; 1913 if (bd0) 1914 goto retfree; 1915 goto ret; 1916 } 1917 word0(rv) = Tiny0; 1918 word1(rv) = Tiny1; 1919 /* The refinement below will clean 1920 * this approximation up. 1921 */ 1922 } 1923 } 1924 } 1925 1926 /* Now the hard part -- adjusting rv to the correct value.*/ 1927 1928 /* Put digits into bd: true value = bd * 10^e */ 1929 1930 bd0 = s2b(s0, nd0, nd, y); 1931 1932 for(;;) { 1933 bd = Balloc(bd0->k); 1934 Bcopy(bd, bd0); 1935 bb = d2b(value(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 1936 bs = i2b(1); 1937 1938 if (e >= 0) { 1939 bb2 = bb5 = 0; 1940 bd2 = bd5 = e; 1941 } 1942 else { 1943 bb2 = bb5 = -e; 1944 bd2 = bd5 = 0; 1945 } 1946 if (bbe >= 0) 1947 bb2 += bbe; 1948 else 1949 bd2 -= bbe; 1950 bs2 = bb2; 1951 #ifdef Sudden_Underflow 1952 #ifdef IBM 1953 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 1954 #else 1955 j = P + 1 - bbbits; 1956 #endif 1957 #else 1958 i = bbe + bbbits - 1; /* logb(rv) */ 1959 if (i < Emin) /* denormal */ 1960 j = bbe + (P-Emin); 1961 else 1962 j = P + 1 - bbbits; 1963 #endif 1964 bb2 += j; 1965 bd2 += j; 1966 i = bb2 < bd2 ? bb2 : bd2; 1967 if (i > bs2) 1968 i = bs2; 1969 if (i > 0) { 1970 bb2 -= i; 1971 bd2 -= i; 1972 bs2 -= i; 1973 } 1974 if (bb5 > 0) { 1975 bs = pow5mult(bs, bb5); 1976 bb1 = mult(bs, bb); 1977 Bfree(bb); 1978 bb = bb1; 1979 } 1980 if (bb2 > 0) 1981 bb = lshift(bb, bb2); 1982 if (bd5 > 0) 1983 bd = pow5mult(bd, bd5); 1984 if (bd2 > 0) 1985 bd = lshift(bd, bd2); 1986 if (bs2 > 0) 1987 bs = lshift(bs, bs2); 1988 delta = diff(bb, bd); 1989 dsign = delta->sign; 1990 delta->sign = 0; 1991 i = cmp(delta, bs); 1992 if (i < 0) { 1993 /* Error is less than half an ulp -- check for 1994 * special case of mantissa a power of two. 1995 */ 1996 if (dsign || word1(rv) || word0(rv) & Bndry_mask) 1997 break; 1998 delta = lshift(delta,Log2P); 1999 if (cmp(delta, bs) > 0) 2000 goto drop_down; 2001 break; 2002 } 2003 if (i == 0) { 2004 /* exactly half-way between */ 2005 if (dsign) { 2006 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 2007 && word1(rv) == 0xffffffff) { 2008 /*boundary case -- increment exponent*/ 2009 word0(rv) = (word0(rv) & Exp_mask) 2010 + Exp_msk1 2011 #ifdef IBM 2012 | Exp_msk1 >> 4 2013 #endif 2014 ; 2015 word1(rv) = 0; 2016 break; 2017 } 2018 } 2019 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 2020 drop_down: 2021 /* boundary case -- decrement exponent */ 2022 #ifdef Sudden_Underflow 2023 L = word0(rv) & Exp_mask; 2024 #ifdef IBM 2025 if (L < Exp_msk1) 2026 #else 2027 if (L <= Exp_msk1) 2028 #endif 2029 goto undfl; 2030 L -= Exp_msk1; 2031 #else 2032 L = (word0(rv) & Exp_mask) - Exp_msk1; 2033 #endif 2034 word0(rv) = L | Bndry_mask1; 2035 word1(rv) = 0xffffffff; 2036 #ifdef IBM 2037 goto cont; 2038 #else 2039 break; 2040 #endif 2041 } 2042 #ifndef ROUND_BIASED 2043 if (!(word1(rv) & LSB)) 2044 break; 2045 #endif 2046 if (dsign) 2047 value(rv) += ulp(value(rv)); 2048 #ifndef ROUND_BIASED 2049 else { 2050 value(rv) -= ulp(value(rv)); 2051 #ifndef Sudden_Underflow 2052 if (!value(rv)) 2053 goto undfl; 2054 #endif 2055 } 2056 #endif 2057 break; 2058 } 2059 if ((aadj = ratio(delta, bs)) <= 2.) { 2060 if (dsign) 2061 aadj = aadj1 = 1.; 2062 else if (word1(rv) || word0(rv) & Bndry_mask) { 2063 #ifndef Sudden_Underflow 2064 if (word1(rv) == Tiny1 && !word0(rv)) 2065 goto undfl; 2066 #endif 2067 aadj = 1.; 2068 aadj1 = -1.; 2069 } 2070 else { 2071 /* special case -- power of FLT_RADIX to be */ 2072 /* rounded down... */ 2073 2074 if (aadj < 2./FLT_RADIX) 2075 aadj = 1./FLT_RADIX; 2076 else 2077 aadj *= 0.5; 2078 aadj1 = -aadj; 2079 } 2080 } 2081 else { 2082 aadj *= 0.5; 2083 aadj1 = dsign ? aadj : -aadj; 2084 #ifdef Check_FLT_ROUNDS 2085 switch(FLT_ROUNDS) { 2086 case 2: /* towards +infinity */ 2087 aadj1 -= 0.5; 2088 break; 2089 case 0: /* towards 0 */ 2090 case 3: /* towards -infinity */ 2091 aadj1 += 0.5; 2092 } 2093 #else 2094 if (FLT_ROUNDS == 0) 2095 aadj1 += 0.5; 2096 #endif 2097 } 2098 y = word0(rv) & Exp_mask; 2099 2100 /* Check for overflow */ 2101 2102 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 2103 value(rv0) = value(rv); 2104 word0(rv) -= P*Exp_msk1; 2105 adj = aadj1 * ulp(value(rv)); 2106 value(rv) += adj; 2107 if ((word0(rv) & Exp_mask) >= 2108 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 2109 if (word0(rv0) == Big0 && word1(rv0) == Big1) 2110 goto ovfl; 2111 word0(rv) = Big0; 2112 word1(rv) = Big1; 2113 goto cont; 2114 } 2115 else 2116 word0(rv) += P*Exp_msk1; 2117 } 2118 else { 2119 #ifdef Sudden_Underflow 2120 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 2121 value(rv0) = value(rv); 2122 word0(rv) += P*Exp_msk1; 2123 adj = aadj1 * ulp(value(rv)); 2124 value(rv) += adj; 2125 #ifdef IBM 2126 if ((word0(rv) & Exp_mask) < P*Exp_msk1) 2127 #else 2128 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 2129 #endif 2130 { 2131 if (word0(rv0) == Tiny0 2132 && word1(rv0) == Tiny1) 2133 goto undfl; 2134 word0(rv) = Tiny0; 2135 word1(rv) = Tiny1; 2136 goto cont; 2137 } 2138 else 2139 word0(rv) -= P*Exp_msk1; 2140 } 2141 else { 2142 adj = aadj1 * ulp(value(rv)); 2143 value(rv) += adj; 2144 } 2145 #else 2146 /* Compute adj so that the IEEE rounding rules will 2147 * correctly round rv + adj in some half-way cases. 2148 * If rv * ulp(rv) is denormalized (i.e., 2149 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 2150 * trouble from bits lost to denormalization; 2151 * example: 1.2e-307 . 2152 */ 2153 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) { 2154 aadj1 = (double)(int)(aadj + 0.5); 2155 if (!dsign) 2156 aadj1 = -aadj1; 2157 } 2158 adj = aadj1 * ulp(value(rv)); 2159 value(rv) += adj; 2160 #endif 2161 } 2162 z = word0(rv) & Exp_mask; 2163 if (y == z) { 2164 /* Can we stop now? */ 2165 L = aadj; 2166 aadj -= L; 2167 /* The tolerances below are conservative. */ 2168 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 2169 if (aadj < .4999999 || aadj > .5000001) 2170 break; 2171 } 2172 else if (aadj < .4999999/FLT_RADIX) 2173 break; 2174 } 2175 cont: 2176 Bfree(bb); 2177 Bfree(bd); 2178 Bfree(bs); 2179 Bfree(delta); 2180 } 2181 retfree: 2182 Bfree(bb); 2183 Bfree(bd); 2184 Bfree(bs); 2185 Bfree(bd0); 2186 Bfree(delta); 2187 ret: 2188 if (se) 2189 /* LINTED interface specification */ 2190 *se = (char *)s; 2191 return sign ? -value(rv) : value(rv); 2192 } 2193 2194 static int 2195 quorem 2196 #ifdef KR_headers 2197 (b, S) Bigint *b, *S; 2198 #else 2199 (Bigint *b, Bigint *S) 2200 #endif 2201 { 2202 int n; 2203 Long borrow, y; 2204 ULong carry, q, ys; 2205 ULong *bx, *bxe, *sx, *sxe; 2206 #ifdef Pack_32 2207 Long z; 2208 ULong si, zs; 2209 #endif 2210 2211 if (b == BIGINT_INVALID || S == BIGINT_INVALID) 2212 return 0; 2213 2214 n = S->wds; 2215 #ifdef DEBUG 2216 /*debug*/ if (b->wds > n) 2217 /*debug*/ Bug("oversize b in quorem"); 2218 #endif 2219 if (b->wds < n) 2220 return 0; 2221 sx = S->x; 2222 sxe = sx + --n; 2223 bx = b->x; 2224 bxe = bx + n; 2225 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 2226 #ifdef DEBUG 2227 /*debug*/ if (q > 9) 2228 /*debug*/ Bug("oversized quotient in quorem"); 2229 #endif 2230 if (q) { 2231 borrow = 0; 2232 carry = 0; 2233 do { 2234 #ifdef Pack_32 2235 si = *sx++; 2236 ys = (si & 0xffff) * q + carry; 2237 zs = (si >> 16) * q + (ys >> 16); 2238 carry = zs >> 16; 2239 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 2240 borrow = (ULong)y >> 16; 2241 Sign_Extend(borrow, y); 2242 z = (*bx >> 16) - (zs & 0xffff) + borrow; 2243 borrow = (ULong)z >> 16; 2244 Sign_Extend(borrow, z); 2245 Storeinc(bx, z, y); 2246 #else 2247 ys = *sx++ * q + carry; 2248 carry = ys >> 16; 2249 y = *bx - (ys & 0xffff) + borrow; 2250 borrow = y >> 16; 2251 Sign_Extend(borrow, y); 2252 *bx++ = y & 0xffff; 2253 #endif 2254 } 2255 while(sx <= sxe); 2256 if (!*bxe) { 2257 bx = b->x; 2258 while(--bxe > bx && !*bxe) 2259 --n; 2260 b->wds = n; 2261 } 2262 } 2263 if (cmp(b, S) >= 0) { 2264 q++; 2265 borrow = 0; 2266 carry = 0; 2267 bx = b->x; 2268 sx = S->x; 2269 do { 2270 #ifdef Pack_32 2271 si = *sx++; 2272 ys = (si & 0xffff) + carry; 2273 zs = (si >> 16) + (ys >> 16); 2274 carry = zs >> 16; 2275 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 2276 borrow = (ULong)y >> 16; 2277 Sign_Extend(borrow, y); 2278 z = (*bx >> 16) - (zs & 0xffff) + borrow; 2279 borrow = (ULong)z >> 16; 2280 Sign_Extend(borrow, z); 2281 Storeinc(bx, z, y); 2282 #else 2283 ys = *sx++ + carry; 2284 carry = ys >> 16; 2285 y = *bx - (ys & 0xffff) + borrow; 2286 borrow = y >> 16; 2287 Sign_Extend(borrow, y); 2288 *bx++ = y & 0xffff; 2289 #endif 2290 } 2291 while(sx <= sxe); 2292 bx = b->x; 2293 bxe = bx + n; 2294 if (!*bxe) { 2295 while(--bxe > bx && !*bxe) 2296 --n; 2297 b->wds = n; 2298 } 2299 } 2300 return q; 2301 } 2302 2303 /* freedtoa(s) must be used to free values s returned by dtoa 2304 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 2305 * but for consistency with earlier versions of dtoa, it is optional 2306 * when MULTIPLE_THREADS is not defined. 2307 */ 2308 2309 void 2310 #ifdef KR_headers 2311 freedtoa(s) char *s; 2312 #else 2313 freedtoa(char *s) 2314 #endif 2315 { 2316 free(s); 2317 } 2318 2319 2320 2321 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 2322 * 2323 * Inspired by "How to Print Floating-Point Numbers Accurately" by 2324 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101]. 2325 * 2326 * Modifications: 2327 * 1. Rather than iterating, we use a simple numeric overestimate 2328 * to determine k = floor(log10(d)). We scale relevant 2329 * quantities using O(log2(k)) rather than O(k) multiplications. 2330 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 2331 * try to generate digits strictly left to right. Instead, we 2332 * compute with fewer bits and propagate the carry if necessary 2333 * when rounding the final digit up. This is often faster. 2334 * 3. Under the assumption that input will be rounded nearest, 2335 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 2336 * That is, we allow equality in stopping tests when the 2337 * round-nearest rule will give the same floating-point value 2338 * as would satisfaction of the stopping test with strict 2339 * inequality. 2340 * 4. We remove common factors of powers of 2 from relevant 2341 * quantities. 2342 * 5. When converting floating-point integers less than 1e16, 2343 * we use floating-point arithmetic rather than resorting 2344 * to multiple-precision integers. 2345 * 6. When asked to produce fewer than 15 digits, we first try 2346 * to get by with floating-point arithmetic; we resort to 2347 * multiple-precision integer arithmetic only if we cannot 2348 * guarantee that the floating-point calculation has given 2349 * the correctly rounded result. For k requested digits and 2350 * "uniformly" distributed input, the probability is 2351 * something like 10^(k-15) that we must resort to the Long 2352 * calculation. 2353 */ 2354 2355 __LIBC_HIDDEN__ char * 2356 __dtoa 2357 #ifdef KR_headers 2358 (_d, mode, ndigits, decpt, sign, rve) 2359 double _d; int mode, ndigits, *decpt, *sign; char **rve; 2360 #else 2361 (double _d, int mode, int ndigits, int *decpt, int *sign, char **rve) 2362 #endif 2363 { 2364 /* Arguments ndigits, decpt, sign are similar to those 2365 of ecvt and fcvt; trailing zeros are suppressed from 2366 the returned string. If not null, *rve is set to point 2367 to the end of the return value. If d is +-Infinity or NaN, 2368 then *decpt is set to 9999. 2369 2370 mode: 2371 0 ==> shortest string that yields d when read in 2372 and rounded to nearest. 2373 1 ==> like 0, but with Steele & White stopping rule; 2374 e.g. with IEEE P754 arithmetic , mode 0 gives 2375 1e23 whereas mode 1 gives 9.999999999999999e22. 2376 2 ==> max(1,ndigits) significant digits. This gives a 2377 return value similar to that of ecvt, except 2378 that trailing zeros are suppressed. 2379 3 ==> through ndigits past the decimal point. This 2380 gives a return value similar to that from fcvt, 2381 except that trailing zeros are suppressed, and 2382 ndigits can be negative. 2383 4-9 should give the same return values as 2-3, i.e., 2384 4 <= mode <= 9 ==> same return as mode 2385 2 + (mode & 1). These modes are mainly for 2386 debugging; often they run slower but sometimes 2387 faster than modes 2-3. 2388 4,5,8,9 ==> left-to-right digit generation. 2389 6-9 ==> don't try fast floating-point estimate 2390 (if applicable). 2391 2392 Values of mode other than 0-9 are treated as mode 0. 2393 2394 Sufficient space is allocated to the return value 2395 to hold the suppressed trailing zeros. 2396 */ 2397 2398 int bbits, b2, b5, be, dig, i, ieps, ilim0, 2399 j, jj1, k, k0, k_check, leftright, m2, m5, s2, s5, 2400 try_quick; 2401 int ilim = 0, ilim1 = 0, spec_case = 0; /* pacify gcc */ 2402 Long L; 2403 #ifndef Sudden_Underflow 2404 int denorm; 2405 ULong x; 2406 #endif 2407 Bigint *b, *b1, *delta, *mhi, *S; 2408 Bigint *mlo = NULL; /* pacify gcc */ 2409 double ds; 2410 char *s, *s0; 2411 Bigint *result = NULL; 2412 int result_k = 0; 2413 _double d, d2, eps; 2414 2415 value(d) = _d; 2416 2417 if (word0(d) & Sign_bit) { 2418 /* set sign for everything, including 0's and NaNs */ 2419 *sign = 1; 2420 word0(d) &= ~Sign_bit; /* clear sign bit */ 2421 } 2422 else 2423 *sign = 0; 2424 2425 #if defined(IEEE_Arith) + defined(VAX) 2426 #ifdef IEEE_Arith 2427 if ((word0(d) & Exp_mask) == Exp_mask) 2428 #else 2429 if (word0(d) == 0x8000) 2430 #endif 2431 { 2432 /* Infinity or NaN */ 2433 *decpt = 9999; 2434 s = 2435 #ifdef IEEE_Arith 2436 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" : 2437 #endif 2438 "NaN"; 2439 result = Balloc(strlen(s)+1); 2440 if (result == BIGINT_INVALID) 2441 return NULL; 2442 s0 = (char *)(void *)result; 2443 strcpy(s0, s); 2444 if (rve) 2445 *rve = 2446 #ifdef IEEE_Arith 2447 s0[3] ? s0 + 8 : 2448 #endif 2449 s0 + 3; 2450 return s0; 2451 } 2452 #endif 2453 #ifdef IBM 2454 value(d) += 0; /* normalize */ 2455 #endif 2456 if (!value(d)) { 2457 *decpt = 1; 2458 result = Balloc(2); 2459 if (result == BIGINT_INVALID) 2460 return NULL; 2461 s0 = (char *)(void *)result; 2462 strcpy(s0, "0"); 2463 if (rve) 2464 *rve = s0 + 1; 2465 return s0; 2466 } 2467 2468 b = d2b(value(d), &be, &bbits); 2469 #ifdef Sudden_Underflow 2470 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); 2471 #else 2472 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) { 2473 #endif 2474 value(d2) = value(d); 2475 word0(d2) &= Frac_mask1; 2476 word0(d2) |= Exp_11; 2477 #ifdef IBM 2478 if (j = 11 - hi0bits(word0(d2) & Frac_mask)) 2479 value(d2) /= 1 << j; 2480 #endif 2481 2482 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 2483 * log10(x) = log(x) / log(10) 2484 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 2485 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 2486 * 2487 * This suggests computing an approximation k to log10(d) by 2488 * 2489 * k = (i - Bias)*0.301029995663981 2490 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 2491 * 2492 * We want k to be too large rather than too small. 2493 * The error in the first-order Taylor series approximation 2494 * is in our favor, so we just round up the constant enough 2495 * to compensate for any error in the multiplication of 2496 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 2497 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 2498 * adding 1e-13 to the constant term more than suffices. 2499 * Hence we adjust the constant term to 0.1760912590558. 2500 * (We could get a more accurate k by invoking log10, 2501 * but this is probably not worthwhile.) 2502 */ 2503 2504 i -= Bias; 2505 #ifdef IBM 2506 i <<= 2; 2507 i += j; 2508 #endif 2509 #ifndef Sudden_Underflow 2510 denorm = 0; 2511 } 2512 else { 2513 /* d is denormalized */ 2514 2515 i = bbits + be + (Bias + (P-1) - 1); 2516 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32) 2517 : word1(d) << (32 - i); 2518 value(d2) = x; 2519 word0(d2) -= 31*Exp_msk1; /* adjust exponent */ 2520 i -= (Bias + (P-1) - 1) + 1; 2521 denorm = 1; 2522 } 2523 #endif 2524 ds = (value(d2)-1.5)*0.289529654602168 + 0.1760912590558 + 2525 i*0.301029995663981; 2526 k = (int)ds; 2527 if (ds < 0. && ds != k) 2528 k--; /* want k = floor(ds) */ 2529 k_check = 1; 2530 if (k >= 0 && k <= Ten_pmax) { 2531 if (value(d) < tens[k]) 2532 k--; 2533 k_check = 0; 2534 } 2535 j = bbits - i - 1; 2536 if (j >= 0) { 2537 b2 = 0; 2538 s2 = j; 2539 } 2540 else { 2541 b2 = -j; 2542 s2 = 0; 2543 } 2544 if (k >= 0) { 2545 b5 = 0; 2546 s5 = k; 2547 s2 += k; 2548 } 2549 else { 2550 b2 -= k; 2551 b5 = -k; 2552 s5 = 0; 2553 } 2554 if (mode < 0 || mode > 9) 2555 mode = 0; 2556 try_quick = 1; 2557 if (mode > 5) { 2558 mode -= 4; 2559 try_quick = 0; 2560 } 2561 leftright = 1; 2562 switch(mode) { 2563 case 0: 2564 case 1: 2565 ilim = ilim1 = -1; 2566 i = 18; 2567 ndigits = 0; 2568 break; 2569 case 2: 2570 leftright = 0; 2571 /* FALLTHROUGH */ 2572 case 4: 2573 if (ndigits <= 0) 2574 ndigits = 1; 2575 ilim = ilim1 = i = ndigits; 2576 break; 2577 case 3: 2578 leftright = 0; 2579 /* FALLTHROUGH */ 2580 case 5: 2581 i = ndigits + k + 1; 2582 ilim = i; 2583 ilim1 = i - 1; 2584 if (i <= 0) 2585 i = 1; 2586 } 2587 j = sizeof(ULong); 2588 for(result_k = 0; (int)(sizeof(Bigint) - sizeof(ULong)) + j <= i; 2589 j <<= 1) result_k++; 2590 // this is really a ugly hack, the code uses Balloc 2591 // instead of malloc, but casts the result into a char* 2592 // it seems the only reason to do that is due to the 2593 // complicated way the block size need to be computed 2594 // buuurk.... 2595 result = Balloc(result_k); 2596 if (result == BIGINT_INVALID) { 2597 Bfree(b); 2598 return NULL; 2599 } 2600 s = s0 = (char *)(void *)result; 2601 2602 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 2603 2604 /* Try to get by with floating-point arithmetic. */ 2605 2606 i = 0; 2607 value(d2) = value(d); 2608 k0 = k; 2609 ilim0 = ilim; 2610 ieps = 2; /* conservative */ 2611 if (k > 0) { 2612 ds = tens[k&0xf]; 2613 j = (unsigned int)k >> 4; 2614 if (j & Bletch) { 2615 /* prevent overflows */ 2616 j &= Bletch - 1; 2617 value(d) /= bigtens[n_bigtens-1]; 2618 ieps++; 2619 } 2620 for(; j; j = (unsigned int)j >> 1, i++) 2621 if (j & 1) { 2622 ieps++; 2623 ds *= bigtens[i]; 2624 } 2625 value(d) /= ds; 2626 } 2627 else if ((jj1 = -k) != 0) { 2628 value(d) *= tens[jj1 & 0xf]; 2629 for(j = (unsigned int)jj1 >> 4; j; 2630 j = (unsigned int)j >> 1, i++) 2631 if (j & 1) { 2632 ieps++; 2633 value(d) *= bigtens[i]; 2634 } 2635 } 2636 if (k_check && value(d) < 1. && ilim > 0) { 2637 if (ilim1 <= 0) 2638 goto fast_failed; 2639 ilim = ilim1; 2640 k--; 2641 value(d) *= 10.; 2642 ieps++; 2643 } 2644 value(eps) = ieps*value(d) + 7.; 2645 word0(eps) -= (P-1)*Exp_msk1; 2646 if (ilim == 0) { 2647 S = mhi = 0; 2648 value(d) -= 5.; 2649 if (value(d) > value(eps)) 2650 goto one_digit; 2651 if (value(d) < -value(eps)) 2652 goto no_digits; 2653 goto fast_failed; 2654 } 2655 #ifndef No_leftright 2656 if (leftright) { 2657 /* Use Steele & White method of only 2658 * generating digits needed. 2659 */ 2660 value(eps) = 0.5/tens[ilim-1] - value(eps); 2661 for(i = 0;;) { 2662 L = value(d); 2663 value(d) -= L; 2664 *s++ = '0' + (int)L; 2665 if (value(d) < value(eps)) 2666 goto ret1; 2667 if (1. - value(d) < value(eps)) 2668 goto bump_up; 2669 if (++i >= ilim) 2670 break; 2671 value(eps) *= 10.; 2672 value(d) *= 10.; 2673 } 2674 } 2675 else { 2676 #endif 2677 /* Generate ilim digits, then fix them up. */ 2678 value(eps) *= tens[ilim-1]; 2679 for(i = 1;; i++, value(d) *= 10.) { 2680 L = value(d); 2681 value(d) -= L; 2682 *s++ = '0' + (int)L; 2683 if (i == ilim) { 2684 if (value(d) > 0.5 + value(eps)) 2685 goto bump_up; 2686 else if (value(d) < 0.5 - value(eps)) { 2687 while(*--s == '0'); 2688 s++; 2689 goto ret1; 2690 } 2691 break; 2692 } 2693 } 2694 #ifndef No_leftright 2695 } 2696 #endif 2697 fast_failed: 2698 s = s0; 2699 value(d) = value(d2); 2700 k = k0; 2701 ilim = ilim0; 2702 } 2703 2704 /* Do we have a "small" integer? */ 2705 2706 if (be >= 0 && k <= Int_max) { 2707 /* Yes. */ 2708 ds = tens[k]; 2709 if (ndigits < 0 && ilim <= 0) { 2710 S = mhi = 0; 2711 if (ilim < 0 || value(d) <= 5*ds) 2712 goto no_digits; 2713 goto one_digit; 2714 } 2715 for(i = 1;; i++) { 2716 L = value(d) / ds; 2717 value(d) -= L*ds; 2718 #ifdef Check_FLT_ROUNDS 2719 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 2720 if (value(d) < 0) { 2721 L--; 2722 value(d) += ds; 2723 } 2724 #endif 2725 *s++ = '0' + (int)L; 2726 if (i == ilim) { 2727 value(d) += value(d); 2728 if (value(d) > ds || (value(d) == ds && L & 1)) { 2729 bump_up: 2730 while(*--s == '9') 2731 if (s == s0) { 2732 k++; 2733 *s = '0'; 2734 break; 2735 } 2736 ++*s++; 2737 } 2738 break; 2739 } 2740 if (!(value(d) *= 10.)) 2741 break; 2742 } 2743 goto ret1; 2744 } 2745 2746 m2 = b2; 2747 m5 = b5; 2748 mhi = mlo = 0; 2749 if (leftright) { 2750 if (mode < 2) { 2751 i = 2752 #ifndef Sudden_Underflow 2753 denorm ? be + (Bias + (P-1) - 1 + 1) : 2754 #endif 2755 #ifdef IBM 2756 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 2757 #else 2758 1 + P - bbits; 2759 #endif 2760 } 2761 else { 2762 j = ilim - 1; 2763 if (m5 >= j) 2764 m5 -= j; 2765 else { 2766 s5 += j -= m5; 2767 b5 += j; 2768 m5 = 0; 2769 } 2770 if ((i = ilim) < 0) { 2771 m2 -= i; 2772 i = 0; 2773 } 2774 } 2775 b2 += i; 2776 s2 += i; 2777 mhi = i2b(1); 2778 } 2779 if (m2 > 0 && s2 > 0) { 2780 i = m2 < s2 ? m2 : s2; 2781 b2 -= i; 2782 m2 -= i; 2783 s2 -= i; 2784 } 2785 if (b5 > 0) { 2786 if (leftright) { 2787 if (m5 > 0) { 2788 mhi = pow5mult(mhi, m5); 2789 b1 = mult(mhi, b); 2790 Bfree(b); 2791 b = b1; 2792 } 2793 if ((j = b5 - m5) != 0) 2794 b = pow5mult(b, j); 2795 } 2796 else 2797 b = pow5mult(b, b5); 2798 } 2799 S = i2b(1); 2800 if (s5 > 0) 2801 S = pow5mult(S, s5); 2802 2803 /* Check for special case that d is a normalized power of 2. */ 2804 2805 if (mode < 2) { 2806 if (!word1(d) && !(word0(d) & Bndry_mask) 2807 #ifndef Sudden_Underflow 2808 && word0(d) & Exp_mask 2809 #endif 2810 ) { 2811 /* The special case */ 2812 b2 += Log2P; 2813 s2 += Log2P; 2814 spec_case = 1; 2815 } 2816 else 2817 spec_case = 0; 2818 } 2819 2820 /* Arrange for convenient computation of quotients: 2821 * shift left if necessary so divisor has 4 leading 0 bits. 2822 * 2823 * Perhaps we should just compute leading 28 bits of S once 2824 * and for all and pass them and a shift to quorem, so it 2825 * can do shifts and ors to compute the numerator for q. 2826 */ 2827 if (S == BIGINT_INVALID) { 2828 i = 0; 2829 } else { 2830 #ifdef Pack_32 2831 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0) 2832 i = 32 - i; 2833 #else 2834 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) 2835 i = 16 - i; 2836 #endif 2837 } 2838 2839 if (i > 4) { 2840 i -= 4; 2841 b2 += i; 2842 m2 += i; 2843 s2 += i; 2844 } 2845 else if (i < 4) { 2846 i += 28; 2847 b2 += i; 2848 m2 += i; 2849 s2 += i; 2850 } 2851 if (b2 > 0) 2852 b = lshift(b, b2); 2853 if (s2 > 0) 2854 S = lshift(S, s2); 2855 if (k_check) { 2856 if (cmp(b,S) < 0) { 2857 k--; 2858 b = multadd(b, 10, 0); /* we botched the k estimate */ 2859 if (leftright) 2860 mhi = multadd(mhi, 10, 0); 2861 ilim = ilim1; 2862 } 2863 } 2864 if (ilim <= 0 && mode > 2) { 2865 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 2866 /* no digits, fcvt style */ 2867 no_digits: 2868 k = -1 - ndigits; 2869 goto ret; 2870 } 2871 one_digit: 2872 *s++ = '1'; 2873 k++; 2874 goto ret; 2875 } 2876 if (leftright) { 2877 if (m2 > 0) 2878 mhi = lshift(mhi, m2); 2879 2880 /* Compute mlo -- check for special case 2881 * that d is a normalized power of 2. 2882 */ 2883 2884 mlo = mhi; 2885 if (spec_case) { 2886 mhi = Balloc(mhi->k); 2887 Bcopy(mhi, mlo); 2888 mhi = lshift(mhi, Log2P); 2889 } 2890 2891 for(i = 1;;i++) { 2892 dig = quorem(b,S) + '0'; 2893 /* Do we yet have the shortest decimal string 2894 * that will round to d? 2895 */ 2896 j = cmp(b, mlo); 2897 delta = diff(S, mhi); 2898 jj1 = delta->sign ? 1 : cmp(b, delta); 2899 Bfree(delta); 2900 #ifndef ROUND_BIASED 2901 if (jj1 == 0 && !mode && !(word1(d) & 1)) { 2902 if (dig == '9') 2903 goto round_9_up; 2904 if (j > 0) 2905 dig++; 2906 *s++ = dig; 2907 goto ret; 2908 } 2909 #endif 2910 if (j < 0 || (j == 0 && !mode 2911 #ifndef ROUND_BIASED 2912 && !(word1(d) & 1) 2913 #endif 2914 )) { 2915 if (jj1 > 0) { 2916 b = lshift(b, 1); 2917 jj1 = cmp(b, S); 2918 if ((jj1 > 0 || (jj1 == 0 && dig & 1)) 2919 && dig++ == '9') 2920 goto round_9_up; 2921 } 2922 *s++ = dig; 2923 goto ret; 2924 } 2925 if (jj1 > 0) { 2926 if (dig == '9') { /* possible if i == 1 */ 2927 round_9_up: 2928 *s++ = '9'; 2929 goto roundoff; 2930 } 2931 *s++ = dig + 1; 2932 goto ret; 2933 } 2934 *s++ = dig; 2935 if (i == ilim) 2936 break; 2937 b = multadd(b, 10, 0); 2938 if (mlo == mhi) 2939 mlo = mhi = multadd(mhi, 10, 0); 2940 else { 2941 mlo = multadd(mlo, 10, 0); 2942 mhi = multadd(mhi, 10, 0); 2943 } 2944 } 2945 } 2946 else 2947 for(i = 1;; i++) { 2948 *s++ = dig = quorem(b,S) + '0'; 2949 if (i >= ilim) 2950 break; 2951 b = multadd(b, 10, 0); 2952 } 2953 2954 /* Round off last digit */ 2955 2956 b = lshift(b, 1); 2957 j = cmp(b, S); 2958 if (j > 0 || (j == 0 && dig & 1)) { 2959 roundoff: 2960 while(*--s == '9') 2961 if (s == s0) { 2962 k++; 2963 *s++ = '1'; 2964 goto ret; 2965 } 2966 ++*s++; 2967 } 2968 else { 2969 while(*--s == '0'); 2970 s++; 2971 } 2972 ret: 2973 Bfree(S); 2974 if (mhi) { 2975 if (mlo && mlo != mhi) 2976 Bfree(mlo); 2977 Bfree(mhi); 2978 } 2979 ret1: 2980 Bfree(b); 2981 if (s == s0) { /* don't return empty string */ 2982 *s++ = '0'; 2983 k = 0; 2984 } 2985 *s = 0; 2986 *decpt = k + 1; 2987 if (rve) 2988 *rve = s; 2989 return s0; 2990 } 2991 2992 #include <limits.h> 2993 2994 char * 2995 rv_alloc(int i) 2996 { 2997 int j, k, *r; 2998 2999 j = sizeof(ULong); 3000 for(k = 0; 3001 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i; 3002 j <<= 1) 3003 k++; 3004 r = (int*)Balloc(k); 3005 *r = k; 3006 return (char *)(r+1); 3007 } 3008 3009 char * 3010 nrv_alloc(char *s, char **rve, int n) 3011 { 3012 char *rv, *t; 3013 3014 t = rv = rv_alloc(n); 3015 while((*t = *s++) !=0) 3016 t++; 3017 if (rve) 3018 *rve = t; 3019 return rv; 3020 } 3021 3022 3023 /* Strings values used by dtoa() */ 3024 #define INFSTR "Infinity" 3025 #define NANSTR "NaN" 3026 3027 #define DBL_ADJ (DBL_MAX_EXP - 2 + ((DBL_MANT_DIG - 1) % 4)) 3028 #define LDBL_ADJ (LDBL_MAX_EXP - 2 + ((LDBL_MANT_DIG - 1) % 4)) 3029 3030 /* 3031 * Round up the given digit string. If the digit string is fff...f, 3032 * this procedure sets it to 100...0 and returns 1 to indicate that 3033 * the exponent needs to be bumped. Otherwise, 0 is returned. 3034 */ 3035 static int 3036 roundup(char *s0, int ndigits) 3037 { 3038 char *s; 3039 3040 for (s = s0 + ndigits - 1; *s == 0xf; s--) { 3041 if (s == s0) { 3042 *s = 1; 3043 return (1); 3044 } 3045 *s = 0; 3046 } 3047 ++*s; 3048 return (0); 3049 } 3050 3051 /* 3052 * Round the given digit string to ndigits digits according to the 3053 * current rounding mode. Note that this could produce a string whose 3054 * value is not representable in the corresponding floating-point 3055 * type. The exponent pointed to by decpt is adjusted if necessary. 3056 */ 3057 static void 3058 dorounding(char *s0, int ndigits, int sign, int *decpt) 3059 { 3060 int adjust = 0; /* do we need to adjust the exponent? */ 3061 3062 switch (FLT_ROUNDS) { 3063 case 0: /* toward zero */ 3064 default: /* implementation-defined */ 3065 break; 3066 case 1: /* to nearest, halfway rounds to even */ 3067 if ((s0[ndigits] > 8) || 3068 (s0[ndigits] == 8 && s0[ndigits + 1] & 1)) 3069 adjust = roundup(s0, ndigits); 3070 break; 3071 case 2: /* toward +inf */ 3072 if (sign == 0) 3073 adjust = roundup(s0, ndigits); 3074 break; 3075 case 3: /* toward -inf */ 3076 if (sign != 0) 3077 adjust = roundup(s0, ndigits); 3078 break; 3079 } 3080 3081 if (adjust) 3082 *decpt += 4; 3083 } 3084 3085 /* 3086 * This procedure converts a double-precision number in IEEE format 3087 * into a string of hexadecimal digits and an exponent of 2. Its 3088 * behavior is bug-for-bug compatible with dtoa() in mode 2, with the 3089 * following exceptions: 3090 * 3091 * - An ndigits < 0 causes it to use as many digits as necessary to 3092 * represent the number exactly. 3093 * - The additional xdigs argument should point to either the string 3094 * "0123456789ABCDEF" or the string "0123456789abcdef", depending on 3095 * which case is desired. 3096 * - This routine does not repeat dtoa's mistake of setting decpt 3097 * to 9999 in the case of an infinity or NaN. INT_MAX is used 3098 * for this purpose instead. 3099 * 3100 * Note that the C99 standard does not specify what the leading digit 3101 * should be for non-zero numbers. For instance, 0x1.3p3 is the same 3102 * as 0x2.6p2 is the same as 0x4.cp3. This implementation chooses the 3103 * first digit so that subsequent digits are aligned on nibble 3104 * boundaries (before rounding). 3105 * 3106 * Inputs: d, xdigs, ndigits 3107 * Outputs: decpt, sign, rve 3108 */ 3109 char * 3110 __hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign, 3111 char **rve) 3112 { 3113 static const int sigfigs = (DBL_MANT_DIG + 3) / 4; 3114 union IEEEd2bits u; 3115 char *s, *s0; 3116 int bufsize, f; 3117 3118 u.d = d; 3119 *sign = u.bits.sign; 3120 3121 switch (f = fpclassify(d)) { 3122 case FP_NORMAL: 3123 *decpt = u.bits.exp - DBL_ADJ; 3124 break; 3125 case FP_ZERO: 3126 return_zero: 3127 *decpt = 1; 3128 return (nrv_alloc("0", rve, 1)); 3129 case FP_SUBNORMAL: 3130 /* 3131 * For processors that treat subnormals as zero, comparison 3132 * with zero will be equal, so we jump to the FP_ZERO case. 3133 */ 3134 if(u.d == 0.0) goto return_zero; 3135 u.d *= 0x1p514; 3136 *decpt = u.bits.exp - (514 + DBL_ADJ); 3137 break; 3138 case FP_INFINITE: 3139 *decpt = INT_MAX; 3140 return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1)); 3141 case FP_NAN: 3142 *decpt = INT_MAX; 3143 return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1)); 3144 default: 3145 #ifdef DEBUG 3146 BugPrintf("fpclassify returned %d\n", f); 3147 #endif 3148 return 0; // FIXME?? 3149 } 3150 3151 /* FP_NORMAL or FP_SUBNORMAL */ 3152 3153 if (ndigits == 0) /* dtoa() compatibility */ 3154 ndigits = 1; 3155 3156 /* 3157 * For simplicity, we generate all the digits even if the 3158 * caller has requested fewer. 3159 */ 3160 bufsize = (sigfigs > ndigits) ? sigfigs : ndigits; 3161 s0 = rv_alloc(bufsize); 3162 3163 /* 3164 * We work from right to left, first adding any requested zero 3165 * padding, then the least significant portion of the 3166 * mantissa, followed by the most significant. The buffer is 3167 * filled with the byte values 0x0 through 0xf, which are 3168 * converted to xdigs[0x0] through xdigs[0xf] after the 3169 * rounding phase. 3170 */ 3171 for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--) 3172 *s = 0; 3173 for (; s > s0 + sigfigs - (DBL_MANL_SIZE / 4) - 1 && s > s0; s--) { 3174 *s = u.bits.manl & 0xf; 3175 u.bits.manl >>= 4; 3176 } 3177 for (; s > s0; s--) { 3178 *s = u.bits.manh & 0xf; 3179 u.bits.manh >>= 4; 3180 } 3181 3182 /* 3183 * At this point, we have snarfed all the bits in the 3184 * mantissa, with the possible exception of the highest-order 3185 * (partial) nibble, which is dealt with by the next 3186 * statement. We also tack on the implicit normalization bit. 3187 */ 3188 *s = u.bits.manh | (1U << ((DBL_MANT_DIG - 1) % 4)); 3189 3190 /* If ndigits < 0, we are expected to auto-size the precision. */ 3191 if (ndigits < 0) { 3192 for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--) 3193 ; 3194 } 3195 3196 if (sigfigs > ndigits && s0[ndigits] != 0) 3197 dorounding(s0, ndigits, u.bits.sign, decpt); 3198 3199 s = s0 + ndigits; 3200 if (rve != NULL) 3201 *rve = s; 3202 *s-- = '\0'; 3203 for (; s >= s0; s--) 3204 *s = xdigs[(unsigned int)*s]; 3205 3206 return (s0); 3207 } 3208 3209 #ifndef NO_HEX_FP /*{*/ 3210 3211 static int 3212 gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign, locale_t loc) 3213 { 3214 Bigint *b; 3215 CONST unsigned char *decpt, *s0, *s, *s1; 3216 unsigned char *strunc; 3217 int big, esign, havedig, irv, j, k, n, n0, nbits, up, zret; 3218 ULong L, lostbits, *x; 3219 Long e, e1; 3220 #ifdef USE_LOCALE 3221 int i; 3222 NORMALIZE_LOCALE(loc); 3223 #ifdef NO_LOCALE_CACHE 3224 const unsigned char *decimalpoint = (unsigned char*)localeconv_l(loc)->decimal_point; 3225 #else 3226 const unsigned char *decimalpoint; 3227 static unsigned char *decimalpoint_cache; 3228 if (!(s0 = decimalpoint_cache)) { 3229 s0 = (unsigned char*)localeconv_l(loc)->decimal_point; 3230 if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) { 3231 strcpy(decimalpoint_cache, s0); 3232 s0 = decimalpoint_cache; 3233 } 3234 } 3235 decimalpoint = s0; 3236 #endif 3237 #endif 3238 3239 #ifndef ANDROID_CHANGES 3240 if (!hexdig['0']) 3241 hexdig_init_D2A(); 3242 #endif 3243 3244 *bp = 0; 3245 havedig = 0; 3246 s0 = *(CONST unsigned char **)sp + 2; 3247 while(s0[havedig] == '0') 3248 havedig++; 3249 s0 += havedig; 3250 s = s0; 3251 decpt = 0; 3252 zret = 0; 3253 e = 0; 3254 if (hexdig[*s]) 3255 havedig++; 3256 else { 3257 zret = 1; 3258 #ifdef USE_LOCALE 3259 for(i = 0; decimalpoint[i]; ++i) { 3260 if (s[i] != decimalpoint[i]) 3261 goto pcheck; 3262 } 3263 decpt = s += i; 3264 #else 3265 if (*s != '.') 3266 goto pcheck; 3267 decpt = ++s; 3268 #endif 3269 if (!hexdig[*s]) 3270 goto pcheck; 3271 while(*s == '0') 3272 s++; 3273 if (hexdig[*s]) 3274 zret = 0; 3275 havedig = 1; 3276 s0 = s; 3277 } 3278 while(hexdig[*s]) 3279 s++; 3280 #ifdef USE_LOCALE 3281 if (*s == *decimalpoint && !decpt) { 3282 for(i = 1; decimalpoint[i]; ++i) { 3283 if (s[i] != decimalpoint[i]) 3284 goto pcheck; 3285 } 3286 decpt = s += i; 3287 #else 3288 if (*s == '.' && !decpt) { 3289 decpt = ++s; 3290 #endif 3291 while(hexdig[*s]) 3292 s++; 3293 }/*}*/ 3294 if (decpt) 3295 e = -(((Long)(s-decpt)) << 2); 3296 pcheck: 3297 s1 = s; 3298 big = esign = 0; 3299 switch(*s) { 3300 case 'p': 3301 case 'P': 3302 switch(*++s) { 3303 case '-': 3304 esign = 1; 3305 /* no break */ 3306 case '+': 3307 s++; 3308 } 3309 if ((n = hexdig[*s]) == 0 || n > 0x19) { 3310 s = s1; 3311 break; 3312 } 3313 e1 = n - 0x10; 3314 while((n = hexdig[*++s]) !=0 && n <= 0x19) { 3315 if (e1 & 0xf8000000) 3316 big = 1; 3317 e1 = 10*e1 + n - 0x10; 3318 } 3319 if (esign) 3320 e1 = -e1; 3321 e += e1; 3322 } 3323 *sp = (char*)s; 3324 if (!havedig) 3325 *sp = (char*)s0 - 1; 3326 if (zret) 3327 return STRTOG_Zero; 3328 if (big) { 3329 if (esign) { 3330 switch(fpi->rounding) { 3331 case FPI_Round_up: 3332 if (sign) 3333 break; 3334 goto ret_tiny; 3335 case FPI_Round_down: 3336 if (!sign) 3337 break; 3338 goto ret_tiny; 3339 } 3340 goto retz; 3341 ret_tiny: 3342 b = Balloc(0); 3343 b->wds = 1; 3344 b->x[0] = 1; 3345 goto dret; 3346 } 3347 switch(fpi->rounding) { 3348 case FPI_Round_near: 3349 goto ovfl1; 3350 case FPI_Round_up: 3351 if (!sign) 3352 goto ovfl1; 3353 goto ret_big; 3354 case FPI_Round_down: 3355 if (sign) 3356 goto ovfl1; 3357 goto ret_big; 3358 } 3359 ret_big: 3360 nbits = fpi->nbits; 3361 n0 = n = nbits >> kshift; 3362 if (nbits & kmask) 3363 ++n; 3364 for(j = n, k = 0; j >>= 1; ++k); 3365 *bp = b = Balloc(k); 3366 b->wds = n; 3367 for(j = 0; j < n0; ++j) 3368 b->x[j] = ALL_ON; 3369 if (n > n0) 3370 b->x[j] = ULbits >> (ULbits - (nbits & kmask)); 3371 *exp = fpi->emin; 3372 return STRTOG_Normal | STRTOG_Inexlo; 3373 } 3374 /* 3375 * Truncate the hex string if it is longer than the precision needed, 3376 * to avoid denial-of-service issues with very large strings. Use 3377 * additional digits to insure precision. Scan to-be-truncated digits 3378 * and replace with either '1' or '0' to ensure proper rounding. 3379 */ 3380 { 3381 int maxdigits = ((fpi->nbits + 3) >> 2) + 2; 3382 size_t nd = s1 - s0; 3383 #ifdef USE_LOCALE 3384 int dplen = strlen((const char *)decimalpoint); 3385 #else 3386 int dplen = 1; 3387 #endif 3388 3389 if (decpt && s0 < decpt) 3390 nd -= dplen; 3391 if (nd > maxdigits && (strunc = alloca(maxdigits + dplen + 2)) != NULL) { 3392 ssize_t nd0 = decpt ? decpt - s0 - dplen : nd; 3393 unsigned char *tp = strunc + maxdigits; 3394 int found = 0; 3395 if ((nd0 -= maxdigits) >= 0 || s0 >= decpt) 3396 memcpy(strunc, s0, maxdigits); 3397 else { 3398 memcpy(strunc, s0, maxdigits + dplen); 3399 tp += dplen; 3400 } 3401 s0 += maxdigits; 3402 e += (nd - (maxdigits + 1)) << 2; 3403 if (nd0 > 0) { 3404 while(nd0-- > 0) 3405 if (*s0++ != '0') { 3406 found++; 3407 break; 3408 } 3409 s0 += dplen; 3410 } 3411 if (!found && decpt) { 3412 while(s0 < s1) 3413 if(*s0++ != '0') { 3414 found++; 3415 break; 3416 } 3417 } 3418 *tp++ = found ? '1' : '0'; 3419 *tp = 0; 3420 s0 = strunc; 3421 s1 = tp; 3422 } 3423 } 3424 3425 n = s1 - s0 - 1; 3426 for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1) 3427 k++; 3428 b = Balloc(k); 3429 x = b->x; 3430 n = 0; 3431 L = 0; 3432 #ifdef USE_LOCALE 3433 for(i = 0; decimalpoint[i+1]; ++i); 3434 #endif 3435 while(s1 > s0) { 3436 #ifdef USE_LOCALE 3437 if (*--s1 == decimalpoint[i]) { 3438 s1 -= i; 3439 continue; 3440 } 3441 #else 3442 if (*--s1 == '.') 3443 continue; 3444 #endif 3445 if (n == ULbits) { 3446 *x++ = L; 3447 L = 0; 3448 n = 0; 3449 } 3450 L |= (hexdig[*s1] & 0x0f) << n; 3451 n += 4; 3452 } 3453 *x++ = L; 3454 b->wds = n = x - b->x; 3455 n = ULbits*n - hi0bits(L); 3456 nbits = fpi->nbits; 3457 lostbits = 0; 3458 x = b->x; 3459 if (n > nbits) { 3460 n -= nbits; 3461 if (any_on(b,n)) { 3462 lostbits = 1; 3463 k = n - 1; 3464 if (x[k>>kshift] & 1 << (k & kmask)) { 3465 lostbits = 2; 3466 if (k > 0 && any_on(b,k)) 3467 lostbits = 3; 3468 } 3469 } 3470 rshift(b, n); 3471 e += n; 3472 } 3473 else if (n < nbits) { 3474 n = nbits - n; 3475 b = lshift(b, n); 3476 e -= n; 3477 x = b->x; 3478 } 3479 if (e > fpi->emax) { 3480 ovfl: 3481 Bfree(b); 3482 ovfl1: 3483 #ifndef NO_ERRNO 3484 errno = ERANGE; 3485 #endif 3486 return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; 3487 } 3488 irv = STRTOG_Normal; 3489 if (e < fpi->emin) { 3490 irv = STRTOG_Denormal; 3491 n = fpi->emin - e; 3492 if (n >= nbits) { 3493 switch (fpi->rounding) { 3494 case FPI_Round_near: 3495 if (n == nbits && (n < 2 || any_on(b,n-1))) 3496 goto one_bit; 3497 break; 3498 case FPI_Round_up: 3499 if (!sign) 3500 goto one_bit; 3501 break; 3502 case FPI_Round_down: 3503 if (sign) { 3504 one_bit: 3505 x[0] = b->wds = 1; 3506 dret: 3507 *bp = b; 3508 *exp = fpi->emin; 3509 #ifndef NO_ERRNO 3510 errno = ERANGE; 3511 #endif 3512 return STRTOG_Denormal | STRTOG_Inexhi 3513 | STRTOG_Underflow; 3514 } 3515 } 3516 Bfree(b); 3517 retz: 3518 #ifndef NO_ERRNO 3519 errno = ERANGE; 3520 #endif 3521 return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow; 3522 } 3523 k = n - 1; 3524 if (lostbits) 3525 lostbits = 1; 3526 else if (k > 0) 3527 lostbits = any_on(b,k); 3528 if (x[k>>kshift] & 1 << (k & kmask)) 3529 lostbits |= 2; 3530 nbits -= n; 3531 rshift(b,n); 3532 e = fpi->emin; 3533 } 3534 if (lostbits) { 3535 up = 0; 3536 switch(fpi->rounding) { 3537 case FPI_Round_zero: 3538 break; 3539 case FPI_Round_near: 3540 if (lostbits & 2 3541 && (lostbits | x[0]) & 1) 3542 up = 1; 3543 break; 3544 case FPI_Round_up: 3545 up = 1 - sign; 3546 break; 3547 case FPI_Round_down: 3548 up = sign; 3549 } 3550 if (up) { 3551 k = b->wds; 3552 b = increment(b); 3553 x = b->x; 3554 if (irv == STRTOG_Denormal) { 3555 if (nbits == fpi->nbits - 1 3556 && x[nbits >> kshift] & 1 << (nbits & kmask)) 3557 irv = STRTOG_Normal; 3558 } 3559 else if (b->wds > k 3560 || ((n = nbits & kmask) !=0 3561 && hi0bits(x[k-1]) < 32-n)) { 3562 rshift(b,1); 3563 if (++e > fpi->emax) 3564 goto ovfl; 3565 } 3566 irv |= STRTOG_Inexhi; 3567 } 3568 else 3569 irv |= STRTOG_Inexlo; 3570 } 3571 *bp = b; 3572 *exp = e; 3573 return irv; 3574 } 3575 3576 #endif /*}*/ 3577 3578 #ifdef __cplusplus 3579 } 3580 #endif 3581
