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