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