dtoa.c (24850B)
1 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */ 2 /**************************************************************** 3 * 4 * The author of this software is David M. Gay. 5 * 6 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. 7 * 8 * Permission to use, copy, modify, and distribute this software for any 9 * purpose without fee is hereby granted, provided that this entire notice 10 * is included in all copies of any software which is or includes a copy 11 * or modification of this software and in all copies of the supporting 12 * documentation for such software. 13 * 14 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 15 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY 16 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 17 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 18 * 19 ***************************************************************/ 20 21 /* Please send bug reports to David M. Gay (dmg at acm dot org, 22 * with " at " changed at "@" and " dot " changed to "."). */ 23 24 /* On a machine with IEEE extended-precision registers, it is 25 * necessary to specify double-precision (53-bit) rounding precision 26 * before invoking strtod or dtoa. If the machine uses (the equivalent 27 * of) Intel 80x87 arithmetic, the call 28 * _control87(PC_53, MCW_PC); 29 * does this with many compilers. Whether this or another call is 30 * appropriate depends on the compiler; for this to work, it may be 31 * necessary to #include "float.h" or another system-dependent header 32 * file. 33 */ 34 35 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 36 * 37 * This strtod returns a nearest machine number to the input decimal 38 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 39 * broken by the IEEE round-even rule. Otherwise ties are broken by 40 * biased rounding (add half and chop). 41 * 42 * Inspired loosely by William D. Clinger's paper "How to Read Floating 43 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 44 * 45 * Modifications: 46 * 47 * 1. We only require IEEE, IBM, or VAX double-precision 48 * arithmetic (not IEEE double-extended). 49 * 2. We get by with floating-point arithmetic in a case that 50 * Clinger missed -- when we're computing d * 10^n 51 * for a small integer d and the integer n is not too 52 * much larger than 22 (the maximum integer k for which 53 * we can represent 10^k exactly), we may be able to 54 * compute (d*10^k) * 10^(e-k) with just one roundoff. 55 * 3. Rather than a bit-at-a-time adjustment of the binary 56 * result in the hard case, we use floating-point 57 * arithmetic to determine the adjustment to within 58 * one bit; only in really hard cases do we need to 59 * compute a second residual. 60 * 4. Because of 3., we don't need a large table of powers of 10 61 * for ten-to-e (just some small tables, e.g. of 10^k 62 * for 0 <= k <= 22). 63 */ 64 65 /* 66 * #define IEEE_8087 for IEEE-arithmetic machines where the least 67 * significant byte has the lowest address. 68 * #define IEEE_MC68k for IEEE-arithmetic machines where the most 69 * significant byte has the lowest address. 70 * #define Long int on machines with 32-bit ints and 64-bit longs. 71 * #define IBM for IBM mainframe-style floating-point arithmetic. 72 * #define VAX for VAX-style floating-point arithmetic (D_floating). 73 * #define No_leftright to omit left-right logic in fast floating-point 74 * computation of dtoa. 75 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 76 * and strtod and dtoa should round accordingly. 77 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 78 * and Honor_FLT_ROUNDS is not #defined. 79 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 80 * that use extended-precision instructions to compute rounded 81 * products and quotients) with IBM. 82 * #define ROUND_BIASED for IEEE-format with biased rounding. 83 * #define Inaccurate_Divide for IEEE-format with correctly rounded 84 * products but inaccurate quotients, e.g., for Intel i860. 85 * #define NO_LONG_LONG on machines that do not have a "long long" 86 * integer type (of >= 64 bits). On such machines, you can 87 * #define Just_16 to store 16 bits per 32-bit Long when doing 88 * high-precision integer arithmetic. Whether this speeds things 89 * up or slows things down depends on the machine and the number 90 * being converted. If long long is available and the name is 91 * something other than "long long", #define Llong to be the name, 92 * and if "unsigned Llong" does not work as an unsigned version of 93 * Llong, #define #ULLong to be the corresponding unsigned type. 94 * #define KR_headers for old-style C function headers. 95 * #define Bad_float_h if your system lacks a float.h or if it does not 96 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, 97 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. 98 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) 99 * if memory is available and otherwise does something you deem 100 * appropriate. If MALLOC is undefined, malloc will be invoked 101 * directly -- and assumed always to succeed. Similarly, if you 102 * want something other than the system's free() to be called to 103 * recycle memory acquired from MALLOC, #define FREE to be the 104 * name of the alternate routine. (Unless you #define 105 * NO_GLOBAL_STATE and call destroydtoa, FREE or free is only 106 * called in pathological cases, e.g., in a dtoa call after a dtoa 107 * return in mode 3 with thousands of digits requested.) 108 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making 109 * memory allocations from a private pool of memory when possible. 110 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, 111 * unless #defined to be a different length. This default length 112 * suffices to get rid of MALLOC calls except for unusual cases, 113 * such as decimal-to-binary conversion of a very long string of 114 * digits. The longest string dtoa can return is about 751 bytes 115 * long. For conversions by strtod of strings of 800 digits and 116 * all dtoa conversions in single-threaded executions with 8-byte 117 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte 118 * pointers, PRIVATE_MEM >= 7112 appears adequate. 119 * #define MULTIPLE_THREADS if the system offers preemptively scheduled 120 * multiple threads. In this case, you must provide (or suitably 121 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed 122 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed 123 * in pow5mult, ensures lazy evaluation of only one copy of high 124 * powers of 5; omitting this lock would introduce a small 125 * probability of wasting memory, but would otherwise be harmless.) 126 * You must also invoke freedtoa(s) to free the value s returned by 127 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. 128 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that 129 * avoids underflows on inputs whose result does not underflow. 130 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format 131 * floating-point numbers and flushes underflows to zero rather 132 * than implementing gradual underflow, then you must also #define 133 * Sudden_Underflow. 134 * #define USE_LOCALE to use the current locale's decimal_point value. 135 * #define SET_INEXACT if IEEE arithmetic is being used and extra 136 * computation should be done to set the inexact flag when the 137 * result is inexact and avoid setting inexact when the result 138 * is exact. In this case, dtoa.c must be compiled in 139 * an environment, perhaps provided by #include "dtoa.c" in a 140 * suitable wrapper, that defines two functions, 141 * int get_inexact(void); 142 * void clear_inexact(void); 143 * such that get_inexact() returns a nonzero value if the 144 * inexact bit is already set, and clear_inexact() sets the 145 * inexact bit to 0. When SET_INEXACT is #defined, strtod 146 * also does extra computations to set the underflow and overflow 147 * flags when appropriate (i.e., when the result is tiny and 148 * inexact or when it is a numeric value rounded to +-infinity). 149 * #define NO_ERRNO if strtod should not assign errno = ERANGE when 150 * the result overflows to +-Infinity or underflows to 0. 151 * #define NO_GLOBAL_STATE to avoid defining any non-const global or 152 * static variables. Instead the necessary state is stored in an 153 * opaque struct, DtoaState, a pointer to which must be passed to 154 * every entry point. Two new functions are added to the API: 155 * DtoaState *newdtoa(void); 156 * void destroydtoa(DtoaState *); 157 */ 158 159 #ifndef Long 160 #define Long long 161 #endif 162 #ifndef ULong 163 typedef unsigned Long ULong; 164 #endif 165 166 #ifdef DEBUG 167 #include <stdio.h> 168 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 169 #endif 170 171 #include <stdlib.h> 172 #include <string.h> 173 174 #ifdef USE_LOCALE 175 #include <locale.h> 176 #endif 177 178 #ifdef MALLOC 179 #ifdef KR_headers 180 extern char *MALLOC(); 181 #else 182 extern void *MALLOC(size_t); 183 #endif 184 #else 185 #define MALLOC malloc 186 #endif 187 188 #ifndef FREE 189 #define FREE free 190 #endif 191 192 #ifndef Omit_Private_Memory 193 #ifndef PRIVATE_MEM 194 #define PRIVATE_MEM 2304 195 #endif 196 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) 197 #endif 198 199 #undef IEEE_Arith 200 #undef Avoid_Underflow 201 #ifdef IEEE_MC68k 202 #define IEEE_Arith 203 #endif 204 #ifdef IEEE_8087 205 #define IEEE_Arith 206 #endif 207 208 #include <errno.h> 209 210 #ifdef Bad_float_h 211 212 #ifdef IEEE_Arith 213 #define DBL_DIG 15 214 #define DBL_MAX_10_EXP 308 215 #define DBL_MAX_EXP 1024 216 #define FLT_RADIX 2 217 #endif /*IEEE_Arith*/ 218 219 #ifdef IBM 220 #define DBL_DIG 16 221 #define DBL_MAX_10_EXP 75 222 #define DBL_MAX_EXP 63 223 #define FLT_RADIX 16 224 #define DBL_MAX 7.2370055773322621e+75 225 #endif 226 227 #ifdef VAX 228 #define DBL_DIG 16 229 #define DBL_MAX_10_EXP 38 230 #define DBL_MAX_EXP 127 231 #define FLT_RADIX 2 232 #define DBL_MAX 1.7014118346046923e+38 233 #endif 234 235 #ifndef LONG_MAX 236 #define LONG_MAX 2147483647 237 #endif 238 239 #else /* ifndef Bad_float_h */ 240 #include <float.h> 241 #endif /* Bad_float_h */ 242 243 #ifndef __MATH_H__ 244 #include <math.h> 245 #endif 246 247 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 248 #error "Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined." 249 #endif 250 251 typedef union { double d; ULong L[2]; } U; 252 253 #define dval(x) ((x).d) 254 #ifdef IEEE_8087 255 #define word0(x) ((x).L[1]) 256 #define word1(x) ((x).L[0]) 257 #else 258 #define word0(x) ((x).L[0]) 259 #define word1(x) ((x).L[1]) 260 #endif 261 262 /* The following definition of Storeinc is appropriate for MIPS processors. 263 * An alternative that might be better on some machines is 264 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 265 */ 266 #if defined(IEEE_8087) + defined(VAX) 267 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 268 ((unsigned short *)a)[0] = (unsigned short)c, a++) 269 #else 270 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 271 ((unsigned short *)a)[1] = (unsigned short)c, a++) 272 #endif 273 274 /* #define P DBL_MANT_DIG */ 275 /* Ten_pmax = floor(P*log(2)/log(5)) */ 276 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 277 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 278 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 279 280 #ifdef IEEE_Arith 281 #define Exp_shift 20 282 #define Exp_shift1 20 283 #define Exp_msk1 0x100000 284 #define Exp_msk11 0x100000 285 #define Exp_mask 0x7ff00000 286 #define P 53 287 #define Bias 1023 288 #define Emin (-1022) 289 #define Exp_1 0x3ff00000 290 #define Exp_11 0x3ff00000 291 #define Ebits 11 292 #define Frac_mask 0xfffff 293 #define Frac_mask1 0xfffff 294 #define Ten_pmax 22 295 #define Bletch 0x10 296 #define Bndry_mask 0xfffff 297 #define Bndry_mask1 0xfffff 298 #define LSB 1 299 #define Sign_bit 0x80000000 300 #define Log2P 1 301 #define Tiny0 0 302 #define Tiny1 1 303 #define Quick_max 14 304 #define Int_max 14 305 #ifndef NO_IEEE_Scale 306 #define Avoid_Underflow 307 #ifdef Flush_Denorm /* debugging option */ 308 #undef Sudden_Underflow 309 #endif 310 #endif 311 312 #ifndef Flt_Rounds 313 #ifdef FLT_ROUNDS 314 #define Flt_Rounds FLT_ROUNDS 315 #else 316 #define Flt_Rounds 1 317 #endif 318 #endif /*Flt_Rounds*/ 319 320 #ifdef Honor_FLT_ROUNDS 321 #define Rounding rounding 322 #undef Check_FLT_ROUNDS 323 #define Check_FLT_ROUNDS 324 #else 325 #define Rounding Flt_Rounds 326 #endif 327 328 #else /* ifndef IEEE_Arith */ 329 #undef Check_FLT_ROUNDS 330 #undef Honor_FLT_ROUNDS 331 #undef SET_INEXACT 332 #undef Sudden_Underflow 333 #define Sudden_Underflow 334 #ifdef IBM 335 #undef Flt_Rounds 336 #define Flt_Rounds 0 337 #define Exp_shift 24 338 #define Exp_shift1 24 339 #define Exp_msk1 0x1000000 340 #define Exp_msk11 0x1000000 341 #define Exp_mask 0x7f000000 342 #define P 14 343 #define Bias 65 344 #define Exp_1 0x41000000 345 #define Exp_11 0x41000000 346 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 347 #define Frac_mask 0xffffff 348 #define Frac_mask1 0xffffff 349 #define Bletch 4 350 #define Ten_pmax 22 351 #define Bndry_mask 0xefffff 352 #define Bndry_mask1 0xffffff 353 #define LSB 1 354 #define Sign_bit 0x80000000 355 #define Log2P 4 356 #define Tiny0 0x100000 357 #define Tiny1 0 358 #define Quick_max 14 359 #define Int_max 15 360 #else /* VAX */ 361 #undef Flt_Rounds 362 #define Flt_Rounds 1 363 #define Exp_shift 23 364 #define Exp_shift1 7 365 #define Exp_msk1 0x80 366 #define Exp_msk11 0x800000 367 #define Exp_mask 0x7f80 368 #define P 56 369 #define Bias 129 370 #define Exp_1 0x40800000 371 #define Exp_11 0x4080 372 #define Ebits 8 373 #define Frac_mask 0x7fffff 374 #define Frac_mask1 0xffff007f 375 #define Ten_pmax 24 376 #define Bletch 2 377 #define Bndry_mask 0xffff007f 378 #define Bndry_mask1 0xffff007f 379 #define LSB 0x10000 380 #define Sign_bit 0x8000 381 #define Log2P 1 382 #define Tiny0 0x80 383 #define Tiny1 0 384 #define Quick_max 15 385 #define Int_max 15 386 #endif /* IBM, VAX */ 387 #endif /* IEEE_Arith */ 388 389 #ifndef IEEE_Arith 390 #define ROUND_BIASED 391 #endif 392 393 #ifdef RND_PRODQUOT 394 #define rounded_product(a,b) a = rnd_prod(a, b) 395 #define rounded_quotient(a,b) a = rnd_quot(a, b) 396 #ifdef KR_headers 397 extern double rnd_prod(), rnd_quot(); 398 #else 399 extern double rnd_prod(double, double), rnd_quot(double, double); 400 #endif 401 #else 402 #define rounded_product(a,b) a *= b 403 #define rounded_quotient(a,b) a /= b 404 #endif 405 406 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 407 #define Big1 0xffffffff 408 409 #ifndef Pack_32 410 #define Pack_32 411 #endif 412 413 #ifdef KR_headers 414 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff) 415 #else 416 #define FFFFFFFF 0xffffffffUL 417 #endif 418 419 #ifdef NO_LONG_LONG 420 #undef ULLong 421 #ifdef Just_16 422 #undef Pack_32 423 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long. 424 * This makes some inner loops simpler and sometimes saves work 425 * during multiplications, but it often seems to make things slightly 426 * slower. Hence the default is now to store 32 bits per Long. 427 */ 428 #endif 429 #else /* long long available */ 430 #ifndef Llong 431 #define Llong long long 432 #endif 433 #ifndef ULLong 434 #define ULLong unsigned Llong 435 #endif 436 #endif /* NO_LONG_LONG */ 437 438 #ifndef MULTIPLE_THREADS 439 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/ 440 #define FREE_DTOA_LOCK(n) /*nothing*/ 441 #endif 442 443 #define Kmax 7 444 445 struct 446 Bigint { 447 struct Bigint *next; 448 int k, maxwds, sign, wds; 449 ULong x[1]; 450 }; 451 452 typedef struct Bigint Bigint; 453 454 #ifdef NO_GLOBAL_STATE 455 #ifdef MULTIPLE_THREADS 456 #error "cannot have both NO_GLOBAL_STATE and MULTIPLE_THREADS" 457 #endif 458 struct 459 DtoaState { 460 #define DECLARE_GLOBAL_STATE /* nothing */ 461 #else 462 #define DECLARE_GLOBAL_STATE static 463 #endif 464 465 DECLARE_GLOBAL_STATE Bigint *freelist[Kmax+1]; 466 DECLARE_GLOBAL_STATE Bigint *p5s; 467 #ifndef Omit_Private_Memory 468 DECLARE_GLOBAL_STATE double private_mem[PRIVATE_mem]; 469 DECLARE_GLOBAL_STATE double *pmem_next 470 #ifndef NO_GLOBAL_STATE 471 = private_mem 472 #endif 473 ; 474 #endif 475 #ifdef NO_GLOBAL_STATE 476 }; 477 typedef struct DtoaState DtoaState; 478 #ifdef KR_headers 479 #define STATE_PARAM state, 480 #define STATE_PARAM_DECL DtoaState *state; 481 #else 482 #define STATE_PARAM DtoaState *state, 483 #endif 484 #define PASS_STATE state, 485 #define GET_STATE(field) (state->field) 486 487 static DtoaState * 488 newdtoa(void) 489 { 490 DtoaState *state = (DtoaState *) MALLOC(sizeof(DtoaState)); 491 if (state) { 492 memset(state, 0, sizeof(DtoaState)); 493 #ifndef Omit_Private_Memory 494 state->pmem_next = state->private_mem; 495 #endif 496 } 497 return state; 498 } 499 500 static void 501 destroydtoa 502 #ifdef KR_headers 503 (state) STATE_PARAM_DECL 504 #else 505 (DtoaState *state) 506 #endif 507 { 508 int i; 509 Bigint *v, *next; 510 511 for (i = 0; i <= Kmax; i++) { 512 for (v = GET_STATE(freelist)[i]; v; v = next) { 513 next = v->next; 514 #ifndef Omit_Private_Memory 515 if ((double*)v < GET_STATE(private_mem) || 516 (double*)v >= GET_STATE(private_mem) + PRIVATE_mem) 517 #endif 518 FREE((void*)v); 519 } 520 } 521 #ifdef Omit_Private_Memory 522 Bigint* p5 = GET_STATE(p5s); 523 while (p5) { 524 Bigint* tmp = p5; 525 p5 = p5->next; 526 FREE(tmp); 527 } 528 #endif 529 FREE((void *)state); 530 } 531 532 #else 533 #define STATE_PARAM /* nothing */ 534 #define STATE_PARAM_DECL /* nothing */ 535 #define PASS_STATE /* nothing */ 536 #define GET_STATE(name) name 537 #endif 538 539 static Bigint * 540 Balloc 541 #ifdef KR_headers 542 (STATE_PARAM k) STATE_PARAM_DECL int k; 543 #else 544 (STATE_PARAM int k) 545 #endif 546 { 547 int x; 548 Bigint *rv; 549 #ifndef Omit_Private_Memory 550 size_t len; 551 #endif 552 553 ACQUIRE_DTOA_LOCK(0); 554 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */ 555 /* but this case seems very unlikely. */ 556 if (k <= Kmax && (rv = GET_STATE(freelist)[k])) 557 GET_STATE(freelist)[k] = rv->next; 558 else { 559 x = 1 << k; 560 #ifdef Omit_Private_Memory 561 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); 562 #else 563 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 564 /sizeof(double); 565 if (k <= Kmax && GET_STATE(pmem_next) - GET_STATE(private_mem) + len <= PRIVATE_mem) { 566 rv = (Bigint*)GET_STATE(pmem_next); 567 GET_STATE(pmem_next) += len; 568 } 569 else 570 rv = (Bigint*)MALLOC(len*sizeof(double)); 571 #endif 572 rv->k = k; 573 rv->maxwds = x; 574 } 575 FREE_DTOA_LOCK(0); 576 rv->sign = rv->wds = 0; 577 return rv; 578 } 579 580 static void 581 Bfree 582 #ifdef KR_headers 583 (STATE_PARAM v) STATE_PARAM_DECL Bigint *v; 584 #else 585 (STATE_PARAM Bigint *v) 586 #endif 587 { 588 if (v) { 589 if (v->k > Kmax) 590 FREE((void*)v); 591 else { 592 ACQUIRE_DTOA_LOCK(0); 593 v->next = GET_STATE(freelist)[v->k]; 594 GET_STATE(freelist)[v->k] = v; 595 FREE_DTOA_LOCK(0); 596 } 597 } 598 } 599 600 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 601 y->wds*sizeof(Long) + 2*sizeof(int)) 602 603 static Bigint * 604 multadd 605 #ifdef KR_headers 606 (STATE_PARAM b, m, a) STATE_PARAM_DECL Bigint *b; int m, a; 607 #else 608 (STATE_PARAM Bigint *b, int m, int a) /* multiply by m and add a */ 609 #endif 610 { 611 int i, wds; 612 #ifdef ULLong 613 ULong *x; 614 ULLong carry, y; 615 #else 616 ULong carry, *x, y; 617 #ifdef Pack_32 618 ULong xi, z; 619 #endif 620 #endif 621 Bigint *b1; 622 623 wds = b->wds; 624 x = b->x; 625 i = 0; 626 carry = a; 627 do { 628 #ifdef ULLong 629 y = *x * (ULLong)m + carry; 630 carry = y >> 32; 631 *x++ = (ULong) y & FFFFFFFF; 632 #else 633 #ifdef Pack_32 634 xi = *x; 635 y = (xi & 0xffff) * m + carry; 636 z = (xi >> 16) * m + (y >> 16); 637 carry = z >> 16; 638 *x++ = (z << 16) + (y & 0xffff); 639 #else 640 y = *x * m + carry; 641 carry = y >> 16; 642 *x++ = y & 0xffff; 643 #endif 644 #endif 645 } 646 while(++i < wds); 647 if (carry) { 648 if (wds >= b->maxwds) { 649 b1 = Balloc(PASS_STATE b->k+1); 650 Bcopy(b1, b); 651 Bfree(PASS_STATE b); 652 b = b1; 653 } 654 b->x[wds++] = (ULong) carry; 655 b->wds = wds; 656 } 657 return b; 658 } 659 660 static int 661 hi0bits 662 #ifdef KR_headers 663 (x) ULong x; 664 #else 665 (ULong x) 666 #endif 667 { 668 int k = 0; 669 670 if (!(x & 0xffff0000)) { 671 k = 16; 672 x <<= 16; 673 } 674 if (!(x & 0xff000000)) { 675 k += 8; 676 x <<= 8; 677 } 678 if (!(x & 0xf0000000)) { 679 k += 4; 680 x <<= 4; 681 } 682 if (!(x & 0xc0000000)) { 683 k += 2; 684 x <<= 2; 685 } 686 if (!(x & 0x80000000)) { 687 k++; 688 if (!(x & 0x40000000)) 689 return 32; 690 } 691 return k; 692 } 693 694 static int 695 lo0bits 696 #ifdef KR_headers 697 (y) ULong *y; 698 #else 699 (ULong *y) 700 #endif 701 { 702 int k; 703 ULong x = *y; 704 705 if (x & 7) { 706 if (x & 1) 707 return 0; 708 if (x & 2) { 709 *y = x >> 1; 710 return 1; 711 } 712 *y = x >> 2; 713 return 2; 714 } 715 k = 0; 716 if (!(x & 0xffff)) { 717 k = 16; 718 x >>= 16; 719 } 720 if (!(x & 0xff)) { 721 k += 8; 722 x >>= 8; 723 } 724 if (!(x & 0xf)) { 725 k += 4; 726 x >>= 4; 727 } 728 if (!(x & 0x3)) { 729 k += 2; 730 x >>= 2; 731 } 732 if (!(x & 1)) { 733 k++; 734 x >>= 1; 735 if (!x) 736 return 32; 737 } 738 *y = x; 739 return k; 740 } 741 742 static Bigint * 743 i2b 744 #ifdef KR_headers 745 (STATE_PARAM i) STATE_PARAM_DECL int i; 746 #else 747 (STATE_PARAM int i) 748 #endif 749 { 750 Bigint *b; 751 752 b = Balloc(PASS_STATE 1); 753 b->x[0] = i; 754 b->wds = 1; 755 return b; 756 } 757 758 static Bigint * 759 lshift 760 #ifdef KR_headers 761 (STATE_PARAM b, k) STATE_PARAM_DECL Bigint *b; int k; 762 #else 763 (STATE_PARAM Bigint *b, int k) 764 #endif 765 { 766 int i, k1, n, n1; 767 Bigint *b1; 768 ULong *x, *x1, *xe, z; 769 770 #ifdef Pack_32 771 n = k >> 5; 772 #else 773 n = k >> 4; 774 #endif 775 k1 = b->k; 776 n1 = n + b->wds + 1; 777 for(i = b->maxwds; n1 > i; i <<= 1) 778 k1++; 779 b1 = Balloc(PASS_STATE k1); 780 x1 = b1->x; 781 for(i = 0; i < n; i++) 782 *x1++ = 0; 783 x = b->x; 784 xe = x + b->wds; 785 #ifdef Pack_32 786 if (k &= 0x1f) { 787 k1 = 32 - k; 788 z = 0; 789 do { 790 *x1++ = *x << k | z; 791 z = *x++ >> k1; 792 } 793 while(x < xe); 794 if ((*x1 = z)) 795 ++n1; 796 } 797 #else 798 if (k &= 0xf) { 799 k1 = 16 - k; 800 z = 0; 801 do { 802 *x1++ = *x << k & 0xffff | z; 803 z = *x++ >> k1; 804 } 805 while(x < xe); 806 if (*x1 = z) 807 ++n1; 808 } 809 #endif 810 else do 811 *x1++ = *x++; 812 while(x < xe); 813 b1->wds = n1 - 1; 814 Bfree(PASS_STATE b); 815 return b1; 816 } 817 818 static int 819 cmp 820 #ifdef KR_headers 821 (a, b) Bigint *a, *b; 822 #else 823 (Bigint *a, Bigint *b) 824 #endif 825 { 826 ULong *xa, *xa0, *xb, *xb0; 827 int i, j; 828 829 i = a->wds; 830 j = b->wds; 831 #ifdef DEBUG 832 if (i > 1 && !a->x[i-1]) 833 Bug("cmp called with a->x[a->wds-1] == 0"); 834 if (j > 1 && !b->x[j-1]) 835 Bug("cmp called with b->x[b->wds-1] == 0"); 836 #endif 837 if (i -= j) 838 return i; 839 xa0 = a->x; 840 xa = xa0 + j; 841 xb0 = b->x; 842 xb = xb0 + j; 843 for(;;) { 844 if (*--xa != *--xb) 845 return *xa < *xb ? -1 : 1; 846 if (xa <= xa0) 847 break; 848 } 849 return 0; 850 } 851 852 static Bigint * 853 diff 854 #ifdef KR_headers 855 (STATE_PARAM a, b) STATE_PARAM_DECL Bigint *a, *b; 856 #else 857 (STATE_PARAM Bigint *a, Bigint *b) 858 #endif 859 { 860 Bigint *c; 861 int i, wa, wb; 862 ULong *xa, *xae, *xb, *xbe, *xc; 863 #ifdef ULLong 864 ULLong borrow, y; 865 #else 866 ULong borrow, y; 867 #ifdef Pack_32 868 ULong z; 869 #endif 870 #endif 871 872 i = cmp(a,b); 873 if (!i) { 874 c = Balloc(PASS_STATE 0); 875 c->wds = 1; 876 c->x[0] = 0; 877 return c; 878 } 879 if (i < 0) { 880 c = a; 881 a = b; 882 b = c; 883 i = 1; 884 } 885 else 886 i = 0; 887 c = Balloc(PASS_STATE a->k); 888 c->sign = i; 889 wa = a->wds; 890 xa = a->x; 891 xae = xa + wa; 892 wb = b->wds; 893 xb = b->x; 894 xbe = xb + wb; 895 xc = c->x; 896 borrow = 0; 897 #ifdef ULLong 898 do { 899 y = (ULLong)*xa++ - *xb++ - borrow; 900 borrow = y >> 32 & (ULong)1; 901 *xc++ = (ULong) y & FFFFFFFF; 902 } 903 while(xb < xbe); 904 while(xa < xae) { 905 y = *xa++ - borrow; 906 borrow = y >> 32 & (ULong)1; 907 *xc++ = (ULong) y & FFFFFFFF; 908 } 909 #else 910 #ifdef Pack_32 911 do { 912 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 913 borrow = (y & 0x10000) >> 16; 914 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 915 borrow = (z & 0x10000) >> 16; 916 Storeinc(xc, z, y); 917 } 918 while(xb < xbe); 919 while(xa < xae) { 920 y = (*xa & 0xffff) - borrow; 921 borrow = (y & 0x10000) >> 16; 922 z = (*xa++ >> 16) - borrow; 923 borrow = (z & 0x10000) >> 16; 924 Storeinc(xc, z, y); 925 } 926 #else 927 do { 928 y = *xa++ - *xb++ - borrow; 929 borrow = (y & 0x10000) >> 16; 930 *xc++ = y & 0xffff; 931 } 932 while(xb < xbe); 933 while(xa < xae) { 934 y = *xa++ - borrow; 935 borrow = (y & 0x10000) >> 16; 936 *xc++ = y & 0xffff; 937 } 938 #endif 939 #endif 940 while(!*--xc) 941 wa--; 942 c->wds = wa; 943 return c; 944 } 945 946 static Bigint * 947 d2b 948 #ifdef KR_headers 949 (STATE_PARAM d, e, bits) STATE_PARAM_DECL U d; int *e, *bits; 950 #else 951 (STATE_PARAM U d, int *e, int *bits) 952 #endif 953 { 954 Bigint *b; 955 int de, k; 956 ULong *x, y, z; 957 #ifndef Sudden_Underflow 958 int i; 959 #endif 960 #ifdef VAX 961 ULong d0, d1; 962 d0 = word0(d) >> 16 | word0(d) << 16; 963 d1 = word1(d) >> 16 | word1(d) << 16; 964 #else 965 #define d0 word0(d) 966 #define d1 word1(d) 967 #endif 968 969 #ifdef Pack_32 970 b = Balloc(PASS_STATE 1); 971 #else 972 b = Balloc(PASS_STATE 2); 973 #endif 974 x = b->x; 975 976 z = d0 & Frac_mask; 977 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 978 #ifdef Sudden_Underflow 979 de = (int)(d0 >> Exp_shift); 980 #ifndef IBM 981 z |= Exp_msk11; 982 #endif 983 #else 984 if ((de = (int)(d0 >> Exp_shift))) 985 z |= Exp_msk1; 986 #endif 987 #ifdef Pack_32 988 if ((y = d1)) { 989 if ((k = lo0bits(&y))) { 990 x[0] = y | z << (32 - k); 991 z >>= k; 992 } 993 else 994 x[0] = y; 995 #ifndef Sudden_Underflow 996 i = 997 #endif 998 b->wds = (x[1] = z) ? 2 : 1; 999 } 1000 else { 1001 k = lo0bits(&z); 1002 x[0] = z; 1003 #ifndef Sudden_Underflow 1004 i = 1005 #endif 1006 b->wds = 1; 1007 k += 32; 1008 } 1009 #else 1010 if (y = d1) { 1011 if (k = lo0bits(&y)) 1012 if (k >= 16) { 1013 x[0] = y | z << 32 - k & 0xffff; 1014 x[1] = z >> k - 16 & 0xffff; 1015 x[2] = z >> k; 1016 i = 2; 1017 } 1018 else { 1019 x[0] = y & 0xffff; 1020 x[1] = y >> 16 | z << 16 - k & 0xffff; 1021 x[2] = z >> k & 0xffff; 1022 x[3] = z >> k+16; 1023 i = 3; 1024 } 1025 else { 1026 x[0] = y & 0xffff; 1027 x[1] = y >> 16; 1028 x[2] = z & 0xffff; 1029 x[3] = z >> 16; 1030 i = 3; 1031 } 1032 } 1033 else { 1034 #ifdef DEBUG 1035 if (!z) 1036 Bug("Zero passed to d2b"); 1037 #endif 1038 k = lo0bits(&z); 1039 if (k >= 16) { 1040 x[0] = z; 1041 i = 0; 1042 } 1043 else { 1044 x[0] = z & 0xffff; 1045 x[1] = z >> 16; 1046 i = 1; 1047 } 1048 k += 32; 1049 } 1050 while(!x[i]) 1051 --i; 1052 b->wds = i + 1; 1053 #endif 1054 #ifndef Sudden_Underflow 1055 if (de) { 1056 #endif 1057 #ifdef IBM 1058 *e = (de - Bias - (P-1) << 2) + k; 1059 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 1060 #else 1061 *e = de - Bias - (P-1) + k; 1062 *bits = P - k; 1063 #endif 1064 #ifndef Sudden_Underflow 1065 } 1066 else { 1067 *e = de - Bias - (P-1) + 1 + k; 1068 #ifdef Pack_32 1069 *bits = 32*i - hi0bits(x[i-1]); 1070 #else 1071 *bits = (i+2)*16 - hi0bits(x[i]); 1072 #endif 1073 } 1074 #endif 1075 return b; 1076 } 1077 #undef d0 1078 #undef d1