00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172 #include "dtoa.h"
00173 #include <config.h>
00174
00175 #include "global.h"
00176
00177
00178
00179
00180 #define IEEE_8087
00181
00182 #define INFNAN_CHECK
00183
00184
00185
00186 #ifndef Long
00187 #define Long int
00188 #endif
00189 #ifndef ULong
00190 typedef unsigned Long ULong;
00191 #endif
00192
00193 #ifdef DEBUG
00194 #include <stdio.h>
00195 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
00196 #endif
00197
00198 #include <stdlib.h>
00199 #include <string.h>
00200
00201 #ifdef USE_LOCALE
00202 #include <locale.h>
00203 #endif
00204
00205 #ifdef MALLOC
00206 extern void *MALLOC(size_t);
00207 #else
00208 #define MALLOC malloc
00209 #endif
00210
00211 #ifndef Omit_Private_Memory
00212 #ifndef PRIVATE_MEM
00213 #define PRIVATE_MEM 2304
00214 #endif
00215 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
00216 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
00217 #endif
00218
00219 #undef IEEE_Arith
00220 #undef Avoid_Underflow
00221 #ifdef IEEE_MC68k
00222 #define IEEE_Arith
00223 #endif
00224 #ifdef IEEE_8087
00225 #define IEEE_Arith
00226 #endif
00227
00228 #include <errno.h>
00229
00230 #ifdef Bad_float_h
00231
00232 #ifdef IEEE_Arith
00233 #define DBL_DIG 15
00234 #define DBL_MAX_10_EXP 308
00235 #define DBL_MAX_EXP 1024
00236 #define FLT_RADIX 2
00237 #endif
00238
00239 #ifdef IBM
00240 #define DBL_DIG 16
00241 #define DBL_MAX_10_EXP 75
00242 #define DBL_MAX_EXP 63
00243 #define FLT_RADIX 16
00244 #define DBL_MAX 7.2370055773322621e+75
00245 #endif
00246
00247 #ifdef VAX
00248 #define DBL_DIG 16
00249 #define DBL_MAX_10_EXP 38
00250 #define DBL_MAX_EXP 127
00251 #define FLT_RADIX 2
00252 #define DBL_MAX 1.7014118346046923e+38
00253 #endif
00254
00255 #ifndef LONG_MAX
00256 #define LONG_MAX 2147483647
00257 #endif
00258
00259 #else
00260 #include <float.h>
00261 #endif
00262
00263 #ifndef __MATH_H__
00264 #include <math.h>
00265 #endif
00266
00267 #define strtod kjs_strtod
00268 #define dtoa kjs_dtoa
00269 #define freedtoa kjs_freedtoa
00270
00271 #ifdef __cplusplus
00272 extern "C" {
00273 #endif
00274
00275
00276 #define CONST const
00277
00278
00279 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
00280 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
00281 #endif
00282
00283 typedef union { double d; ULong L[2]; } U;
00284
00285 #define dval(x) (x).d
00286 #ifdef IEEE_8087
00287 #define word0(x) (x).L[1]
00288 #define word1(x) (x).L[0]
00289 #else
00290 #define word0(x) (x).L[0]
00291 #define word1(x) (x).L[1]
00292 #endif
00293
00294
00295
00296
00297 #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
00298
00299
00300
00301
00302
00303
00304
00305 #ifdef IEEE_Arith
00306 #define Exp_shift 20
00307 #define Exp_shift1 20
00308 #define Exp_msk1 0x100000
00309 #define Exp_msk11 0x100000
00310 #define Exp_mask 0x7ff00000
00311 #define P 53
00312 #define Bias 1023
00313 #define Emin (-1022)
00314 #define Exp_1 0x3ff00000
00315 #define Exp_11 0x3ff00000
00316 #define Ebits 11
00317 #define Frac_mask 0xfffff
00318 #define Frac_mask1 0xfffff
00319 #define Ten_pmax 22
00320 #define Bletch 0x10
00321 #define Bndry_mask 0xfffff
00322 #define Bndry_mask1 0xfffff
00323 #define LSB 1
00324 #define Sign_bit 0x80000000
00325 #define Log2P 1
00326 #define Tiny0 0
00327 #define Tiny1 1
00328 #define Quick_max 14
00329 #define Int_max 14
00330 #ifndef NO_IEEE_Scale
00331 #define Avoid_Underflow
00332 #ifdef Flush_Denorm
00333 #undef Sudden_Underflow
00334 #endif
00335 #endif
00336
00337 #ifndef Flt_Rounds
00338 #ifdef FLT_ROUNDS
00339 #define Flt_Rounds FLT_ROUNDS
00340 #else
00341 #define Flt_Rounds 1
00342 #endif
00343 #endif
00344
00345 #ifdef Honor_FLT_ROUNDS
00346 #define Rounding rounding
00347 #undef Check_FLT_ROUNDS
00348 #define Check_FLT_ROUNDS
00349 #else
00350 #define Rounding Flt_Rounds
00351 #endif
00352
00353 #else
00354 #undef Check_FLT_ROUNDS
00355 #undef Honor_FLT_ROUNDS
00356 #undef SET_INEXACT
00357 #undef Sudden_Underflow
00358 #define Sudden_Underflow
00359 #ifdef IBM
00360 #undef Flt_Rounds
00361 #define Flt_Rounds 0
00362 #define Exp_shift 24
00363 #define Exp_shift1 24
00364 #define Exp_msk1 0x1000000
00365 #define Exp_msk11 0x1000000
00366 #define Exp_mask 0x7f000000
00367 #define P 14
00368 #define Bias 65
00369 #define Exp_1 0x41000000
00370 #define Exp_11 0x41000000
00371 #define Ebits 8
00372 #define Frac_mask 0xffffff
00373 #define Frac_mask1 0xffffff
00374 #define Bletch 4
00375 #define Ten_pmax 22
00376 #define Bndry_mask 0xefffff
00377 #define Bndry_mask1 0xffffff
00378 #define LSB 1
00379 #define Sign_bit 0x80000000
00380 #define Log2P 4
00381 #define Tiny0 0x100000
00382 #define Tiny1 0
00383 #define Quick_max 14
00384 #define Int_max 15
00385 #else
00386 #undef Flt_Rounds
00387 #define Flt_Rounds 1
00388 #define Exp_shift 23
00389 #define Exp_shift1 7
00390 #define Exp_msk1 0x80
00391 #define Exp_msk11 0x800000
00392 #define Exp_mask 0x7f80
00393 #define P 56
00394 #define Bias 129
00395 #define Exp_1 0x40800000
00396 #define Exp_11 0x4080
00397 #define Ebits 8
00398 #define Frac_mask 0x7fffff
00399 #define Frac_mask1 0xffff007f
00400 #define Ten_pmax 24
00401 #define Bletch 2
00402 #define Bndry_mask 0xffff007f
00403 #define Bndry_mask1 0xffff007f
00404 #define LSB 0x10000
00405 #define Sign_bit 0x8000
00406 #define Log2P 1
00407 #define Tiny0 0x80
00408 #define Tiny1 0
00409 #define Quick_max 15
00410 #define Int_max 15
00411 #endif
00412 #endif
00413
00414 #ifndef IEEE_Arith
00415 #define ROUND_BIASED
00416 #endif
00417
00418 #ifdef RND_PRODQUOT
00419 #define rounded_product(a,b) a = rnd_prod(a, b)
00420 #define rounded_quotient(a,b) a = rnd_quot(a, b)
00421 extern double rnd_prod(double, double), rnd_quot(double, double);
00422 #else
00423 #define rounded_product(a,b) a *= b
00424 #define rounded_quotient(a,b) a /= b
00425 #endif
00426
00427 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
00428 #define Big1 0xffffffff
00429
00430 #ifndef Pack_32
00431 #define Pack_32
00432 #endif
00433
00434 #define FFFFFFFF 0xffffffffUL
00435
00436 #ifdef NO_LONG_LONG
00437 #undef ULLong
00438 #ifdef Just_16
00439 #undef Pack_32
00440
00441
00442
00443
00444
00445 #endif
00446 #else
00447 #ifndef Llong
00448 #define Llong long long
00449 #endif
00450 #ifndef ULLong
00451 #define ULLong unsigned Llong
00452 #endif
00453 #endif
00454
00455 #ifndef MULTIPLE_THREADS
00456 #define ACQUIRE_DTOA_LOCK(n)
00457 #define FREE_DTOA_LOCK(n)
00458 #endif
00459
00460 #define Kmax (sizeof(size_t) << 3)
00461
00462 struct
00463 Bigint {
00464 struct Bigint *next;
00465 int k, maxwds, sign, wds;
00466 ULong x[1];
00467 };
00468
00469 typedef struct Bigint Bigint;
00470
00471 static Bigint *freelist[Kmax+1];
00472
00473 static Bigint *
00474 Balloc
00475 (int k)
00476 {
00477 int x;
00478 Bigint *rv;
00479 #ifndef Omit_Private_Memory
00480 unsigned int len;
00481 #endif
00482
00483 ACQUIRE_DTOA_LOCK(0);
00484 if ((rv = freelist[k])) {
00485 freelist[k] = rv->next;
00486 }
00487 else {
00488 x = 1 << k;
00489 #ifdef Omit_Private_Memory
00490 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
00491 #else
00492 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
00493 /sizeof(double);
00494 if (pmem_next - private_mem + len <= (unsigned)PRIVATE_mem) {
00495 rv = (Bigint*)pmem_next;
00496 pmem_next += len;
00497 }
00498 else
00499 rv = (Bigint*)MALLOC(len*sizeof(double));
00500 #endif
00501 rv->k = k;
00502 rv->maxwds = x;
00503 }
00504 FREE_DTOA_LOCK(0);
00505 rv->sign = rv->wds = 0;
00506 return rv;
00507 }
00508
00509 static void
00510 Bfree
00511 (Bigint *v)
00512 {
00513 if (v) {
00514 ACQUIRE_DTOA_LOCK(0);
00515 v->next = freelist[v->k];
00516 freelist[v->k] = v;
00517 FREE_DTOA_LOCK(0);
00518 }
00519 }
00520
00521 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
00522 y->wds*sizeof(Long) + 2*sizeof(int))
00523
00524 static Bigint *
00525 multadd
00526 (Bigint *b, int m, int a)
00527 {
00528 int i, wds;
00529 #ifdef ULLong
00530 ULong *x;
00531 ULLong carry, y;
00532 #else
00533 ULong carry, *x, y;
00534 #ifdef Pack_32
00535 ULong xi, z;
00536 #endif
00537 #endif
00538 Bigint *b1;
00539
00540 wds = b->wds;
00541 x = b->x;
00542 i = 0;
00543 carry = a;
00544 do {
00545 #ifdef ULLong
00546 y = *x * (ULLong)m + carry;
00547 carry = y >> 32;
00548 *x++ = (ULong)y & FFFFFFFF;
00549 #else
00550 #ifdef Pack_32
00551 xi = *x;
00552 y = (xi & 0xffff) * m + carry;
00553 z = (xi >> 16) * m + (y >> 16);
00554 carry = z >> 16;
00555 *x++ = (z << 16) + (y & 0xffff);
00556 #else
00557 y = *x * m + carry;
00558 carry = y >> 16;
00559 *x++ = y & 0xffff;
00560 #endif
00561 #endif
00562 }
00563 while(++i < wds);
00564 if (carry) {
00565 if (wds >= b->maxwds) {
00566 b1 = Balloc(b->k+1);
00567 Bcopy(b1, b);
00568 Bfree(b);
00569 b = b1;
00570 }
00571 b->x[wds++] = (ULong)carry;
00572 b->wds = wds;
00573 }
00574 return b;
00575 }
00576
00577 static Bigint *
00578 s2b
00579 (CONST char *s, int nd0, int nd, ULong y9)
00580 {
00581 Bigint *b;
00582 int i, k;
00583 Long x, y;
00584
00585 x = (nd + 8) / 9;
00586 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
00587 #ifdef Pack_32
00588 b = Balloc(k);
00589 b->x[0] = y9;
00590 b->wds = 1;
00591 #else
00592 b = Balloc(k+1);
00593 b->x[0] = y9 & 0xffff;
00594 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
00595 #endif
00596
00597 i = 9;
00598 if (9 < nd0) {
00599 s += 9;
00600 do b = multadd(b, 10, *s++ - '0');
00601 while(++i < nd0);
00602 s++;
00603 }
00604 else
00605 s += 10;
00606 for(; i < nd; i++)
00607 b = multadd(b, 10, *s++ - '0');
00608 return b;
00609 }
00610
00611 static int
00612 hi0bits
00613 (register ULong x)
00614 {
00615 register int k = 0;
00616
00617 if (!(x & 0xffff0000)) {
00618 k = 16;
00619 x <<= 16;
00620 }
00621 if (!(x & 0xff000000)) {
00622 k += 8;
00623 x <<= 8;
00624 }
00625 if (!(x & 0xf0000000)) {
00626 k += 4;
00627 x <<= 4;
00628 }
00629 if (!(x & 0xc0000000)) {
00630 k += 2;
00631 x <<= 2;
00632 }
00633 if (!(x & 0x80000000)) {
00634 k++;
00635 if (!(x & 0x40000000))
00636 return 32;
00637 }
00638 return k;
00639 }
00640
00641 static int
00642 lo0bits
00643 (ULong *y)
00644 {
00645 register int k;
00646 register ULong x = *y;
00647
00648 if (x & 7) {
00649 if (x & 1)
00650 return 0;
00651 if (x & 2) {
00652 *y = x >> 1;
00653 return 1;
00654 }
00655 *y = x >> 2;
00656 return 2;
00657 }
00658 k = 0;
00659 if (!(x & 0xffff)) {
00660 k = 16;
00661 x >>= 16;
00662 }
00663 if (!(x & 0xff)) {
00664 k += 8;
00665 x >>= 8;
00666 }
00667 if (!(x & 0xf)) {
00668 k += 4;
00669 x >>= 4;
00670 }
00671 if (!(x & 0x3)) {
00672 k += 2;
00673 x >>= 2;
00674 }
00675 if (!(x & 1)) {
00676 k++;
00677 x >>= 1;
00678 if (!x & 1)
00679 return 32;
00680 }
00681 *y = x;
00682 return k;
00683 }
00684
00685 static Bigint *
00686 i2b
00687 (int i)
00688 {
00689 Bigint *b;
00690
00691 b = Balloc(1);
00692 b->x[0] = i;
00693 b->wds = 1;
00694 return b;
00695 }
00696
00697 static Bigint *
00698 mult
00699 (Bigint *a, Bigint *b)
00700 {
00701 Bigint *c;
00702 int k, wa, wb, wc;
00703 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
00704 ULong y;
00705 #ifdef ULLong
00706 ULLong carry, z;
00707 #else
00708 ULong carry, z;
00709 #ifdef Pack_32
00710 ULong z2;
00711 #endif
00712 #endif
00713
00714 if (a->wds < b->wds) {
00715 c = a;
00716 a = b;
00717 b = c;
00718 }
00719 k = a->k;
00720 wa = a->wds;
00721 wb = b->wds;
00722 wc = wa + wb;
00723 if (wc > a->maxwds)
00724 k++;
00725 c = Balloc(k);
00726 for(x = c->x, xa = x + wc; x < xa; x++)
00727 *x = 0;
00728 xa = a->x;
00729 xae = xa + wa;
00730 xb = b->x;
00731 xbe = xb + wb;
00732 xc0 = c->x;
00733 #ifdef ULLong
00734 for(; xb < xbe; xc0++) {
00735 if ((y = *xb++)) {
00736 x = xa;
00737 xc = xc0;
00738 carry = 0;
00739 do {
00740 z = *x++ * (ULLong)y + *xc + carry;
00741 carry = z >> 32;
00742 *xc++ = (ULong)z & FFFFFFFF;
00743 }
00744 while(x < xae);
00745 *xc = (ULong)carry;
00746 }
00747 }
00748 #else
00749 #ifdef Pack_32
00750 for(; xb < xbe; xb++, xc0++) {
00751 if (y = *xb & 0xffff) {
00752 x = xa;
00753 xc = xc0;
00754 carry = 0;
00755 do {
00756 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
00757 carry = z >> 16;
00758 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
00759 carry = z2 >> 16;
00760 Storeinc(xc, z2, z);
00761 }
00762 while(x < xae);
00763 *xc = carry;
00764 }
00765 if (y = *xb >> 16) {
00766 x = xa;
00767 xc = xc0;
00768 carry = 0;
00769 z2 = *xc;
00770 do {
00771 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
00772 carry = z >> 16;
00773 Storeinc(xc, z, z2);
00774 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
00775 carry = z2 >> 16;
00776 }
00777 while(x < xae);
00778 *xc = z2;
00779 }
00780 }
00781 #else
00782 for(; xb < xbe; xc0++) {
00783 if (y = *xb++) {
00784 x = xa;
00785 xc = xc0;
00786 carry = 0;
00787 do {
00788 z = *x++ * y + *xc + carry;
00789 carry = z >> 16;
00790 *xc++ = z & 0xffff;
00791 }
00792 while(x < xae);
00793 *xc = carry;
00794 }
00795 }
00796 #endif
00797 #endif
00798 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
00799 c->wds = wc;
00800 return c;
00801 }
00802
00803 static Bigint *p5s;
00804
00805 static Bigint *
00806 pow5mult
00807 (Bigint *b, int k)
00808 {
00809 Bigint *b1, *p5, *p51;
00810 int i;
00811 static int p05[3] = { 5, 25, 125 };
00812
00813 if ((i = k & 3))
00814 b = multadd(b, p05[i-1], 0);
00815
00816 if (!(k >>= 2))
00817 return b;
00818 if (!(p5 = p5s)) {
00819
00820 #ifdef MULTIPLE_THREADS
00821 ACQUIRE_DTOA_LOCK(1);
00822 if (!(p5 = p5s)) {
00823 p5 = p5s = i2b(625);
00824 p5->next = 0;
00825 }
00826 FREE_DTOA_LOCK(1);
00827 #else
00828 p5 = p5s = i2b(625);
00829 p5->next = 0;
00830 #endif
00831 }
00832 for(;;) {
00833 if (k & 1) {
00834 b1 = mult(b, p5);
00835 Bfree(b);
00836 b = b1;
00837 }
00838 if (!(k >>= 1))
00839 break;
00840 if (!(p51 = p5->next)) {
00841 #ifdef MULTIPLE_THREADS
00842 ACQUIRE_DTOA_LOCK(1);
00843 if (!(p51 = p5->next)) {
00844 p51 = p5->next = mult(p5,p5);
00845 p51->next = 0;
00846 }
00847 FREE_DTOA_LOCK(1);
00848 #else
00849 p51 = p5->next = mult(p5,p5);
00850 p51->next = 0;
00851 #endif
00852 }
00853 p5 = p51;
00854 }
00855 return b;
00856 }
00857
00858 static Bigint *
00859 lshift
00860 (Bigint *b, int k)
00861 {
00862 int i, k1, n, n1;
00863 Bigint *b1;
00864 ULong *x, *x1, *xe, z;
00865
00866 #ifdef Pack_32
00867 n = k >> 5;
00868 #else
00869 n = k >> 4;
00870 #endif
00871 k1 = b->k;
00872 n1 = n + b->wds + 1;
00873 for(i = b->maxwds; n1 > i; i <<= 1)
00874 k1++;
00875 b1 = Balloc(k1);
00876 x1 = b1->x;
00877 for(i = 0; i < n; i++)
00878 *x1++ = 0;
00879 x = b->x;
00880 xe = x + b->wds;
00881 #ifdef Pack_32
00882 if (k &= 0x1f) {
00883 k1 = 32 - k;
00884 z = 0;
00885 do {
00886 *x1++ = *x << k | z;
00887 z = *x++ >> k1;
00888 }
00889 while(x < xe);
00890 if ((*x1 = z))
00891 ++n1;
00892 }
00893 #else
00894 if (k &= 0xf) {
00895 k1 = 16 - k;
00896 z = 0;
00897 do {
00898 *x1++ = *x << k & 0xffff | z;
00899 z = *x++ >> k1;
00900 }
00901 while(x < xe);
00902 if (*x1 = z)
00903 ++n1;
00904 }
00905 #endif
00906 else do
00907 *x1++ = *x++;
00908 while(x < xe);
00909 b1->wds = n1 - 1;
00910 Bfree(b);
00911 return b1;
00912 }
00913
00914 static int
00915 cmp
00916 (Bigint *a, Bigint *b)
00917 {
00918 ULong *xa, *xa0, *xb, *xb0;
00919 int i, j;
00920
00921 i = a->wds;
00922 j = b->wds;
00923 #ifdef DEBUG
00924 if (i > 1 && !a->x[i-1])
00925 Bug("cmp called with a->x[a->wds-1] == 0");
00926 if (j > 1 && !b->x[j-1])
00927 Bug("cmp called with b->x[b->wds-1] == 0");
00928 #endif
00929 if (i -= j)
00930 return i;
00931 xa0 = a->x;
00932 xa = xa0 + j;
00933 xb0 = b->x;
00934 xb = xb0 + j;
00935 for(;;) {
00936 if (*--xa != *--xb)
00937 return *xa < *xb ? -1 : 1;
00938 if (xa <= xa0)
00939 break;
00940 }
00941 return 0;
00942 }
00943
00944 static Bigint *
00945 diff
00946 (Bigint *a, Bigint *b)
00947 {
00948 Bigint *c;
00949 int i, wa, wb;
00950 ULong *xa, *xae, *xb, *xbe, *xc;
00951 #ifdef ULLong
00952 ULLong borrow, y;
00953 #else
00954 ULong borrow, y;
00955 #ifdef Pack_32
00956 ULong z;
00957 #endif
00958 #endif
00959
00960 i = cmp(a,b);
00961 if (!i) {
00962 c = Balloc(0);
00963 c->wds = 1;
00964 c->x[0] = 0;
00965 return c;
00966 }
00967 if (i < 0) {
00968 c = a;
00969 a = b;
00970 b = c;
00971 i = 1;
00972 }
00973 else
00974 i = 0;
00975 c = Balloc(a->k);
00976 c->sign = i;
00977 wa = a->wds;
00978 xa = a->x;
00979 xae = xa + wa;
00980 wb = b->wds;
00981 xb = b->x;
00982 xbe = xb + wb;
00983 xc = c->x;
00984 borrow = 0;
00985 #ifdef ULLong
00986 do {
00987 y = (ULLong)*xa++ - *xb++ - borrow;
00988 borrow = y >> 32 & (ULong)1;
00989 *xc++ = (ULong)y & FFFFFFFF;
00990 }
00991 while(xb < xbe);
00992 while(xa < xae) {
00993 y = *xa++ - borrow;
00994 borrow = y >> 32 & (ULong)1;
00995 *xc++ = (ULong)y & FFFFFFFF;
00996 }
00997 #else
00998 #ifdef Pack_32
00999 do {
01000 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
01001 borrow = (y & 0x10000) >> 16;
01002 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
01003 borrow = (z & 0x10000) >> 16;
01004 Storeinc(xc, z, y);
01005 }
01006 while(xb < xbe);
01007 while(xa < xae) {
01008 y = (*xa & 0xffff) - borrow;
01009 borrow = (y & 0x10000) >> 16;
01010 z = (*xa++ >> 16) - borrow;
01011 borrow = (z & 0x10000) >> 16;
01012 Storeinc(xc, z, y);
01013 }
01014 #else
01015 do {
01016 y = *xa++ - *xb++ - borrow;
01017 borrow = (y & 0x10000) >> 16;
01018 *xc++ = y & 0xffff;
01019 }
01020 while(xb < xbe);
01021 while(xa < xae) {
01022 y = *xa++ - borrow;
01023 borrow = (y & 0x10000) >> 16;
01024 *xc++ = y & 0xffff;
01025 }
01026 #endif
01027 #endif
01028 while(!*--xc)
01029 wa--;
01030 c->wds = wa;
01031 return c;
01032 }
01033
01034 static double
01035 ulp
01036 (double dx)
01037 {
01038 register Long L;
01039 U x, a;
01040
01041 dval(x) = dx;
01042 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
01043 #ifndef Avoid_Underflow
01044 #ifndef Sudden_Underflow
01045 if (L > 0) {
01046 #endif
01047 #endif
01048 #ifdef IBM
01049 L |= Exp_msk1 >> 4;
01050 #endif
01051 word0(a) = L;
01052 word1(a) = 0;
01053 #ifndef Avoid_Underflow
01054 #ifndef Sudden_Underflow
01055 }
01056 else {
01057 L = -L >> Exp_shift;
01058 if (L < Exp_shift) {
01059 word0(a) = 0x80000 >> L;
01060 word1(a) = 0;
01061 }
01062 else {
01063 word0(a) = 0;
01064 L -= Exp_shift;
01065 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
01066 }
01067 }
01068 #endif
01069 #endif
01070 return dval(a);
01071 }
01072
01073 static double
01074 b2d
01075 (Bigint *a, int *e)
01076 {
01077 ULong *xa, *xa0, w, y, z;
01078 int k;
01079 U d;
01080 #ifdef VAX
01081 ULong d0, d1;
01082 #else
01083 #define d0 word0(d)
01084 #define d1 word1(d)
01085 #endif
01086
01087 xa0 = a->x;
01088 xa = xa0 + a->wds;
01089 y = *--xa;
01090 #ifdef DEBUG
01091 if (!y) Bug("zero y in b2d");
01092 #endif
01093 k = hi0bits(y);
01094 *e = 32 - k;
01095 #ifdef Pack_32
01096 if (k < Ebits) {
01097 d0 = Exp_1 | y >> Ebits - k;
01098 w = xa > xa0 ? *--xa : 0;
01099 d1 = y << (32-Ebits) + k | w >> Ebits - k;
01100 goto ret_d;
01101 }
01102 z = xa > xa0 ? *--xa : 0;
01103 if (k -= Ebits) {
01104 d0 = Exp_1 | y << k | z >> 32 - k;
01105 y = xa > xa0 ? *--xa : 0;
01106 d1 = z << k | y >> 32 - k;
01107 }
01108 else {
01109 d0 = Exp_1 | y;
01110 d1 = z;
01111 }
01112 #else
01113 if (k < Ebits + 16) {
01114 z = xa > xa0 ? *--xa : 0;
01115 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
01116 w = xa > xa0 ? *--xa : 0;
01117 y = xa > xa0 ? *--xa : 0;
01118 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
01119 goto ret_d;
01120 }
01121 z = xa > xa0 ? *--xa : 0;
01122 w = xa > xa0 ? *--xa : 0;
01123 k -= Ebits + 16;
01124 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
01125 y = xa > xa0 ? *--xa : 0;
01126 d1 = w << k + 16 | y << k;
01127 #endif
01128 ret_d:
01129 #ifdef VAX
01130 word0(d) = d0 >> 16 | d0 << 16;
01131 word1(d) = d1 >> 16 | d1 << 16;
01132 #else
01133 #undef d0
01134 #undef d1
01135 #endif
01136 return dval(d);
01137 }
01138
01139 static Bigint *
01140 d2b
01141 (double dd, int *e, int *bits)
01142 {
01143 U d;
01144 Bigint *b;
01145 int de, k;
01146 ULong *x, y, z;
01147 #ifndef Sudden_Underflow
01148 int i;
01149 #endif
01150 #ifdef VAX
01151 ULong d0, d1;
01152 #endif
01153 dval(d) = dd;
01154 #ifdef VAX
01155 d0 = word0(d) >> 16 | word0(d) << 16;
01156 d1 = word1(d) >> 16 | word1(d) << 16;
01157 #else
01158 #define d0 word0(d)
01159 #define d1 word1(d)
01160 #endif
01161
01162 #ifdef Pack_32
01163 b = Balloc(1);
01164 #else
01165 b = Balloc(2);
01166 #endif
01167 x = b->x;
01168
01169 z = d0 & Frac_mask;
01170 d0 &= 0x7fffffff;
01171 #ifdef Sudden_Underflow
01172 de = (int)(d0 >> Exp_shift);
01173 #ifndef IBM
01174 z |= Exp_msk11;
01175 #endif
01176 #else
01177 if ((de = (int)(d0 >> Exp_shift)))
01178 z |= Exp_msk1;
01179 #endif
01180 #ifdef Pack_32
01181 if ((y = d1)) {
01182 if ((k = lo0bits(&y))) {
01183 x[0] = y | z << 32 - k;
01184 z >>= k;
01185 }
01186 else
01187 x[0] = y;
01188 #ifndef Sudden_Underflow
01189 i =
01190 #endif
01191 b->wds = (x[1] = z) ? 2 : 1;
01192 }
01193 else {
01194 #ifdef DEBUG
01195 if (!z)
01196 Bug("Zero passed to d2b");
01197 #endif
01198 k = lo0bits(&z);
01199 x[0] = z;
01200 #ifndef Sudden_Underflow
01201 i =
01202 #endif
01203 b->wds = 1;
01204 k += 32;
01205 }
01206 #else
01207 if (y = d1) {
01208 if (k = lo0bits(&y))
01209 if (k >= 16) {
01210 x[0] = y | z << 32 - k & 0xffff;
01211 x[1] = z >> k - 16 & 0xffff;
01212 x[2] = z >> k;
01213 i = 2;
01214 }
01215 else {
01216 x[0] = y & 0xffff;
01217 x[1] = y >> 16 | z << 16 - k & 0xffff;
01218 x[2] = z >> k & 0xffff;
01219 x[3] = z >> k+16;
01220 i = 3;
01221 }
01222 else {
01223 x[0] = y & 0xffff;
01224 x[1] = y >> 16;
01225 x[2] = z & 0xffff;
01226 x[3] = z >> 16;
01227 i = 3;
01228 }
01229 }
01230 else {
01231 #ifdef DEBUG
01232 if (!z)
01233 Bug("Zero passed to d2b");
01234 #endif
01235 k = lo0bits(&z);
01236 if (k >= 16) {
01237 x[0] = z;
01238 i = 0;
01239 }
01240 else {
01241 x[0] = z & 0xffff;
01242 x[1] = z >> 16;
01243 i = 1;
01244 }
01245 k += 32;
01246 }
01247 while(!x[i])
01248 --i;
01249 b->wds = i + 1;
01250 #endif
01251 #ifndef Sudden_Underflow
01252 if (de) {
01253 #endif
01254 #ifdef IBM
01255 *e = (de - Bias - (P-1) << 2) + k;
01256 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
01257 #else
01258 *e = de - Bias - (P-1) + k;
01259 *bits = P - k;
01260 #endif
01261 #ifndef Sudden_Underflow
01262 }
01263 else {
01264 *e = de - Bias - (P-1) + 1 + k;
01265 #ifdef Pack_32
01266 *bits = 32*i - hi0bits(x[i-1]);
01267 #else
01268 *bits = (i+2)*16 - hi0bits(x[i]);
01269 #endif
01270 }
01271 #endif
01272 return b;
01273 }
01274 #undef d0
01275 #undef d1
01276
01277 static double
01278 ratio
01279 (Bigint *a, Bigint *b)
01280 {
01281 U da, db;
01282 int k, ka, kb;
01283
01284 dval(da) = b2d(a, &ka);
01285 dval(db) = b2d(b, &kb);
01286 #ifdef Pack_32
01287 k = ka - kb + 32*(a->wds - b->wds);
01288 #else
01289 k = ka - kb + 16*(a->wds - b->wds);
01290 #endif
01291 #ifdef IBM
01292 if (k > 0) {
01293 word0(da) += (k >> 2)*Exp_msk1;
01294 if (k &= 3)
01295 dval(da) *= 1 << k;
01296 }
01297 else {
01298 k = -k;
01299 word0(db) += (k >> 2)*Exp_msk1;
01300 if (k &= 3)
01301 dval(db) *= 1 << k;
01302 }
01303 #else
01304 if (k > 0)
01305 word0(da) += k*Exp_msk1;
01306 else {
01307 k = -k;
01308 word0(db) += k*Exp_msk1;
01309 }
01310 #endif
01311 return dval(da) / dval(db);
01312 }
01313
01314 static CONST double
01315 tens[] = {
01316 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
01317 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
01318 1e20, 1e21, 1e22
01319 #ifdef VAX
01320 , 1e23, 1e24
01321 #endif
01322 };
01323
01324 static CONST double
01325 #ifdef IEEE_Arith
01326 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
01327 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
01328 #ifdef Avoid_Underflow
01329 9007199254740992.*9007199254740992.e-256
01330
01331 #else
01332 1e-256
01333 #endif
01334 };
01335
01336
01337 #define Scale_Bit 0x10
01338 #define n_bigtens 5
01339 #else
01340 #ifdef IBM
01341 bigtens[] = { 1e16, 1e32, 1e64 };
01342 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
01343 #define n_bigtens 3
01344 #else
01345 bigtens[] = { 1e16, 1e32 };
01346 static CONST double tinytens[] = { 1e-16, 1e-32 };
01347 #define n_bigtens 2
01348 #endif
01349 #endif
01350
01351 #ifndef IEEE_Arith
01352 #undef INFNAN_CHECK
01353 #endif
01354
01355 #ifdef INFNAN_CHECK
01356
01357 #ifndef NAN_WORD0
01358 #define NAN_WORD0 0x7ff80000
01359 #endif
01360
01361 #ifndef NAN_WORD1
01362 #define NAN_WORD1 0
01363 #endif
01364
01365 static int
01366 match
01367 (CONST char **sp, CONST char *t)
01368 {
01369 int c, d;
01370 CONST char *s = *sp;
01371
01372 while((d = *t++)) {
01373 if ((c = *++s) >= 'A' && c <= 'Z')
01374 c += 'a' - 'A';
01375 if (c != d)
01376 return 0;
01377 }
01378 *sp = s + 1;
01379 return 1;
01380 }
01381
01382 #ifndef No_Hex_NaN
01383 static void
01384 hexnan
01385 (U *rvp, CONST char **sp)
01386 {
01387 ULong c, x[2];
01388 CONST char *s;
01389 int havedig, udx0, xshift;
01390
01391 x[0] = x[1] = 0;
01392 havedig = xshift = 0;
01393 udx0 = 1;
01394 s = *sp;
01395 while((c = *(CONST unsigned char*)++s)) {
01396 if (c >= '0' && c <= '9')
01397 c -= '0';
01398 else if (c >= 'a' && c <= 'f')
01399 c += 10 - 'a';
01400 else if (c >= 'A' && c <= 'F')
01401 c += 10 - 'A';
01402 else if (c <= ' ') {
01403 if (udx0 && havedig) {
01404 udx0 = 0;
01405 xshift = 1;
01406 }
01407 continue;
01408 }
01409 else if ( c == ')' && havedig) {
01410 *sp = s + 1;
01411 break;
01412 }
01413 else
01414 return;
01415 havedig = 1;
01416 if (xshift) {
01417 xshift = 0;
01418 x[0] = x[1];
01419 x[1] = 0;
01420 }
01421 if (udx0)
01422 x[0] = (x[0] << 4) | (x[1] >> 28);
01423 x[1] = (x[1] << 4) | c;
01424 }
01425 if ((x[0] &= 0xfffff) || x[1]) {
01426 word0(*rvp) = Exp_mask | x[0];
01427 word1(*rvp) = x[1];
01428 }
01429 }
01430 #endif
01431 #endif
01432
01433 double
01434 strtod
01435 (CONST char *s00, char **se)
01436 {
01437 #ifdef Avoid_Underflow
01438 int scale;
01439 #endif
01440 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
01441 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
01442 CONST char *s, *s0, *s1;
01443 double aadj, aadj1, adj;
01444 U aadj2, rv, rv0;
01445 Long L;
01446 ULong y, z;
01447 Bigint *bb = NULL, *bb1 = NULL, *bd = NULL, *bd0 = NULL, *bs = NULL, *delta = NULL;
01448 #ifdef SET_INEXACT
01449 int inexact, oldinexact;
01450 #endif
01451 #ifdef Honor_FLT_ROUNDS
01452 int rounding;
01453 #endif
01454 #ifdef USE_LOCALE
01455 CONST char *s2;
01456 #endif
01457
01458 sign = nz0 = nz = 0;
01459 dval(rv) = 0.;
01460 for(s = s00;;s++) switch(*s) {
01461 case '-':
01462 sign = 1;
01463
01464 case '+':
01465 if (*++s)
01466 goto break2;
01467
01468 case 0:
01469 goto ret0;
01470 case '\t':
01471 case '\n':
01472 case '\v':
01473 case '\f':
01474 case '\r':
01475 case ' ':
01476 continue;
01477 default:
01478 goto break2;
01479 }
01480 break2:
01481 if (*s == '0') {
01482 nz0 = 1;
01483 while(*++s == '0') ;
01484 if (!*s)
01485 goto ret;
01486 }
01487 s0 = s;
01488 y = z = 0;
01489 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
01490 if (nd < 9)
01491 y = 10*y + c - '0';
01492 else if (nd < 16)
01493 z = 10*z + c - '0';
01494 nd0 = nd;
01495 #ifdef USE_LOCALE
01496 s1 = localeconv()->decimal_point;
01497 if (c == *s1) {
01498 c = '.';
01499 if (*++s1) {
01500 s2 = s;
01501 for(;;) {
01502 if (*++s2 != *s1) {
01503 c = 0;
01504 break;
01505 }
01506 if (!*++s1) {
01507 s = s2;
01508 break;
01509 }
01510 }
01511 }
01512 }
01513 #endif
01514 if (c == '.') {
01515 c = *++s;
01516 if (!nd) {
01517 for(; c == '0'; c = *++s)
01518 nz++;
01519 if (c > '0' && c <= '9') {
01520 s0 = s;
01521 nf += nz;
01522 nz = 0;
01523 goto have_dig;
01524 }
01525 goto dig_done;
01526 }
01527 for(; c >= '0' && c <= '9'; c = *++s) {
01528 have_dig:
01529 nz++;
01530 if (c -= '0') {
01531 nf += nz;
01532 for(i = 1; i < nz; i++)
01533 if (nd++ < 9)
01534 y *= 10;
01535 else if (nd <= DBL_DIG + 1)
01536 z *= 10;
01537 if (nd++ < 9)
01538 y = 10*y + c;
01539 else if (nd <= DBL_DIG + 1)
01540 z = 10*z + c;
01541 nz = 0;
01542 }
01543 }
01544 }
01545 dig_done:
01546 e = 0;
01547 if (c == 'e' || c == 'E') {
01548 if (!nd && !nz && !nz0) {
01549 goto ret0;
01550 }
01551 s00 = s;
01552 esign = 0;
01553 switch(c = *++s) {
01554 case '-':
01555 esign = 1;
01556 case '+':
01557 c = *++s;
01558 }
01559 if (c >= '0' && c <= '9') {
01560 while(c == '0')
01561 c = *++s;
01562 if (c > '0' && c <= '9') {
01563 L = c - '0';
01564 s1 = s;
01565 while((c = *++s) >= '0' && c <= '9')
01566 L = 10*L + c - '0';
01567 if (s - s1 > 8 || L > 19999)
01568
01569
01570
01571 e = 19999;
01572 else
01573 e = (int)L;
01574 if (esign)
01575 e = -e;
01576 }
01577 else
01578 e = 0;
01579 }
01580 else
01581 s = s00;
01582 }
01583 if (!nd) {
01584 if (!nz && !nz0) {
01585 #ifdef INFNAN_CHECK
01586
01587 switch(c) {
01588 case 'i':
01589 case 'I':
01590 if (match(&s,"nf")) {
01591 --s;
01592 if (!match(&s,"inity"))
01593 ++s;
01594 word0(rv) = 0x7ff00000;
01595 word1(rv) = 0;
01596 goto ret;
01597 }
01598 break;
01599 case 'n':
01600 case 'N':
01601 if (match(&s, "an")) {
01602 word0(rv) = NAN_WORD0;
01603 word1(rv) = NAN_WORD1;
01604 #ifndef No_Hex_NaN
01605 if (*s == '(')
01606 hexnan(&rv, &s);
01607 #endif
01608 goto ret;
01609 }
01610 }
01611 #endif
01612 ret0:
01613 s = s00;
01614 sign = 0;
01615 }
01616 goto ret;
01617 }
01618 e1 = e -= nf;
01619
01620
01621
01622
01623
01624
01625 if (!nd0)
01626 nd0 = nd;
01627 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
01628 dval(rv) = y;
01629 if (k > 9) {
01630 #ifdef SET_INEXACT
01631 if (k > DBL_DIG)
01632 oldinexact = get_inexact();
01633 #endif
01634 dval(rv) = tens[k - 9] * dval(rv) + z;
01635 }
01636 bd0 = 0;
01637 if (nd <= DBL_DIG
01638 #ifndef RND_PRODQUOT
01639 #ifndef Honor_FLT_ROUNDS
01640 && Flt_Rounds == 1
01641 #endif
01642 #endif
01643 ) {
01644 if (!e)
01645 goto ret;
01646 if (e > 0) {
01647 if (e <= Ten_pmax) {
01648 #ifdef VAX
01649 goto vax_ovfl_check;
01650 #else
01651 #ifdef Honor_FLT_ROUNDS
01652
01653 if (sign) {
01654 rv = -rv;
01655 sign = 0;
01656 }
01657 #endif
01658 rounded_product(dval(rv), tens[e]);
01659 goto ret;
01660 #endif
01661 }
01662 i = DBL_DIG - nd;
01663 if (e <= Ten_pmax + i) {
01664
01665
01666
01667 #ifdef Honor_FLT_ROUNDS
01668
01669 if (sign) {
01670 rv = -rv;
01671 sign = 0;
01672 }
01673 #endif
01674 e -= i;
01675 dval(rv) *= tens[i];
01676 #ifdef VAX
01677
01678
01679
01680 vax_ovfl_check:
01681 word0(rv) -= P*Exp_msk1;
01682 rounded_product(dval(rv), tens[e]);
01683 if ((word0(rv) & Exp_mask)
01684 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
01685 goto ovfl;
01686 word0(rv) += P*Exp_msk1;
01687 #else
01688 rounded_product(dval(rv), tens[e]);
01689 #endif
01690 goto ret;
01691 }
01692 }
01693 #ifndef Inaccurate_Divide
01694 else if (e >= -Ten_pmax) {
01695 #ifdef Honor_FLT_ROUNDS
01696
01697 if (sign) {
01698 rv = -rv;
01699 sign = 0;
01700 }
01701 #endif
01702 rounded_quotient(dval(rv), tens[-e]);
01703 goto ret;
01704 }
01705 #endif
01706 }
01707 e1 += nd - k;
01708
01709 #ifdef IEEE_Arith
01710 #ifdef SET_INEXACT
01711 inexact = 1;
01712 if (k <= DBL_DIG)
01713 oldinexact = get_inexact();
01714 #endif
01715 #ifdef Avoid_Underflow
01716 scale = 0;
01717 #endif
01718 #ifdef Honor_FLT_ROUNDS
01719 if ((rounding = Flt_Rounds) >= 2) {
01720 if (sign)
01721 rounding = rounding == 2 ? 0 : 2;
01722 else
01723 if (rounding != 2)
01724 rounding = 0;
01725 }
01726 #endif
01727 #endif
01728
01729
01730
01731 if (e1 > 0) {
01732 if ((i = e1 & 15))
01733 dval(rv) *= tens[i];
01734 if (e1 &= ~15) {
01735 if (e1 > DBL_MAX_10_EXP) {
01736 ovfl:
01737 #ifndef NO_ERRNO
01738 errno = ERANGE;
01739 #endif
01740
01741 #ifdef IEEE_Arith
01742 #ifdef Honor_FLT_ROUNDS
01743 switch(rounding) {
01744 case 0:
01745 case 3:
01746 word0(rv) = Big0;
01747 word1(rv) = Big1;
01748 break;
01749 default:
01750 word0(rv) = Exp_mask;
01751 word1(rv) = 0;
01752 }
01753 #else
01754 word0(rv) = Exp_mask;
01755 word1(rv) = 0;
01756 #endif
01757 #ifdef SET_INEXACT
01758
01759 dval(rv0) = 1e300;
01760 dval(rv0) *= dval(rv0);
01761 #endif
01762 #else
01763 word0(rv) = Big0;
01764 word1(rv) = Big1;
01765 #endif
01766 if (bd0)
01767 goto retfree;
01768 goto ret;
01769 }
01770 e1 >>= 4;
01771 for(j = 0; e1 > 1; j++, e1 >>= 1)
01772 if (e1 & 1)
01773 dval(rv) *= bigtens[j];
01774
01775 word0(rv) -= P*Exp_msk1;
01776 dval(rv) *= bigtens[j];
01777 if ((z = word0(rv) & Exp_mask)
01778 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
01779 goto ovfl;
01780 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
01781
01782
01783 word0(rv) = Big0;
01784 word1(rv) = Big1;
01785 }
01786 else
01787 word0(rv) += P*Exp_msk1;
01788 }
01789 }
01790 else if (e1 < 0) {
01791 e1 = -e1;
01792 if ((i = e1 & 15))
01793 dval(rv) /= tens[i];
01794 if (e1 >>= 4) {
01795 if (e1 >= 1 << n_bigtens)
01796 goto undfl;
01797 #ifdef Avoid_Underflow
01798 if (e1 & Scale_Bit)
01799 scale = 2*P;
01800 for(j = 0; e1 > 0; j++, e1 >>= 1)
01801 if (e1 & 1)
01802 dval(rv) *= tinytens[j];
01803 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
01804 >> Exp_shift)) > 0) {
01805
01806 if (j >= 32) {
01807 word1(rv) = 0;
01808 if (j >= 53)
01809 word0(rv) = (P+2)*Exp_msk1;
01810 else
01811 word0(rv) &= 0xffffffff << j-32;
01812 }
01813 else
01814 word1(rv) &= 0xffffffff << j;
01815 }
01816 #else
01817 for(j = 0; e1 > 1; j++, e1 >>= 1)
01818 if (e1 & 1)
01819 dval(rv) *= tinytens[j];
01820
01821 dval(rv0) = dval(rv);
01822 dval(rv) *= tinytens[j];
01823 if (!dval(rv)) {
01824 dval(rv) = 2.*dval(rv0);
01825 dval(rv) *= tinytens[j];
01826 #endif
01827 if (!dval(rv)) {
01828 undfl:
01829 dval(rv) = 0.;
01830 #ifndef NO_ERRNO
01831 errno = ERANGE;
01832 #endif
01833 if (bd0)
01834 goto retfree;
01835 goto ret;
01836 }
01837 #ifndef Avoid_Underflow
01838 word0(rv) = Tiny0;
01839 word1(rv) = Tiny1;
01840
01841
01842
01843 }
01844 #endif
01845 }
01846 }
01847
01848
01849
01850
01851
01852 bd0 = s2b(s0, nd0, nd, y);
01853
01854 for(;;) {
01855 bd = Balloc(bd0->k);
01856 Bcopy(bd, bd0);
01857 bb = d2b(dval(rv), &bbe, &bbbits);
01858 bs = i2b(1);
01859
01860 if (e >= 0) {
01861 bb2 = bb5 = 0;
01862 bd2 = bd5 = e;
01863 }
01864 else {
01865 bb2 = bb5 = -e;
01866 bd2 = bd5 = 0;
01867 }
01868 if (bbe >= 0)
01869 bb2 += bbe;
01870 else
01871 bd2 -= bbe;
01872 bs2 = bb2;
01873 #ifdef Honor_FLT_ROUNDS
01874 if (rounding != 1)
01875 bs2++;
01876 #endif
01877 #ifdef Avoid_Underflow
01878 j = bbe - scale;
01879 i = j + bbbits - 1;
01880 if (i < Emin)
01881 j += P - Emin;
01882 else
01883 j = P + 1 - bbbits;
01884 #else
01885 #ifdef Sudden_Underflow
01886 #ifdef IBM
01887 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
01888 #else
01889 j = P + 1 - bbbits;
01890 #endif
01891 #else
01892 j = bbe;
01893 i = j + bbbits - 1;
01894 if (i < Emin)
01895 j += P - Emin;
01896 else
01897 j = P + 1 - bbbits;
01898 #endif
01899 #endif
01900 bb2 += j;
01901 bd2 += j;
01902 #ifdef Avoid_Underflow
01903 bd2 += scale;
01904 #endif
01905 i = bb2 < bd2 ? bb2 : bd2;
01906 if (i > bs2)
01907 i = bs2;
01908 if (i > 0) {
01909 bb2 -= i;
01910 bd2 -= i;
01911 bs2 -= i;
01912 }
01913 if (bb5 > 0) {
01914 bs = pow5mult(bs, bb5);
01915 bb1 = mult(bs, bb);
01916 Bfree(bb);
01917 bb = bb1;
01918 }
01919 if (bb2 > 0)
01920 bb = lshift(bb, bb2);
01921 if (bd5 > 0)
01922 bd = pow5mult(bd, bd5);
01923 if (bd2 > 0)
01924 bd = lshift(bd, bd2);
01925 if (bs2 > 0)
01926 bs = lshift(bs, bs2);
01927 delta = diff(bb, bd);
01928 dsign = delta->sign;
01929 delta->sign = 0;
01930 i = cmp(delta, bs);
01931 #ifdef Honor_FLT_ROUNDS
01932 if (rounding != 1) {
01933 if (i < 0) {
01934
01935 if (!delta->x[0] && delta->wds <= 1) {
01936
01937 #ifdef SET_INEXACT
01938 inexact = 0;
01939 #endif
01940 break;
01941 }
01942 if (rounding) {
01943 if (dsign) {
01944 adj = 1.;
01945 goto apply_adj;
01946 }
01947 }
01948 else if (!dsign) {
01949 adj = -1.;
01950 if (!word1(rv)
01951 && !(word0(rv) & Frac_mask)) {
01952 y = word0(rv) & Exp_mask;
01953 #ifdef Avoid_Underflow
01954 if (!scale || y > 2*P*Exp_msk1)
01955 #else
01956 if (y)
01957 #endif
01958 {
01959 delta = lshift(delta,Log2P);
01960 if (cmp(delta, bs) <= 0)
01961 adj = -0.5;
01962 }
01963 }
01964 apply_adj:
01965 #ifdef Avoid_Underflow
01966 if (scale && (y = word0(rv) & Exp_mask)
01967 <= 2*P*Exp_msk1)
01968 word0(adj) += (2*P+1)*Exp_msk1 - y;
01969 #else
01970 #ifdef Sudden_Underflow
01971 if ((word0(rv) & Exp_mask) <=
01972 P*Exp_msk1) {
01973 word0(rv) += P*Exp_msk1;
01974 dval(rv) += adj*ulp(dval(rv));
01975 word0(rv) -= P*Exp_msk1;
01976 }
01977 else
01978 #endif
01979 #endif
01980 dval(rv) += adj*ulp(dval(rv));
01981 }
01982 break;
01983 }
01984 adj = ratio(delta, bs);
01985 if (adj < 1.)
01986 adj = 1.;
01987 if (adj <= 0x7ffffffe) {
01988
01989 y = adj;
01990 if (y != adj) {
01991 if (!((rounding>>1) ^ dsign))
01992 y++;
01993 adj = y;
01994 }
01995 }
01996 #ifdef Avoid_Underflow
01997 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
01998 word0(adj) += (2*P+1)*Exp_msk1 - y;
01999 #else
02000 #ifdef Sudden_Underflow
02001 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02002 word0(rv) += P*Exp_msk1;
02003 adj *= ulp(dval(rv));
02004 if (dsign)
02005 dval(rv) += adj;
02006 else
02007 dval(rv) -= adj;
02008 word0(rv) -= P*Exp_msk1;
02009 goto cont;
02010 }
02011 #endif
02012 #endif
02013 adj *= ulp(dval(rv));
02014 if (dsign)
02015 dval(rv) += adj;
02016 else
02017 dval(rv) -= adj;
02018 goto cont;
02019 }
02020 #endif
02021
02022 if (i < 0) {
02023
02024
02025
02026 if (dsign || word1(rv) || word0(rv) & Bndry_mask
02027 #ifdef IEEE_Arith
02028 #ifdef Avoid_Underflow
02029 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
02030 #else
02031 || (word0(rv) & Exp_mask) <= Exp_msk1
02032 #endif
02033 #endif
02034 ) {
02035 #ifdef SET_INEXACT
02036 if (!delta->x[0] && delta->wds <= 1)
02037 inexact = 0;
02038 #endif
02039 break;
02040 }
02041 if (!delta->x[0] && delta->wds <= 1) {
02042
02043 #ifdef SET_INEXACT
02044 inexact = 0;
02045 #endif
02046 break;
02047 }
02048 delta = lshift(delta,Log2P);
02049 if (cmp(delta, bs) > 0)
02050 goto drop_down;
02051 break;
02052 }
02053 if (i == 0) {
02054
02055 if (dsign) {
02056 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
02057 && word1(rv) == (
02058 #ifdef Avoid_Underflow
02059 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02060 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
02061 #endif
02062 0xffffffff)) {
02063
02064 word0(rv) = (word0(rv) & Exp_mask)
02065 + Exp_msk1
02066 #ifdef IBM
02067 | Exp_msk1 >> 4
02068 #endif
02069 ;
02070 word1(rv) = 0;
02071 #ifdef Avoid_Underflow
02072 dsign = 0;
02073 #endif
02074 break;
02075 }
02076 }
02077 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
02078 drop_down:
02079
02080 #ifdef Sudden_Underflow
02081 L = word0(rv) & Exp_mask;
02082 #ifdef IBM
02083 if (L < Exp_msk1)
02084 #else
02085 #ifdef Avoid_Underflow
02086 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
02087 #else
02088 if (L <= Exp_msk1)
02089 #endif
02090 #endif
02091 goto undfl;
02092 L -= Exp_msk1;
02093 #else
02094 #ifdef Avoid_Underflow
02095 if (scale) {
02096 L = word0(rv) & Exp_mask;
02097 if (L <= (2*P+1)*Exp_msk1) {
02098 if (L > (P+2)*Exp_msk1)
02099
02100
02101 break;
02102
02103 goto undfl;
02104 }
02105 }
02106 #endif
02107 L = (word0(rv) & Exp_mask) - Exp_msk1;
02108 #endif
02109 word0(rv) = L | Bndry_mask1;
02110 word1(rv) = 0xffffffff;
02111 #ifdef IBM
02112 goto cont;
02113 #else
02114 break;
02115 #endif
02116 }
02117 #ifndef ROUND_BIASED
02118 if (!(word1(rv) & LSB))
02119 break;
02120 #endif
02121 if (dsign)
02122 dval(rv) += ulp(dval(rv));
02123 #ifndef ROUND_BIASED
02124 else {
02125 dval(rv) -= ulp(dval(rv));
02126 #ifndef Sudden_Underflow
02127 if (!dval(rv))
02128 goto undfl;
02129 #endif
02130 }
02131 #ifdef Avoid_Underflow
02132 dsign = 1 - dsign;
02133 #endif
02134 #endif
02135 break;
02136 }
02137 if ((aadj = ratio(delta, bs)) <= 2.) {
02138 if (dsign)
02139 aadj = aadj1 = 1.;
02140 else if (word1(rv) || word0(rv) & Bndry_mask) {
02141 #ifndef Sudden_Underflow
02142 if (word1(rv) == Tiny1 && !word0(rv))
02143 goto undfl;
02144 #endif
02145 aadj = 1.;
02146 aadj1 = -1.;
02147 }
02148 else {
02149
02150
02151
02152 if (aadj < 2./FLT_RADIX)
02153 aadj = 1./FLT_RADIX;
02154 else
02155 aadj *= 0.5;
02156 aadj1 = -aadj;
02157 }
02158 }
02159 else {
02160 aadj *= 0.5;
02161 aadj1 = dsign ? aadj : -aadj;
02162 #ifdef Check_FLT_ROUNDS
02163 switch(Rounding) {
02164 case 2:
02165 aadj1 -= 0.5;
02166 break;
02167 case 0:
02168 case 3:
02169 aadj1 += 0.5;
02170 }
02171 #else
02172 if (Flt_Rounds == 0)
02173 aadj1 += 0.5;
02174 #endif
02175 }
02176 y = word0(rv) & Exp_mask;
02177
02178
02179
02180 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
02181 dval(rv0) = dval(rv);
02182 word0(rv) -= P*Exp_msk1;
02183 adj = aadj1 * ulp(dval(rv));
02184 dval(rv) += adj;
02185 if ((word0(rv) & Exp_mask) >=
02186 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
02187 if (word0(rv0) == Big0 && word1(rv0) == Big1)
02188 goto ovfl;
02189 word0(rv) = Big0;
02190 word1(rv) = Big1;
02191 goto cont;
02192 }
02193 else
02194 word0(rv) += P*Exp_msk1;
02195 }
02196 else {
02197 #ifdef Avoid_Underflow
02198 if (scale && y <= 2*P*Exp_msk1) {
02199 if (aadj <= 0x7fffffff) {
02200 if ((z = (ULong)aadj) <= 0)
02201 z = 1;
02202 aadj = z;
02203 aadj1 = dsign ? aadj : -aadj;
02204 }
02205 dval(aadj2) = aadj1;
02206 word0(aadj2) += (2*P+1)*Exp_msk1 - y;
02207 aadj1 = dval(aadj2);
02208 }
02209 adj = aadj1 * ulp(dval(rv));
02210 dval(rv) += adj;
02211 #else
02212 #ifdef Sudden_Underflow
02213 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02214 dval(rv0) = dval(rv);
02215 word0(rv) += P*Exp_msk1;
02216 adj = aadj1 * ulp(dval(rv));
02217 dval(rv) += adj;
02218 #ifdef IBM
02219 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
02220 #else
02221 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
02222 #endif
02223 {
02224 if (word0(rv0) == Tiny0
02225 && word1(rv0) == Tiny1)
02226 goto undfl;
02227 word0(rv) = Tiny0;
02228 word1(rv) = Tiny1;
02229 goto cont;
02230 }
02231 else
02232 word0(rv) -= P*Exp_msk1;
02233 }
02234 else {
02235 adj = aadj1 * ulp(dval(rv));
02236 dval(rv) += adj;
02237 }
02238 #else
02239
02240
02241
02242
02243
02244
02245
02246 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
02247 aadj1 = (double)(int)(aadj + 0.5);
02248 if (!dsign)
02249 aadj1 = -aadj1;
02250 }
02251 adj = aadj1 * ulp(dval(rv));
02252 dval(rv) += adj;
02253 #endif
02254 #endif
02255 }
02256 z = word0(rv) & Exp_mask;
02257 #ifndef SET_INEXACT
02258 #ifdef Avoid_Underflow
02259 if (!scale)
02260 #endif
02261 if (y == z) {
02262
02263 L = (Long)aadj;
02264 aadj -= L;
02265
02266 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
02267 if (aadj < .4999999 || aadj > .5000001)
02268 break;
02269 }
02270 else if (aadj < .4999999/FLT_RADIX)
02271 break;
02272 }
02273 #endif
02274 cont:
02275 Bfree(bb);
02276 Bfree(bd);
02277 Bfree(bs);
02278 Bfree(delta);
02279 }
02280 #ifdef SET_INEXACT
02281 if (inexact) {
02282 if (!oldinexact) {
02283 word0(rv0) = Exp_1 + (70 << Exp_shift);
02284 word1(rv0) = 0;
02285 dval(rv0) += 1.;
02286 }
02287 }
02288 else if (!oldinexact)
02289 clear_inexact();
02290 #endif
02291 #ifdef Avoid_Underflow
02292 if (scale) {
02293 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
02294 word1(rv0) = 0;
02295 dval(rv) *= dval(rv0);
02296 #ifndef NO_ERRNO
02297
02298 if (word0(rv) == 0 && word1(rv) == 0)
02299 errno = ERANGE;
02300 #endif
02301 }
02302 #endif
02303 #ifdef SET_INEXACT
02304 if (inexact && !(word0(rv) & Exp_mask)) {
02305
02306 dval(rv0) = 1e-300;
02307 dval(rv0) *= dval(rv0);
02308 }
02309 #endif
02310 retfree:
02311 Bfree(bb);
02312 Bfree(bd);
02313 Bfree(bs);
02314 Bfree(bd0);
02315 Bfree(delta);
02316 ret:
02317 if (se)
02318 *se = (char *)s;
02319 return sign ? -dval(rv) : dval(rv);
02320 }
02321
02322 static int
02323 quorem
02324 (Bigint *b, Bigint *S)
02325 {
02326 int n;
02327 ULong *bx, *bxe, q, *sx, *sxe;
02328 #ifdef ULLong
02329 ULLong borrow, carry, y, ys;
02330 #else
02331 ULong borrow, carry, y, ys;
02332 #ifdef Pack_32
02333 ULong si, z, zs;
02334 #endif
02335 #endif
02336
02337 n = S->wds;
02338 #ifdef DEBUG
02339 if (b->wds > n)
02340 Bug("oversize b in quorem");
02341 #endif
02342 if (b->wds < n)
02343 return 0;
02344 sx = S->x;
02345 sxe = sx + --n;
02346 bx = b->x;
02347 bxe = bx + n;
02348 q = *bxe / (*sxe + 1);
02349 #ifdef DEBUG
02350 if (q > 9)
02351 Bug("oversized quotient in quorem");
02352 #endif
02353 if (q) {
02354 borrow = 0;
02355 carry = 0;
02356 do {
02357 #ifdef ULLong
02358 ys = *sx++ * (ULLong)q + carry;
02359 carry = ys >> 32;
02360 y = *bx - (ys & FFFFFFFF) - borrow;
02361 borrow = y >> 32 & (ULong)1;
02362 *bx++ = (ULong)y & FFFFFFFF;
02363 #else
02364 #ifdef Pack_32
02365 si = *sx++;
02366 ys = (si & 0xffff) * q + carry;
02367 zs = (si >> 16) * q + (ys >> 16);
02368 carry = zs >> 16;
02369 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
02370 borrow = (y & 0x10000) >> 16;
02371 z = (*bx >> 16) - (zs & 0xffff) - borrow;
02372 borrow = (z & 0x10000) >> 16;
02373 Storeinc(bx, z, y);
02374 #else
02375 ys = *sx++ * q + carry;
02376 carry = ys >> 16;
02377 y = *bx - (ys & 0xffff) - borrow;
02378 borrow = (y & 0x10000) >> 16;
02379 *bx++ = y & 0xffff;
02380 #endif
02381 #endif
02382 }
02383 while(sx <= sxe);
02384 if (!*bxe) {
02385 bx = b->x;
02386 while(--bxe > bx && !*bxe)
02387 --n;
02388 b->wds = n;
02389 }
02390 }
02391 if (cmp(b, S) >= 0) {
02392 q++;
02393 borrow = 0;
02394 carry = 0;
02395 bx = b->x;
02396 sx = S->x;
02397 do {
02398 #ifdef ULLong
02399 ys = *sx++ + carry;
02400 carry = ys >> 32;
02401 y = *bx - (ys & FFFFFFFF) - borrow;
02402 borrow = y >> 32 & (ULong)1;
02403 *bx++ = (ULong)y & FFFFFFFF;
02404 #else
02405 #ifdef Pack_32
02406 si = *sx++;
02407 ys = (si & 0xffff) + carry;
02408 zs = (si >> 16) + (ys >> 16);
02409 carry = zs >> 16;
02410 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
02411 borrow = (y & 0x10000) >> 16;
02412 z = (*bx >> 16) - (zs & 0xffff) - borrow;
02413 borrow = (z & 0x10000) >> 16;
02414 Storeinc(bx, z, y);
02415 #else
02416 ys = *sx++ + carry;
02417 carry = ys >> 16;
02418 y = *bx - (ys & 0xffff) - borrow;
02419 borrow = (y & 0x10000) >> 16;
02420 *bx++ = y & 0xffff;
02421 #endif
02422 #endif
02423 }
02424 while(sx <= sxe);
02425 bx = b->x;
02426 bxe = bx + n;
02427 if (!*bxe) {
02428 while(--bxe > bx && !*bxe)
02429 --n;
02430 b->wds = n;
02431 }
02432 }
02433 return q;
02434 }
02435
02436 #ifndef MULTIPLE_THREADS
02437 static char *dtoa_result;
02438 #endif
02439
02440 static char *
02441 rv_alloc(int i)
02442 {
02443 int j, k, *r;
02444
02445 j = sizeof(ULong);
02446 for(k = 0;
02447 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i;
02448 j <<= 1)
02449 k++;
02450 r = (int*)Balloc(k);
02451 *r = k;
02452 return
02453 #ifndef MULTIPLE_THREADS
02454 dtoa_result =
02455 #endif
02456 (char *)(r+1);
02457 }
02458
02459 static char *
02460 nrv_alloc(CONST char *s, char **rve, int n)
02461 {
02462 char *rv, *t;
02463
02464 t = rv = rv_alloc(n);
02465 while((*t = *s++)) t++;
02466 if (rve)
02467 *rve = t;
02468 return rv;
02469 }
02470
02471
02472
02473
02474
02475
02476
02477 void
02478 freedtoa(char *s)
02479 {
02480 Bigint *b = (Bigint *)((int *)s - 1);
02481 b->maxwds = 1 << (b->k = *(int*)b);
02482 Bfree(b);
02483 #ifndef MULTIPLE_THREADS
02484 if (s == dtoa_result)
02485 dtoa_result = 0;
02486 #endif
02487 }
02488
02489
02490
02491
02492
02493
02494
02495
02496
02497
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507
02508
02509
02510
02511
02512
02513
02514
02515
02516
02517
02518
02519
02520
02521
02522
02523 char *
02524 dtoa
02525 (double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
02526 {
02527
02528
02529
02530
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560
02561 int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0,
02562 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
02563 spec_case, try_quick;
02564 Long L;
02565 #ifndef Sudden_Underflow
02566 int denorm;
02567 ULong x;
02568 #endif
02569 Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
02570 U d, d2, eps;
02571 double ds;
02572 char *s, *s0;
02573 #ifdef Honor_FLT_ROUNDS
02574 int rounding;
02575 #endif
02576 #ifdef SET_INEXACT
02577 int inexact, oldinexact;
02578 #endif
02579
02580 #ifndef MULTIPLE_THREADS
02581 if (dtoa_result) {
02582 freedtoa(dtoa_result);
02583 dtoa_result = 0;
02584 }
02585 #endif
02586
02587 dval(d) = dd;
02588 if (word0(d) & Sign_bit) {
02589
02590 *sign = 1;
02591 word0(d) &= ~Sign_bit;
02592 }
02593 else
02594 *sign = 0;
02595
02596 #if defined(IEEE_Arith) + defined(VAX)
02597 #ifdef IEEE_Arith
02598 if ((word0(d) & Exp_mask) == Exp_mask)
02599 #else
02600 if (word0(d) == 0x8000)
02601 #endif
02602 {
02603
02604 *decpt = 9999;
02605 #ifdef IEEE_Arith
02606 if (!word1(d) && !(word0(d) & 0xfffff))
02607 return nrv_alloc("Infinity", rve, 8);
02608 #endif
02609 return nrv_alloc("NaN", rve, 3);
02610 }
02611 #endif
02612 #ifdef IBM
02613 dval(d) += 0;
02614 #endif
02615 if (!dval(d)) {
02616 *decpt = 1;
02617 return nrv_alloc("0", rve, 1);
02618 }
02619
02620 #ifdef SET_INEXACT
02621 try_quick = oldinexact = get_inexact();
02622 inexact = 1;
02623 #endif
02624 #ifdef Honor_FLT_ROUNDS
02625 if ((rounding = Flt_Rounds) >= 2) {
02626 if (*sign)
02627 rounding = rounding == 2 ? 0 : 2;
02628 else
02629 if (rounding != 2)
02630 rounding = 0;
02631 }
02632 #endif
02633
02634 b = d2b(dval(d), &be, &bbits);
02635 #ifdef Sudden_Underflow
02636 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
02637 #else
02638 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
02639 #endif
02640 dval(d2) = dval(d);
02641 word0(d2) &= Frac_mask1;
02642 word0(d2) |= Exp_11;
02643 #ifdef IBM
02644 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
02645 dval(d2) /= 1 << j;
02646 #endif
02647
02648
02649
02650
02651
02652
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669
02670 i -= Bias;
02671 #ifdef IBM
02672 i <<= 2;
02673 i += j;
02674 #endif
02675 #ifndef Sudden_Underflow
02676 denorm = 0;
02677 }
02678 else {
02679
02680
02681 i = bbits + be + (Bias + (P-1) - 1);
02682 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
02683 : word1(d) << 32 - i;
02684 dval(d2) = x;
02685 word0(d2) -= 31*Exp_msk1;
02686 i -= (Bias + (P-1) - 1) + 1;
02687 denorm = 1;
02688 }
02689 #endif
02690 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
02691 k = (int)ds;
02692 if (ds < 0. && ds != k)
02693 k--;
02694 k_check = 1;
02695 if (k >= 0 && k <= Ten_pmax) {
02696 if (dval(d) < tens[k])
02697 k--;
02698 k_check = 0;
02699 }
02700 j = bbits - i - 1;
02701 if (j >= 0) {
02702 b2 = 0;
02703 s2 = j;
02704 }
02705 else {
02706 b2 = -j;
02707 s2 = 0;
02708 }
02709 if (k >= 0) {
02710 b5 = 0;
02711 s5 = k;
02712 s2 += k;
02713 }
02714 else {
02715 b2 -= k;
02716 b5 = -k;
02717 s5 = 0;
02718 }
02719 if (mode < 0 || mode > 9)
02720 mode = 0;
02721
02722 #ifndef SET_INEXACT
02723 #ifdef Check_FLT_ROUNDS
02724 try_quick = Rounding == 1;
02725 #else
02726 try_quick = 1;
02727 #endif
02728 #endif
02729
02730 if (mode > 5) {
02731 mode -= 4;
02732 try_quick = 0;
02733 }
02734 leftright = 1;
02735 switch(mode) {
02736 case 0:
02737 case 1:
02738 ilim = ilim1 = -1;
02739 i = 18;
02740 ndigits = 0;
02741 break;
02742 case 2:
02743 leftright = 0;
02744
02745 case 4:
02746 if (ndigits <= 0)
02747 ndigits = 1;
02748 ilim = ilim1 = i = ndigits;
02749 break;
02750 case 3:
02751 leftright = 0;
02752
02753 case 5:
02754 i = ndigits + k + 1;
02755 ilim = i;
02756 ilim1 = i - 1;
02757 if (i <= 0)
02758 i = 1;
02759 }
02760 s = s0 = rv_alloc(i);
02761
02762 #ifdef Honor_FLT_ROUNDS
02763 if (mode > 1 && rounding != 1)
02764 leftright = 0;
02765 #endif
02766
02767 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
02768
02769
02770
02771 i = 0;
02772 dval(d2) = dval(d);
02773 k0 = k;
02774 ilim0 = ilim;
02775 ieps = 2;
02776 if (k > 0) {
02777 ds = tens[k&0xf];
02778 j = k >> 4;
02779 if (j & Bletch) {
02780
02781 j &= Bletch - 1;
02782 dval(d) /= bigtens[n_bigtens-1];
02783 ieps++;
02784 }
02785 for(; j; j >>= 1, i++)
02786 if (j & 1) {
02787 ieps++;
02788 ds *= bigtens[i];
02789 }
02790 dval(d) /= ds;
02791 }
02792 else if ((j1 = -k)) {
02793 dval(d) *= tens[j1 & 0xf];
02794 for(j = j1 >> 4; j; j >>= 1, i++)
02795 if (j & 1) {
02796 ieps++;
02797 dval(d) *= bigtens[i];
02798 }
02799 }
02800 if (k_check && dval(d) < 1. && ilim > 0) {
02801 if (ilim1 <= 0)
02802 goto fast_failed;
02803 ilim = ilim1;
02804 k--;
02805 dval(d) *= 10.;
02806 ieps++;
02807 }
02808 dval(eps) = ieps*dval(d) + 7.;
02809 word0(eps) -= (P-1)*Exp_msk1;
02810 if (ilim == 0) {
02811 S = mhi = 0;
02812 dval(d) -= 5.;
02813 if (dval(d) > dval(eps))
02814 goto one_digit;
02815 if (dval(d) < -dval(eps))
02816 goto no_digits;
02817 goto fast_failed;
02818 }
02819 #ifndef No_leftright
02820 if (leftright) {
02821
02822
02823
02824 dval(eps) = 0.5/tens[ilim-1] - dval(eps);
02825 for(i = 0;;) {
02826 L = (long int)dval(d);
02827 dval(d) -= L;
02828 *s++ = '0' + (int)L;
02829 if (dval(d) < dval(eps))
02830 goto ret1;
02831 if (1. - dval(d) < dval(eps))
02832 goto bump_up;
02833 if (++i >= ilim)
02834 break;
02835 dval(eps) *= 10.;
02836 dval(d) *= 10.;
02837 }
02838 }
02839 else {
02840 #endif
02841
02842 dval(eps) *= tens[ilim-1];
02843 for(i = 1;; i++, dval(d) *= 10.) {
02844 L = (Long)(dval(d));
02845 if (!(dval(d) -= L))
02846 ilim = i;
02847 *s++ = '0' + (int)L;
02848 if (i == ilim) {
02849 if (dval(d) > 0.5 + dval(eps))
02850 goto bump_up;
02851 else if (dval(d) < 0.5 - dval(eps)) {
02852 while(*--s == '0')
02853 ;
02854 s++;
02855 goto ret1;
02856 }
02857 break;
02858 }
02859 }
02860 #ifndef No_leftright
02861 }
02862 #endif
02863 fast_failed:
02864 s = s0;
02865 dval(d) = dval(d2);
02866 k = k0;
02867 ilim = ilim0;
02868 }
02869
02870
02871
02872 if (be >= 0 && k <= Int_max) {
02873
02874 ds = tens[k];
02875 if (ndigits < 0 && ilim <= 0) {
02876 S = mhi = 0;
02877 if (ilim < 0 || dval(d) <= 5*ds)
02878 goto no_digits;
02879 goto one_digit;
02880 }
02881 for(i = 1;; i++, dval(d) *= 10.) {
02882 L = (Long)(dval(d) / ds);
02883 dval(d) -= L*ds;
02884 #ifdef Check_FLT_ROUNDS
02885
02886 if (dval(d) < 0) {
02887 L--;
02888 dval(d) += ds;
02889 }
02890 #endif
02891 *s++ = '0' + (int)L;
02892 if (!dval(d)) {
02893 #ifdef SET_INEXACT
02894 inexact = 0;
02895 #endif
02896 break;
02897 }
02898 if (i == ilim) {
02899 #ifdef Honor_FLT_ROUNDS
02900 if (mode > 1)
02901 switch(rounding) {
02902 case 0: goto ret1;
02903 case 2: goto bump_up;
02904 }
02905 #endif
02906 dval(d) += dval(d);
02907 if (dval(d) > ds || dval(d) == ds && L & 1) {
02908 bump_up:
02909 while(*--s == '9')
02910 if (s == s0) {
02911 k++;
02912 *s = '0';
02913 break;
02914 }
02915 ++*s++;
02916 }
02917 break;
02918 }
02919 }
02920 goto ret1;
02921 }
02922
02923 m2 = b2;
02924 m5 = b5;
02925 mhi = mlo = 0;
02926 if (leftright) {
02927 i =
02928 #ifndef Sudden_Underflow
02929 denorm ? be + (Bias + (P-1) - 1 + 1) :
02930 #endif
02931 #ifdef IBM
02932 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
02933 #else
02934 1 + P - bbits;
02935 #endif
02936 b2 += i;
02937 s2 += i;
02938 mhi = i2b(1);
02939 }
02940 if (m2 > 0 && s2 > 0) {
02941 i = m2 < s2 ? m2 : s2;
02942 b2 -= i;
02943 m2 -= i;
02944 s2 -= i;
02945 }
02946 if (b5 > 0) {
02947 if (leftright) {
02948 if (m5 > 0) {
02949 mhi = pow5mult(mhi, m5);
02950 b1 = mult(mhi, b);
02951 Bfree(b);
02952 b = b1;
02953 }
02954 if ((j = b5 - m5))
02955 b = pow5mult(b, j);
02956 }
02957 else
02958 b = pow5mult(b, b5);
02959 }
02960 S = i2b(1);
02961 if (s5 > 0)
02962 S = pow5mult(S, s5);
02963
02964
02965
02966 spec_case = 0;
02967 if ((mode < 2 || leftright)
02968 #ifdef Honor_FLT_ROUNDS
02969 && rounding == 1
02970 #endif
02971 ) {
02972 if (!word1(d) && !(word0(d) & Bndry_mask)
02973 #ifndef Sudden_Underflow
02974 && word0(d) & (Exp_mask & ~Exp_msk1)
02975 #endif
02976 ) {
02977
02978 b2 += Log2P;
02979 s2 += Log2P;
02980 spec_case = 1;
02981 }
02982 }
02983
02984
02985
02986
02987
02988
02989
02990
02991 #ifdef Pack_32
02992 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
02993 i = 32 - i;
02994 #else
02995 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
02996 i = 16 - i;
02997 #endif
02998 if (i > 4) {
02999 i -= 4;
03000 b2 += i;
03001 m2 += i;
03002 s2 += i;
03003 }
03004 else if (i < 4) {
03005 i += 28;
03006 b2 += i;
03007 m2 += i;
03008 s2 += i;
03009 }
03010 if (b2 > 0)
03011 b = lshift(b, b2);
03012 if (s2 > 0)
03013 S = lshift(S, s2);
03014 if (k_check) {
03015 if (cmp(b,S) < 0) {
03016 k--;
03017 b = multadd(b, 10, 0);
03018 if (leftright)
03019 mhi = multadd(mhi, 10, 0);
03020 ilim = ilim1;
03021 }
03022 }
03023 if (ilim <= 0 && (mode == 3 || mode == 5)) {
03024 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
03025
03026 no_digits:
03027 k = -1 - ndigits;
03028 goto ret;
03029 }
03030 one_digit:
03031 *s++ = '1';
03032 k++;
03033 goto ret;
03034 }
03035 if (leftright) {
03036 if (m2 > 0)
03037 mhi = lshift(mhi, m2);
03038
03039
03040
03041
03042
03043 mlo = mhi;
03044 if (spec_case) {
03045 mhi = Balloc(mhi->k);
03046 Bcopy(mhi, mlo);
03047 mhi = lshift(mhi, Log2P);
03048 }
03049
03050 for(i = 1;;i++) {
03051 dig = quorem(b,S) + '0';
03052
03053
03054
03055 j = cmp(b, mlo);
03056 delta = diff(S, mhi);
03057 j1 = delta->sign ? 1 : cmp(b, delta);
03058 Bfree(delta);
03059 #ifndef ROUND_BIASED
03060 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
03061 #ifdef Honor_FLT_ROUNDS
03062 && rounding >= 1
03063 #endif
03064 ) {
03065 if (dig == '9')
03066 goto round_9_up;
03067 if (j > 0)
03068 dig++;
03069 #ifdef SET_INEXACT
03070 else if (!b->x[0] && b->wds <= 1)
03071 inexact = 0;
03072 #endif
03073 *s++ = dig;
03074 goto ret;
03075 }
03076 #endif
03077 if (j < 0 || j == 0 && mode != 1
03078 #ifndef ROUND_BIASED
03079 && !(word1(d) & 1)
03080 #endif
03081 ) {
03082 if (!b->x[0] && b->wds <= 1) {
03083 #ifdef SET_INEXACT
03084 inexact = 0;
03085 #endif
03086 goto accept_dig;
03087 }
03088 #ifdef Honor_FLT_ROUNDS
03089 if (mode > 1)
03090 switch(rounding) {
03091 case 0: goto accept_dig;
03092 case 2: goto keep_dig;
03093 }
03094 #endif
03095 if (j1 > 0) {
03096 b = lshift(b, 1);
03097 j1 = cmp(b, S);
03098 if ((j1 > 0 || j1 == 0 && dig & 1)
03099 && dig++ == '9')
03100 goto round_9_up;
03101 }
03102 accept_dig:
03103 *s++ = dig;
03104 goto ret;
03105 }
03106 if (j1 > 0) {
03107 #ifdef Honor_FLT_ROUNDS
03108 if (!rounding)
03109 goto accept_dig;
03110 #endif
03111 if (dig == '9') {
03112 round_9_up:
03113 *s++ = '9';
03114 goto roundoff;
03115 }
03116 *s++ = dig + 1;
03117 goto ret;
03118 }
03119 #ifdef Honor_FLT_ROUNDS
03120 keep_dig:
03121 #endif
03122 *s++ = dig;
03123 if (i == ilim)
03124 break;
03125 b = multadd(b, 10, 0);
03126 if (mlo == mhi)
03127 mlo = mhi = multadd(mhi, 10, 0);
03128 else {
03129 mlo = multadd(mlo, 10, 0);
03130 mhi = multadd(mhi, 10, 0);
03131 }
03132 }
03133 }
03134 else
03135 for(i = 1;; i++) {
03136 *s++ = dig = quorem(b,S) + '0';
03137 if (!b->x[0] && b->wds <= 1) {
03138 #ifdef SET_INEXACT
03139 inexact = 0;
03140 #endif
03141 goto ret;
03142 }
03143 if (i >= ilim)
03144 break;
03145 b = multadd(b, 10, 0);
03146 }
03147
03148
03149
03150 #ifdef Honor_FLT_ROUNDS
03151 switch(rounding) {
03152 case 0: goto trimzeros;
03153 case 2: goto roundoff;
03154 }
03155 #endif
03156 b = lshift(b, 1);
03157 j = cmp(b, S);
03158 if (j > 0 || j == 0 && dig & 1) {
03159 roundoff:
03160 while(*--s == '9')
03161 if (s == s0) {
03162 k++;
03163 *s++ = '1';
03164 goto ret;
03165 }
03166 ++*s++;
03167 }
03168 else {
03169 #ifdef Honor_FLT_ROUNDS
03170 trimzeros:
03171 #endif
03172 while(*--s == '0')
03173 ;
03174 s++;
03175 }
03176 ret:
03177 Bfree(S);
03178 if (mhi) {
03179 if (mlo && mlo != mhi)
03180 Bfree(mlo);
03181 Bfree(mhi);
03182 }
03183 ret1:
03184 #ifdef SET_INEXACT
03185 if (inexact) {
03186 if (!oldinexact) {
03187 word0(d) = Exp_1 + (70 << Exp_shift);
03188 word1(d) = 0;
03189 dval(d) += 1.;
03190 }
03191 }
03192 else if (!oldinexact)
03193 clear_inexact();
03194 #endif
03195 Bfree(b);
03196 *s = 0;
03197 *decpt = k + 1;
03198 if (rve)
03199 *rve = s;
03200 return s0;
03201 }
03202 #ifdef __cplusplus
03203 }
03204 #endif