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