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