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