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 #include "config.h"
00074 #ifdef WORDS_BIGENDIAN
00075 #define IEEE_MC68k
00076 #else
00077 #define IEEE_8087
00078 #endif
00079 #ifndef HAVE_LONG_LONG
00080 #define NO_LONG_LONG
00081 #endif
00082 #define Omit_Private_Memory
00083 #include "ckd_alloc.h"
00084 #undef USE_LOCALE
00085
00086
00087 #include "prim_type.h"
00088 #define Long int32
00089 #define ULong uint32
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 #ifndef Long
00196 #define Long long
00197 #endif
00198 #ifndef ULong
00199 typedef unsigned Long ULong;
00200 #endif
00201
00202 #ifdef DEBUG
00203 #include "stdio.h"
00204 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
00205 #endif
00206
00207 #include "stdlib.h"
00208 #include "string.h"
00209
00210 #ifdef USE_LOCALE
00211 #include "locale.h"
00212 #endif
00213
00214
00215
00216 #undef IEEE_Arith
00217 #undef Avoid_Underflow
00218 #ifdef IEEE_MC68k
00219 #define IEEE_Arith
00220 #endif
00221 #ifdef IEEE_8087
00222 #define IEEE_Arith
00223 #endif
00224
00225 #ifdef IEEE_Arith
00226 #ifndef NO_INFNAN_CHECK
00227 #undef INFNAN_CHECK
00228 #define INFNAN_CHECK
00229 #endif
00230 #else
00231 #undef INFNAN_CHECK
00232 #endif
00233
00234 #include "errno.h"
00235
00236 #ifdef Bad_float_h
00237
00238 #ifdef IEEE_Arith
00239 #define DBL_DIG 15
00240 #define DBL_MAX_10_EXP 308
00241 #define DBL_MAX_EXP 1024
00242 #define FLT_RADIX 2
00243 #endif
00244
00245 #ifdef IBM
00246 #define DBL_DIG 16
00247 #define DBL_MAX_10_EXP 75
00248 #define DBL_MAX_EXP 63
00249 #define FLT_RADIX 16
00250 #define DBL_MAX 7.2370055773322621e+75
00251 #endif
00252
00253 #ifdef VAX
00254 #define DBL_DIG 16
00255 #define DBL_MAX_10_EXP 38
00256 #define DBL_MAX_EXP 127
00257 #define FLT_RADIX 2
00258 #define DBL_MAX 1.7014118346046923e+38
00259 #endif
00260
00261 #ifndef LONG_MAX
00262 #define LONG_MAX 2147483647
00263 #endif
00264
00265 #else
00266 #include "float.h"
00267 #endif
00268
00269 #ifndef __MATH_H__
00270 #include "math.h"
00271 #endif
00272
00273 #ifdef __cplusplus
00274 extern "C" {
00275 #endif
00276
00277 #ifndef CONST
00278 #ifdef KR_headers
00279 #define CONST
00280 #else
00281 #define CONST const
00282 #endif
00283 #endif
00284
00285
00286 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
00287 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
00288 #endif
00289
00291 typedef union { double d; ULong L[2]; } U;
00292
00293 #ifdef YES_ALIAS
00294 #define dval(x) x
00295 #ifdef IEEE_8087
00296 #define word0(x) ((ULong *)&x)[1]
00297 #define word1(x) ((ULong *)&x)[0]
00298 #else
00299 #define word0(x) ((ULong *)&x)[0]
00300 #define word1(x) ((ULong *)&x)[1]
00301 #endif
00302 #else
00303 #ifdef IEEE_8087
00304 #define word0(x) ((U*)&x)->L[1]
00305 #define word1(x) ((U*)&x)->L[0]
00306 #else
00307 #define word0(x) ((U*)&x)->L[0]
00308 #define word1(x) ((U*)&x)->L[1]
00309 #endif
00310 #define dval(x) ((U*)&x)->d
00311 #endif
00312
00313
00314
00315
00316
00317 #if defined(IEEE_8087) + defined(VAX)
00318 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
00319 ((unsigned short *)a)[0] = (unsigned short)c, a++)
00320 #else
00321 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
00322 ((unsigned short *)a)[1] = (unsigned short)c, a++)
00323 #endif
00324
00325
00326
00327
00328
00329
00330
00331 #ifdef IEEE_Arith
00332 #define Exp_shift 20
00333 #define Exp_shift1 20
00334 #define Exp_msk1 0x100000
00335 #define Exp_msk11 0x100000
00336 #define Exp_mask 0x7ff00000
00337 #define P 53
00338 #define Bias 1023
00339 #define Emin (-1022)
00340 #define Exp_1 0x3ff00000
00341 #define Exp_11 0x3ff00000
00342 #define Ebits 11
00343 #define Frac_mask 0xfffff
00344 #define Frac_mask1 0xfffff
00345 #define Ten_pmax 22
00346 #define Bletch 0x10
00347 #define Bndry_mask 0xfffff
00348 #define Bndry_mask1 0xfffff
00349 #define LSB 1
00350 #define Sign_bit 0x80000000
00351 #define Log2P 1
00352 #define Tiny0 0
00353 #define Tiny1 1
00354 #define Quick_max 14
00355 #define Int_max 14
00356 #ifndef NO_IEEE_Scale
00357 #define Avoid_Underflow
00358 #ifdef Flush_Denorm
00359 #undef Sudden_Underflow
00360 #endif
00361 #endif
00362
00363 #ifndef Flt_Rounds
00364 #ifdef FLT_ROUNDS
00365 #define Flt_Rounds FLT_ROUNDS
00366 #else
00367 #define Flt_Rounds 1
00368 #endif
00369 #endif
00370
00371 #ifdef Honor_FLT_ROUNDS
00372 #define Rounding rounding
00373 #undef Check_FLT_ROUNDS
00374 #define Check_FLT_ROUNDS
00375 #else
00376 #define Rounding Flt_Rounds
00377 #endif
00378
00379 #else
00380 #undef Check_FLT_ROUNDS
00381 #undef Honor_FLT_ROUNDS
00382 #undef SET_INEXACT
00383 #undef Sudden_Underflow
00384 #define Sudden_Underflow
00385 #ifdef IBM
00386 #undef Flt_Rounds
00387 #define Flt_Rounds 0
00388 #define Exp_shift 24
00389 #define Exp_shift1 24
00390 #define Exp_msk1 0x1000000
00391 #define Exp_msk11 0x1000000
00392 #define Exp_mask 0x7f000000
00393 #define P 14
00394 #define Bias 65
00395 #define Exp_1 0x41000000
00396 #define Exp_11 0x41000000
00397 #define Ebits 8
00398 #define Frac_mask 0xffffff
00399 #define Frac_mask1 0xffffff
00400 #define Bletch 4
00401 #define Ten_pmax 22
00402 #define Bndry_mask 0xefffff
00403 #define Bndry_mask1 0xffffff
00404 #define LSB 1
00405 #define Sign_bit 0x80000000
00406 #define Log2P 4
00407 #define Tiny0 0x100000
00408 #define Tiny1 0
00409 #define Quick_max 14
00410 #define Int_max 15
00411 #else
00412 #undef Flt_Rounds
00413 #define Flt_Rounds 1
00414 #define Exp_shift 23
00415 #define Exp_shift1 7
00416 #define Exp_msk1 0x80
00417 #define Exp_msk11 0x800000
00418 #define Exp_mask 0x7f80
00419 #define P 56
00420 #define Bias 129
00421 #define Exp_1 0x40800000
00422 #define Exp_11 0x4080
00423 #define Ebits 8
00424 #define Frac_mask 0x7fffff
00425 #define Frac_mask1 0xffff007f
00426 #define Ten_pmax 24
00427 #define Bletch 2
00428 #define Bndry_mask 0xffff007f
00429 #define Bndry_mask1 0xffff007f
00430 #define LSB 0x10000
00431 #define Sign_bit 0x8000
00432 #define Log2P 1
00433 #define Tiny0 0x80
00434 #define Tiny1 0
00435 #define Quick_max 15
00436 #define Int_max 15
00437 #endif
00438 #endif
00439
00440 #ifndef IEEE_Arith
00441 #define ROUND_BIASED
00442 #endif
00443
00444 #ifdef RND_PRODQUOT
00445 #define rounded_product(a,b) a = rnd_prod(a, b)
00446 #define rounded_quotient(a,b) a = rnd_quot(a, b)
00447 #ifdef KR_headers
00448 extern double rnd_prod(), rnd_quot();
00449 #else
00450 extern double rnd_prod(double, double), rnd_quot(double, double);
00451 #endif
00452 #else
00453 #define rounded_product(a,b) a *= b
00454 #define rounded_quotient(a,b) a /= b
00455 #endif
00456
00457 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
00458 #define Big1 0xffffffff
00459
00460 #ifndef Pack_32
00461 #define Pack_32
00462 #endif
00463
00464 #ifdef KR_headers
00465 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
00466 #else
00467 #define FFFFFFFF 0xffffffffUL
00468 #endif
00469
00470 #ifdef NO_LONG_LONG
00471 #undef ULLong
00472 #ifdef Just_16
00473 #undef Pack_32
00474
00475
00476
00477
00478
00479 #endif
00480 #else
00481 #ifndef Llong
00482 #define Llong long long
00483 #endif
00484 #ifndef ULLong
00485 #define ULLong unsigned Llong
00486 #endif
00487 #endif
00488
00489 #ifndef MULTIPLE_THREADS
00490 #define ACQUIRE_DTOA_LOCK(n)
00491 #define FREE_DTOA_LOCK(n)
00492 #endif
00493
00494 #define Kmax 15
00495
00496 #ifdef __cplusplus
00497 extern "C" double sb_strtod(const char *s00, char **se);
00498 #endif
00499
00500 struct
00501 Bigint {
00502 struct Bigint *next;
00503 int k, maxwds, sign, wds;
00504 ULong x[1];
00505 };
00506
00507 typedef struct Bigint Bigint;
00508
00509 static Bigint *
00510 Balloc
00511 #ifdef KR_headers
00512 (k) int k;
00513 #else
00514 (int k)
00515 #endif
00516 {
00517 int x;
00518 size_t len;
00519 Bigint *rv;
00520
00521 x = 1 << k;
00522 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
00523 /sizeof(double);
00524 rv = ckd_malloc(len*sizeof(double));
00525 rv->k = k;
00526 rv->maxwds = x;
00527 rv->sign = rv->wds = 0;
00528 return rv;
00529 }
00530
00531 static void
00532 Bfree
00533 #ifdef KR_headers
00534 (v) Bigint *v;
00535 #else
00536 (Bigint *v)
00537 #endif
00538 {
00539 ckd_free(v);
00540 }
00541
00542 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
00543 y->wds*sizeof(Long) + 2*sizeof(int))
00544
00545 static Bigint *
00546 multadd
00547 #ifdef KR_headers
00548 (b, m, a) Bigint *b; int m, a;
00549 #else
00550 (Bigint *b, int m, int a)
00551 #endif
00552 {
00553 int i, wds;
00554 #ifdef ULLong
00555 ULong *x;
00556 ULLong carry, y;
00557 #else
00558 ULong carry, *x, y;
00559 #ifdef Pack_32
00560 ULong xi, z;
00561 #endif
00562 #endif
00563 Bigint *b1;
00564
00565 wds = b->wds;
00566 x = b->x;
00567 i = 0;
00568 carry = a;
00569 do {
00570 #ifdef ULLong
00571 y = *x * (ULLong)m + carry;
00572 carry = y >> 32;
00573 *x++ = y & FFFFFFFF;
00574 #else
00575 #ifdef Pack_32
00576 xi = *x;
00577 y = (xi & 0xffff) * m + carry;
00578 z = (xi >> 16) * m + (y >> 16);
00579 carry = z >> 16;
00580 *x++ = (z << 16) + (y & 0xffff);
00581 #else
00582 y = *x * m + carry;
00583 carry = y >> 16;
00584 *x++ = y & 0xffff;
00585 #endif
00586 #endif
00587 }
00588 while(++i < wds);
00589 if (carry) {
00590 if (wds >= b->maxwds) {
00591 b1 = Balloc(b->k+1);
00592 Bcopy(b1, b);
00593 Bfree(b);
00594 b = b1;
00595 }
00596 b->x[wds++] = carry;
00597 b->wds = wds;
00598 }
00599 return b;
00600 }
00601
00602 static Bigint *
00603 s2b
00604 #ifdef KR_headers
00605 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
00606 #else
00607 (CONST char *s, int nd0, int nd, ULong y9)
00608 #endif
00609 {
00610 Bigint *b;
00611 int i, k;
00612 Long x, y;
00613
00614 x = (nd + 8) / 9;
00615 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
00616 #ifdef Pack_32
00617 b = Balloc(k);
00618 b->x[0] = y9;
00619 b->wds = 1;
00620 #else
00621 b = Balloc(k+1);
00622 b->x[0] = y9 & 0xffff;
00623 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
00624 #endif
00625
00626 i = 9;
00627 if (9 < nd0) {
00628 s += 9;
00629 do b = multadd(b, 10, *s++ - '0');
00630 while(++i < nd0);
00631 s++;
00632 }
00633 else
00634 s += 10;
00635 for(; i < nd; i++)
00636 b = multadd(b, 10, *s++ - '0');
00637 return b;
00638 }
00639
00640 static int
00641 hi0bits
00642 #ifdef KR_headers
00643 (x) register ULong x;
00644 #else
00645 (register ULong x)
00646 #endif
00647 {
00648 register int k = 0;
00649
00650 if (!(x & 0xffff0000)) {
00651 k = 16;
00652 x <<= 16;
00653 }
00654 if (!(x & 0xff000000)) {
00655 k += 8;
00656 x <<= 8;
00657 }
00658 if (!(x & 0xf0000000)) {
00659 k += 4;
00660 x <<= 4;
00661 }
00662 if (!(x & 0xc0000000)) {
00663 k += 2;
00664 x <<= 2;
00665 }
00666 if (!(x & 0x80000000)) {
00667 k++;
00668 if (!(x & 0x40000000))
00669 return 32;
00670 }
00671 return k;
00672 }
00673
00674 static int
00675 lo0bits
00676 #ifdef KR_headers
00677 (y) ULong *y;
00678 #else
00679 (ULong *y)
00680 #endif
00681 {
00682 register int k;
00683 register ULong x = *y;
00684
00685 if (x & 7) {
00686 if (x & 1)
00687 return 0;
00688 if (x & 2) {
00689 *y = x >> 1;
00690 return 1;
00691 }
00692 *y = x >> 2;
00693 return 2;
00694 }
00695 k = 0;
00696 if (!(x & 0xffff)) {
00697 k = 16;
00698 x >>= 16;
00699 }
00700 if (!(x & 0xff)) {
00701 k += 8;
00702 x >>= 8;
00703 }
00704 if (!(x & 0xf)) {
00705 k += 4;
00706 x >>= 4;
00707 }
00708 if (!(x & 0x3)) {
00709 k += 2;
00710 x >>= 2;
00711 }
00712 if (!(x & 1)) {
00713 k++;
00714 x >>= 1;
00715 if (!x)
00716 return 32;
00717 }
00718 *y = x;
00719 return k;
00720 }
00721
00722 static Bigint *
00723 i2b
00724 #ifdef KR_headers
00725 (i) int i;
00726 #else
00727 (int i)
00728 #endif
00729 {
00730 Bigint *b;
00731
00732 b = Balloc(1);
00733 b->x[0] = i;
00734 b->wds = 1;
00735 return b;
00736 }
00737
00738 static Bigint *
00739 mult
00740 #ifdef KR_headers
00741 (a, b) Bigint *a, *b;
00742 #else
00743 (Bigint *a, Bigint *b)
00744 #endif
00745 {
00746 Bigint *c;
00747 int k, wa, wb, wc;
00748 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
00749 ULong y;
00750 #ifdef ULLong
00751 ULLong carry, z;
00752 #else
00753 ULong carry, z;
00754 #ifdef Pack_32
00755 ULong z2;
00756 #endif
00757 #endif
00758
00759 if (a->wds < b->wds) {
00760 c = a;
00761 a = b;
00762 b = c;
00763 }
00764 k = a->k;
00765 wa = a->wds;
00766 wb = b->wds;
00767 wc = wa + wb;
00768 if (wc > a->maxwds)
00769 k++;
00770 c = Balloc(k);
00771 for(x = c->x, xa = x + wc; x < xa; x++)
00772 *x = 0;
00773 xa = a->x;
00774 xae = xa + wa;
00775 xb = b->x;
00776 xbe = xb + wb;
00777 xc0 = c->x;
00778 #ifdef ULLong
00779 for(; xb < xbe; xc0++) {
00780 if ((y = *xb++)) {
00781 x = xa;
00782 xc = xc0;
00783 carry = 0;
00784 do {
00785 z = *x++ * (ULLong)y + *xc + carry;
00786 carry = z >> 32;
00787 *xc++ = z & FFFFFFFF;
00788 }
00789 while(x < xae);
00790 *xc = carry;
00791 }
00792 }
00793 #else
00794 #ifdef Pack_32
00795 for(; xb < xbe; xb++, xc0++) {
00796 if (y = *xb & 0xffff) {
00797 x = xa;
00798 xc = xc0;
00799 carry = 0;
00800 do {
00801 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
00802 carry = z >> 16;
00803 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
00804 carry = z2 >> 16;
00805 Storeinc(xc, z2, z);
00806 }
00807 while(x < xae);
00808 *xc = carry;
00809 }
00810 if (y = *xb >> 16) {
00811 x = xa;
00812 xc = xc0;
00813 carry = 0;
00814 z2 = *xc;
00815 do {
00816 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
00817 carry = z >> 16;
00818 Storeinc(xc, z, z2);
00819 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
00820 carry = z2 >> 16;
00821 }
00822 while(x < xae);
00823 *xc = z2;
00824 }
00825 }
00826 #else
00827 for(; xb < xbe; xc0++) {
00828 if (y = *xb++) {
00829 x = xa;
00830 xc = xc0;
00831 carry = 0;
00832 do {
00833 z = *x++ * y + *xc + carry;
00834 carry = z >> 16;
00835 *xc++ = z & 0xffff;
00836 }
00837 while(x < xae);
00838 *xc = carry;
00839 }
00840 }
00841 #endif
00842 #endif
00843 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
00844 c->wds = wc;
00845 return c;
00846 }
00847
00848 static Bigint *
00849 pow5mult
00850 #ifdef KR_headers
00851 (b, k) Bigint *b; int k;
00852 #else
00853 (Bigint *b, int k)
00854 #endif
00855 {
00856 Bigint *b1, *p5, *p51;
00857 int i;
00858 static int CONST p05[3] = { 5, 25, 125 };
00859
00860 if ((i = k & 3))
00861 b = multadd(b, p05[i-1], 0);
00862
00863 if (!(k >>= 2))
00864 return b;
00865
00866 p5 = i2b(625);
00867 for(;;) {
00868 if (k & 1) {
00869 b1 = mult(b, p5);
00870 Bfree(b);
00871 b = b1;
00872 }
00873 if (!(k >>= 1))
00874 break;
00875 p51 = mult(p5,p5);
00876 Bfree(p5);
00877 p5 = p51;
00878 }
00879 Bfree(p5);
00880 return b;
00881 }
00882
00883 static Bigint *
00884 lshift
00885 #ifdef KR_headers
00886 (b, k) Bigint *b; int k;
00887 #else
00888 (Bigint *b, int k)
00889 #endif
00890 {
00891 int i, k1, n, n1;
00892 Bigint *b1;
00893 ULong *x, *x1, *xe, z;
00894
00895 #ifdef Pack_32
00896 n = k >> 5;
00897 #else
00898 n = k >> 4;
00899 #endif
00900 k1 = b->k;
00901 n1 = n + b->wds + 1;
00902 for(i = b->maxwds; n1 > i; i <<= 1)
00903 k1++;
00904 b1 = Balloc(k1);
00905 x1 = b1->x;
00906 for(i = 0; i < n; i++)
00907 *x1++ = 0;
00908 x = b->x;
00909 xe = x + b->wds;
00910 #ifdef Pack_32
00911 if (k &= 0x1f) {
00912 k1 = 32 - k;
00913 z = 0;
00914 do {
00915 *x1++ = *x << k | z;
00916 z = *x++ >> k1;
00917 }
00918 while(x < xe);
00919 if ((*x1 = z))
00920 ++n1;
00921 }
00922 #else
00923 if (k &= 0xf) {
00924 k1 = 16 - k;
00925 z = 0;
00926 do {
00927 *x1++ = *x << k & 0xffff | z;
00928 z = *x++ >> k1;
00929 }
00930 while(x < xe);
00931 if (*x1 = z)
00932 ++n1;
00933 }
00934 #endif
00935 else do
00936 *x1++ = *x++;
00937 while(x < xe);
00938 b1->wds = n1 - 1;
00939 Bfree(b);
00940 return b1;
00941 }
00942
00943 static int
00944 cmp
00945 #ifdef KR_headers
00946 (a, b) Bigint *a, *b;
00947 #else
00948 (Bigint *a, Bigint *b)
00949 #endif
00950 {
00951 ULong *xa, *xa0, *xb, *xb0;
00952 int i, j;
00953
00954 i = a->wds;
00955 j = b->wds;
00956 #ifdef DEBUG
00957 if (i > 1 && !a->x[i-1])
00958 Bug("cmp called with a->x[a->wds-1] == 0");
00959 if (j > 1 && !b->x[j-1])
00960 Bug("cmp called with b->x[b->wds-1] == 0");
00961 #endif
00962 if (i -= j)
00963 return i;
00964 xa0 = a->x;
00965 xa = xa0 + j;
00966 xb0 = b->x;
00967 xb = xb0 + j;
00968 for(;;) {
00969 if (*--xa != *--xb)
00970 return *xa < *xb ? -1 : 1;
00971 if (xa <= xa0)
00972 break;
00973 }
00974 return 0;
00975 }
00976
00977 static Bigint *
00978 diff
00979 #ifdef KR_headers
00980 (a, b) Bigint *a, *b;
00981 #else
00982 (Bigint *a, Bigint *b)
00983 #endif
00984 {
00985 Bigint *c;
00986 int i, wa, wb;
00987 ULong *xa, *xae, *xb, *xbe, *xc;
00988 #ifdef ULLong
00989 ULLong borrow, y;
00990 #else
00991 ULong borrow, y;
00992 #ifdef Pack_32
00993 ULong z;
00994 #endif
00995 #endif
00996
00997 i = cmp(a,b);
00998 if (!i) {
00999 c = Balloc(0);
01000 c->wds = 1;
01001 c->x[0] = 0;
01002 return c;
01003 }
01004 if (i < 0) {
01005 c = a;
01006 a = b;
01007 b = c;
01008 i = 1;
01009 }
01010 else
01011 i = 0;
01012 c = Balloc(a->k);
01013 c->sign = i;
01014 wa = a->wds;
01015 xa = a->x;
01016 xae = xa + wa;
01017 wb = b->wds;
01018 xb = b->x;
01019 xbe = xb + wb;
01020 xc = c->x;
01021 borrow = 0;
01022 #ifdef ULLong
01023 do {
01024 y = (ULLong)*xa++ - *xb++ - borrow;
01025 borrow = y >> 32 & (ULong)1;
01026 *xc++ = y & FFFFFFFF;
01027 }
01028 while(xb < xbe);
01029 while(xa < xae) {
01030 y = *xa++ - borrow;
01031 borrow = y >> 32 & (ULong)1;
01032 *xc++ = y & FFFFFFFF;
01033 }
01034 #else
01035 #ifdef Pack_32
01036 do {
01037 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
01038 borrow = (y & 0x10000) >> 16;
01039 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
01040 borrow = (z & 0x10000) >> 16;
01041 Storeinc(xc, z, y);
01042 }
01043 while(xb < xbe);
01044 while(xa < xae) {
01045 y = (*xa & 0xffff) - borrow;
01046 borrow = (y & 0x10000) >> 16;
01047 z = (*xa++ >> 16) - borrow;
01048 borrow = (z & 0x10000) >> 16;
01049 Storeinc(xc, z, y);
01050 }
01051 #else
01052 do {
01053 y = *xa++ - *xb++ - borrow;
01054 borrow = (y & 0x10000) >> 16;
01055 *xc++ = y & 0xffff;
01056 }
01057 while(xb < xbe);
01058 while(xa < xae) {
01059 y = *xa++ - borrow;
01060 borrow = (y & 0x10000) >> 16;
01061 *xc++ = y & 0xffff;
01062 }
01063 #endif
01064 #endif
01065 while(!*--xc)
01066 wa--;
01067 c->wds = wa;
01068 return c;
01069 }
01070
01071 static double
01072 ulp
01073 #ifdef KR_headers
01074 (x) double x;
01075 #else
01076 (double x)
01077 #endif
01078 {
01079 register Long L;
01080 double a;
01081
01082 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
01083 #ifndef Avoid_Underflow
01084 #ifndef Sudden_Underflow
01085 if (L > 0) {
01086 #endif
01087 #endif
01088 #ifdef IBM
01089 L |= Exp_msk1 >> 4;
01090 #endif
01091 word0(a) = L;
01092 word1(a) = 0;
01093 #ifndef Avoid_Underflow
01094 #ifndef Sudden_Underflow
01095 }
01096 else {
01097 L = -L >> Exp_shift;
01098 if (L < Exp_shift) {
01099 word0(a) = 0x80000 >> L;
01100 word1(a) = 0;
01101 }
01102 else {
01103 word0(a) = 0;
01104 L -= Exp_shift;
01105 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
01106 }
01107 }
01108 #endif
01109 #endif
01110 return dval(a);
01111 }
01112
01113 static double
01114 b2d
01115 #ifdef KR_headers
01116 (a, e) Bigint *a; int *e;
01117 #else
01118 (Bigint *a, int *e)
01119 #endif
01120 {
01121 ULong *xa, *xa0, w, y, z;
01122 int k;
01123 double d;
01124 #ifdef VAX
01125 ULong d0, d1;
01126 #else
01127 #define d0 word0(d)
01128 #define d1 word1(d)
01129 #endif
01130
01131 xa0 = a->x;
01132 xa = xa0 + a->wds;
01133 y = *--xa;
01134 #ifdef DEBUG
01135 if (!y) Bug("zero y in b2d");
01136 #endif
01137 k = hi0bits(y);
01138 *e = 32 - k;
01139 #ifdef Pack_32
01140 if (k < Ebits) {
01141 d0 = Exp_1 | y >> (Ebits - k);
01142 w = xa > xa0 ? *--xa : 0;
01143 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
01144 goto ret_d;
01145 }
01146 z = xa > xa0 ? *--xa : 0;
01147 if (k -= Ebits) {
01148 d0 = Exp_1 | y << k | z >> (32 - k);
01149 y = xa > xa0 ? *--xa : 0;
01150 d1 = z << k | y >> (32 - k);
01151 }
01152 else {
01153 d0 = Exp_1 | y;
01154 d1 = z;
01155 }
01156 #else
01157 if (k < Ebits + 16) {
01158 z = xa > xa0 ? *--xa : 0;
01159 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
01160 w = xa > xa0 ? *--xa : 0;
01161 y = xa > xa0 ? *--xa : 0;
01162 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
01163 goto ret_d;
01164 }
01165 z = xa > xa0 ? *--xa : 0;
01166 w = xa > xa0 ? *--xa : 0;
01167 k -= Ebits + 16;
01168 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
01169 y = xa > xa0 ? *--xa : 0;
01170 d1 = w << k + 16 | y << k;
01171 #endif
01172 ret_d:
01173 #ifdef VAX
01174 word0(d) = d0 >> 16 | d0 << 16;
01175 word1(d) = d1 >> 16 | d1 << 16;
01176 #else
01177 #undef d0
01178 #undef d1
01179 #endif
01180 return dval(d);
01181 }
01182
01183 static Bigint *
01184 d2b
01185 #ifdef KR_headers
01186 (d, e, bits) double d; int *e, *bits;
01187 #else
01188 (double d, int *e, int *bits)
01189 #endif
01190 {
01191 Bigint *b;
01192 int de, k;
01193 ULong *x, y, z;
01194 #ifndef Sudden_Underflow
01195 int i;
01196 #endif
01197 #ifdef VAX
01198 ULong d0, d1;
01199 d0 = word0(d) >> 16 | word0(d) << 16;
01200 d1 = word1(d) >> 16 | word1(d) << 16;
01201 #else
01202 #define d0 word0(d)
01203 #define d1 word1(d)
01204 #endif
01205
01206 #ifdef Pack_32
01207 b = Balloc(1);
01208 #else
01209 b = Balloc(2);
01210 #endif
01211 x = b->x;
01212
01213 z = d0 & Frac_mask;
01214 d0 &= 0x7fffffff;
01215 #ifdef Sudden_Underflow
01216 de = (int)(d0 >> Exp_shift);
01217 #ifndef IBM
01218 z |= Exp_msk11;
01219 #endif
01220 #else
01221 if ((de = (int)(d0 >> Exp_shift)))
01222 z |= Exp_msk1;
01223 #endif
01224 #ifdef Pack_32
01225 if ((y = d1)) {
01226 if ((k = lo0bits(&y))) {
01227 x[0] = y | z << (32 - k);
01228 z >>= k;
01229 }
01230 else
01231 x[0] = y;
01232 #ifndef Sudden_Underflow
01233 i =
01234 #endif
01235 b->wds = (x[1] = z) ? 2 : 1;
01236 }
01237 else {
01238 #ifdef DEBUG
01239 if (!z)
01240 Bug("Zero passed to d2b");
01241 #endif
01242 k = lo0bits(&z);
01243 x[0] = z;
01244 #ifndef Sudden_Underflow
01245 i =
01246 #endif
01247 b->wds = 1;
01248 k += 32;
01249 }
01250 #else
01251 if (y = d1) {
01252 if (k = lo0bits(&y))
01253 if (k >= 16) {
01254 x[0] = y | z << 32 - k & 0xffff;
01255 x[1] = z >> k - 16 & 0xffff;
01256 x[2] = z >> k;
01257 i = 2;
01258 }
01259 else {
01260 x[0] = y & 0xffff;
01261 x[1] = y >> 16 | z << 16 - k & 0xffff;
01262 x[2] = z >> k & 0xffff;
01263 x[3] = z >> k+16;
01264 i = 3;
01265 }
01266 else {
01267 x[0] = y & 0xffff;
01268 x[1] = y >> 16;
01269 x[2] = z & 0xffff;
01270 x[3] = z >> 16;
01271 i = 3;
01272 }
01273 }
01274 else {
01275 #ifdef DEBUG
01276 if (!z)
01277 Bug("Zero passed to d2b");
01278 #endif
01279 k = lo0bits(&z);
01280 if (k >= 16) {
01281 x[0] = z;
01282 i = 0;
01283 }
01284 else {
01285 x[0] = z & 0xffff;
01286 x[1] = z >> 16;
01287 i = 1;
01288 }
01289 k += 32;
01290 }
01291 while(!x[i])
01292 --i;
01293 b->wds = i + 1;
01294 #endif
01295 #ifndef Sudden_Underflow
01296 if (de) {
01297 #endif
01298 #ifdef IBM
01299 *e = (de - Bias - (P-1) << 2) + k;
01300 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
01301 #else
01302 *e = de - Bias - (P-1) + k;
01303 *bits = P - k;
01304 #endif
01305 #ifndef Sudden_Underflow
01306 }
01307 else {
01308 *e = de - Bias - (P-1) + 1 + k;
01309 #ifdef Pack_32
01310 *bits = 32*i - hi0bits(x[i-1]);
01311 #else
01312 *bits = (i+2)*16 - hi0bits(x[i]);
01313 #endif
01314 }
01315 #endif
01316 return b;
01317 }
01318 #undef d0
01319 #undef d1
01320
01321 static double
01322 ratio
01323 #ifdef KR_headers
01324 (a, b) Bigint *a, *b;
01325 #else
01326 (Bigint *a, Bigint *b)
01327 #endif
01328 {
01329 double da, db;
01330 int k, ka, kb;
01331
01332 dval(da) = b2d(a, &ka);
01333 dval(db) = b2d(b, &kb);
01334 #ifdef Pack_32
01335 k = ka - kb + 32*(a->wds - b->wds);
01336 #else
01337 k = ka - kb + 16*(a->wds - b->wds);
01338 #endif
01339 #ifdef IBM
01340 if (k > 0) {
01341 word0(da) += (k >> 2)*Exp_msk1;
01342 if (k &= 3)
01343 dval(da) *= 1 << k;
01344 }
01345 else {
01346 k = -k;
01347 word0(db) += (k >> 2)*Exp_msk1;
01348 if (k &= 3)
01349 dval(db) *= 1 << k;
01350 }
01351 #else
01352 if (k > 0)
01353 word0(da) += k*Exp_msk1;
01354 else {
01355 k = -k;
01356 word0(db) += k*Exp_msk1;
01357 }
01358 #endif
01359 return dval(da) / dval(db);
01360 }
01361
01362 static CONST double
01363 tens[] = {
01364 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
01365 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
01366 1e20, 1e21, 1e22
01367 #ifdef VAX
01368 , 1e23, 1e24
01369 #endif
01370 };
01371
01372 static CONST double
01373 #ifdef IEEE_Arith
01374 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
01375 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
01376 #ifdef Avoid_Underflow
01377 9007199254740992.*9007199254740992.e-256
01378
01379 #else
01380 1e-256
01381 #endif
01382 };
01383
01384
01385 #define Scale_Bit 0x10
01386 #define n_bigtens 5
01387 #else
01388 #ifdef IBM
01389 bigtens[] = { 1e16, 1e32, 1e64 };
01390 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
01391 #define n_bigtens 3
01392 #else
01393 bigtens[] = { 1e16, 1e32 };
01394 static CONST double tinytens[] = { 1e-16, 1e-32 };
01395 #define n_bigtens 2
01396 #endif
01397 #endif
01398
01399 #ifdef INFNAN_CHECK
01400
01401 #ifndef NAN_WORD0
01402 #define NAN_WORD0 0x7ff80000
01403 #endif
01404
01405 #ifndef NAN_WORD1
01406 #define NAN_WORD1 0
01407 #endif
01408
01409 static int
01410 match
01411 #ifdef KR_headers
01412 (sp, t) char **sp, *t;
01413 #else
01414 (CONST char **sp, char *t)
01415 #endif
01416 {
01417 int c, d;
01418 CONST char *s = *sp;
01419
01420 while((d = *t++)) {
01421 if ((c = *++s) >= 'A' && c <= 'Z')
01422 c += 'a' - 'A';
01423 if (c != d)
01424 return 0;
01425 }
01426 *sp = s + 1;
01427 return 1;
01428 }
01429
01430 #ifndef No_Hex_NaN
01431 static void
01432 hexnan
01433 #ifdef KR_headers
01434 (rvp, sp) double *rvp; CONST char **sp;
01435 #else
01436 (double *rvp, CONST char **sp)
01437 #endif
01438 {
01439 ULong c, x[2];
01440 CONST char *s;
01441 int havedig, udx0, xshift;
01442
01443 x[0] = x[1] = 0;
01444 havedig = xshift = 0;
01445 udx0 = 1;
01446 s = *sp;
01447
01448 while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
01449 ++s;
01450 if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
01451 s += 2;
01452 while((c = *(CONST unsigned char*)++s)) {
01453 if (c >= '0' && c <= '9')
01454 c -= '0';
01455 else if (c >= 'a' && c <= 'f')
01456 c += 10 - 'a';
01457 else if (c >= 'A' && c <= 'F')
01458 c += 10 - 'A';
01459 else if (c <= ' ') {
01460 if (udx0 && havedig) {
01461 udx0 = 0;
01462 xshift = 1;
01463 }
01464 continue;
01465 }
01466 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
01467 else if ( c == ')' && havedig) {
01468 *sp = s + 1;
01469 break;
01470 }
01471 else
01472 return;
01473 #else
01474 else {
01475 do {
01476 if ( c == ')') {
01477 *sp = s + 1;
01478 break;
01479 }
01480 } while((c = *++s));
01481 break;
01482 }
01483 #endif
01484 havedig = 1;
01485 if (xshift) {
01486 xshift = 0;
01487 x[0] = x[1];
01488 x[1] = 0;
01489 }
01490 if (udx0)
01491 x[0] = (x[0] << 4) | (x[1] >> 28);
01492 x[1] = (x[1] << 4) | c;
01493 }
01494 if ((x[0] &= 0xfffff) || x[1]) {
01495 word0(*rvp) = Exp_mask | x[0];
01496 word1(*rvp) = x[1];
01497 }
01498 }
01499 #endif
01500 #endif
01501
01502 double
01503 sb_strtod
01504 #ifdef KR_headers
01505 (s00, se) CONST char *s00; char **se;
01506 #else
01507 (CONST char *s00, char **se)
01508 #endif
01509 {
01510 #ifdef Avoid_Underflow
01511 int scale;
01512 #endif
01513 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
01514 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
01515 CONST char *s, *s0, *s1;
01516 double aadj, aadj1, adj, rv, rv0;
01517 Long L;
01518 ULong y, z;
01519 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
01520 #ifdef SET_INEXACT
01521 int inexact, oldinexact;
01522 #endif
01523 #ifdef Honor_FLT_ROUNDS
01524 int rounding;
01525 #endif
01526 #ifdef USE_LOCALE
01527 CONST char *s2;
01528 #endif
01529
01530 sign = nz0 = nz = 0;
01531 dval(rv) = 0.;
01532 for(s = s00;;s++) switch(*s) {
01533 case '-':
01534 sign = 1;
01535
01536 case '+':
01537 if (*++s)
01538 goto break2;
01539
01540 case 0:
01541 goto ret0;
01542 case '\t':
01543 case '\n':
01544 case '\v':
01545 case '\f':
01546 case '\r':
01547 case ' ':
01548 continue;
01549 default:
01550 goto break2;
01551 }
01552 break2:
01553 if (*s == '0') {
01554 nz0 = 1;
01555 while(*++s == '0') ;
01556 if (!*s)
01557 goto ret;
01558 }
01559 s0 = s;
01560 y = z = 0;
01561 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
01562 if (nd < 9)
01563 y = 10*y + c - '0';
01564 else if (nd < 16)
01565 z = 10*z + c - '0';
01566 nd0 = nd;
01567 #ifdef USE_LOCALE
01568 s1 = localeconv()->decimal_point;
01569 if (c == *s1) {
01570 c = '.';
01571 if (*++s1) {
01572 s2 = s;
01573 for(;;) {
01574 if (*++s2 != *s1) {
01575 c = 0;
01576 break;
01577 }
01578 if (!*++s1) {
01579 s = s2;
01580 break;
01581 }
01582 }
01583 }
01584 }
01585 #endif
01586 if (c == '.') {
01587 c = *++s;
01588 if (!nd) {
01589 for(; c == '0'; c = *++s)
01590 nz++;
01591 if (c > '0' && c <= '9') {
01592 s0 = s;
01593 nf += nz;
01594 nz = 0;
01595 goto have_dig;
01596 }
01597 goto dig_done;
01598 }
01599 for(; c >= '0' && c <= '9'; c = *++s) {
01600 have_dig:
01601 nz++;
01602 if (c -= '0') {
01603 nf += nz;
01604 for(i = 1; i < nz; i++)
01605 if (nd++ < 9)
01606 y *= 10;
01607 else if (nd <= DBL_DIG + 1)
01608 z *= 10;
01609 if (nd++ < 9)
01610 y = 10*y + c;
01611 else if (nd <= DBL_DIG + 1)
01612 z = 10*z + c;
01613 nz = 0;
01614 }
01615 }
01616 }
01617 dig_done:
01618 e = 0;
01619 if (c == 'e' || c == 'E') {
01620 if (!nd && !nz && !nz0) {
01621 goto ret0;
01622 }
01623 s00 = s;
01624 esign = 0;
01625 switch(c = *++s) {
01626 case '-':
01627 esign = 1;
01628 case '+':
01629 c = *++s;
01630 }
01631 if (c >= '0' && c <= '9') {
01632 while(c == '0')
01633 c = *++s;
01634 if (c > '0' && c <= '9') {
01635 L = c - '0';
01636 s1 = s;
01637 while((c = *++s) >= '0' && c <= '9')
01638 L = 10*L + c - '0';
01639 if (s - s1 > 8 || L > 19999)
01640
01641
01642
01643 e = 19999;
01644 else
01645 e = (int)L;
01646 if (esign)
01647 e = -e;
01648 }
01649 else
01650 e = 0;
01651 }
01652 else
01653 s = s00;
01654 }
01655 if (!nd) {
01656 if (!nz && !nz0) {
01657 #ifdef INFNAN_CHECK
01658
01659 switch(c) {
01660 case 'i':
01661 case 'I':
01662 if (match(&s,"nf")) {
01663 --s;
01664 if (!match(&s,"inity"))
01665 ++s;
01666 word0(rv) = 0x7ff00000;
01667 word1(rv) = 0;
01668 goto ret;
01669 }
01670 break;
01671 case 'n':
01672 case 'N':
01673 if (match(&s, "an")) {
01674 word0(rv) = NAN_WORD0;
01675 word1(rv) = NAN_WORD1;
01676 #ifndef No_Hex_NaN
01677 if (*s == '(')
01678 hexnan(&rv, &s);
01679 #endif
01680 goto ret;
01681 }
01682 }
01683 #endif
01684 ret0:
01685 s = s00;
01686 sign = 0;
01687 }
01688 goto ret;
01689 }
01690 e1 = e -= nf;
01691
01692
01693
01694
01695
01696
01697 if (!nd0)
01698 nd0 = nd;
01699 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
01700 dval(rv) = y;
01701 if (k > 9) {
01702 #ifdef SET_INEXACT
01703 if (k > DBL_DIG)
01704 oldinexact = get_inexact();
01705 #endif
01706 dval(rv) = tens[k - 9] * dval(rv) + z;
01707 }
01708 bd0 = 0;
01709 if (nd <= DBL_DIG
01710 #ifndef RND_PRODQUOT
01711 #ifndef Honor_FLT_ROUNDS
01712 && Flt_Rounds == 1
01713 #endif
01714 #endif
01715 ) {
01716 if (!e)
01717 goto ret;
01718 if (e > 0) {
01719 if (e <= Ten_pmax) {
01720 #ifdef VAX
01721 goto vax_ovfl_check;
01722 #else
01723 #ifdef Honor_FLT_ROUNDS
01724
01725 if (sign) {
01726 rv = -rv;
01727 sign = 0;
01728 }
01729 #endif
01730 rounded_product(dval(rv), tens[e]);
01731 goto ret;
01732 #endif
01733 }
01734 i = DBL_DIG - nd;
01735 if (e <= Ten_pmax + i) {
01736
01737
01738
01739 #ifdef Honor_FLT_ROUNDS
01740
01741 if (sign) {
01742 rv = -rv;
01743 sign = 0;
01744 }
01745 #endif
01746 e -= i;
01747 dval(rv) *= tens[i];
01748 #ifdef VAX
01749
01750
01751
01752 vax_ovfl_check:
01753 word0(rv) -= P*Exp_msk1;
01754 rounded_product(dval(rv), tens[e]);
01755 if ((word0(rv) & Exp_mask)
01756 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
01757 goto ovfl;
01758 word0(rv) += P*Exp_msk1;
01759 #else
01760 rounded_product(dval(rv), tens[e]);
01761 #endif
01762 goto ret;
01763 }
01764 }
01765 #ifndef Inaccurate_Divide
01766 else if (e >= -Ten_pmax) {
01767 #ifdef Honor_FLT_ROUNDS
01768
01769 if (sign) {
01770 rv = -rv;
01771 sign = 0;
01772 }
01773 #endif
01774 rounded_quotient(dval(rv), tens[-e]);
01775 goto ret;
01776 }
01777 #endif
01778 }
01779 e1 += nd - k;
01780
01781 #ifdef IEEE_Arith
01782 #ifdef SET_INEXACT
01783 inexact = 1;
01784 if (k <= DBL_DIG)
01785 oldinexact = get_inexact();
01786 #endif
01787 #ifdef Avoid_Underflow
01788 scale = 0;
01789 #endif
01790 #ifdef Honor_FLT_ROUNDS
01791 if ((rounding = Flt_Rounds) >= 2) {
01792 if (sign)
01793 rounding = rounding == 2 ? 0 : 2;
01794 else
01795 if (rounding != 2)
01796 rounding = 0;
01797 }
01798 #endif
01799 #endif
01800
01801
01802
01803 if (e1 > 0) {
01804 if ((i = e1 & 15))
01805 dval(rv) *= tens[i];
01806 if (e1 &= ~15) {
01807 if (e1 > DBL_MAX_10_EXP) {
01808 ovfl:
01809 #ifndef NO_ERRNO
01810 errno = ERANGE;
01811 #endif
01812
01813 #ifdef IEEE_Arith
01814 #ifdef Honor_FLT_ROUNDS
01815 switch(rounding) {
01816 case 0:
01817 case 3:
01818 word0(rv) = Big0;
01819 word1(rv) = Big1;
01820 break;
01821 default:
01822 word0(rv) = Exp_mask;
01823 word1(rv) = 0;
01824 }
01825 #else
01826 word0(rv) = Exp_mask;
01827 word1(rv) = 0;
01828 #endif
01829 #ifdef SET_INEXACT
01830
01831 dval(rv0) = 1e300;
01832 dval(rv0) *= dval(rv0);
01833 #endif
01834 #else
01835 word0(rv) = Big0;
01836 word1(rv) = Big1;
01837 #endif
01838 if (bd0)
01839 goto retfree;
01840 goto ret;
01841 }
01842 e1 >>= 4;
01843 for(j = 0; e1 > 1; j++, e1 >>= 1)
01844 if (e1 & 1)
01845 dval(rv) *= bigtens[j];
01846
01847 word0(rv) -= P*Exp_msk1;
01848 dval(rv) *= bigtens[j];
01849 if ((z = word0(rv) & Exp_mask)
01850 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
01851 goto ovfl;
01852 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
01853
01854
01855 word0(rv) = Big0;
01856 word1(rv) = Big1;
01857 }
01858 else
01859 word0(rv) += P*Exp_msk1;
01860 }
01861 }
01862 else if (e1 < 0) {
01863 e1 = -e1;
01864 if ((i = e1 & 15))
01865 dval(rv) /= tens[i];
01866 if (e1 >>= 4) {
01867 if (e1 >= 1 << n_bigtens)
01868 goto undfl;
01869 #ifdef Avoid_Underflow
01870 if (e1 & Scale_Bit)
01871 scale = 2*P;
01872 for(j = 0; e1 > 0; j++, e1 >>= 1)
01873 if (e1 & 1)
01874 dval(rv) *= tinytens[j];
01875 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
01876 >> Exp_shift)) > 0) {
01877
01878 if (j >= 32) {
01879 word1(rv) = 0;
01880 if (j >= 53)
01881 word0(rv) = (P+2)*Exp_msk1;
01882 else
01883 word0(rv) &= 0xffffffff << (j-32);
01884 }
01885 else
01886 word1(rv) &= 0xffffffff << j;
01887 }
01888 #else
01889 for(j = 0; e1 > 1; j++, e1 >>= 1)
01890 if (e1 & 1)
01891 dval(rv) *= tinytens[j];
01892
01893 dval(rv0) = dval(rv);
01894 dval(rv) *= tinytens[j];
01895 if (!dval(rv)) {
01896 dval(rv) = 2.*dval(rv0);
01897 dval(rv) *= tinytens[j];
01898 #endif
01899 if (!dval(rv)) {
01900 undfl:
01901 dval(rv) = 0.;
01902 #ifndef NO_ERRNO
01903 errno = ERANGE;
01904 #endif
01905 if (bd0)
01906 goto retfree;
01907 goto ret;
01908 }
01909 #ifndef Avoid_Underflow
01910 word0(rv) = Tiny0;
01911 word1(rv) = Tiny1;
01912
01913
01914
01915 }
01916 #endif
01917 }
01918 }
01919
01920
01921
01922
01923
01924 bd0 = s2b(s0, nd0, nd, y);
01925
01926 for(;;) {
01927 bd = Balloc(bd0->k);
01928 Bcopy(bd, bd0);
01929 bb = d2b(dval(rv), &bbe, &bbbits);
01930 bs = i2b(1);
01931
01932 if (e >= 0) {
01933 bb2 = bb5 = 0;
01934 bd2 = bd5 = e;
01935 }
01936 else {
01937 bb2 = bb5 = -e;
01938 bd2 = bd5 = 0;
01939 }
01940 if (bbe >= 0)
01941 bb2 += bbe;
01942 else
01943 bd2 -= bbe;
01944 bs2 = bb2;
01945 #ifdef Honor_FLT_ROUNDS
01946 if (rounding != 1)
01947 bs2++;
01948 #endif
01949 #ifdef Avoid_Underflow
01950 j = bbe - scale;
01951 i = j + bbbits - 1;
01952 if (i < Emin)
01953 j += P - Emin;
01954 else
01955 j = P + 1 - bbbits;
01956 #else
01957 #ifdef Sudden_Underflow
01958 #ifdef IBM
01959 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
01960 #else
01961 j = P + 1 - bbbits;
01962 #endif
01963 #else
01964 j = bbe;
01965 i = j + bbbits - 1;
01966 if (i < Emin)
01967 j += P - Emin;
01968 else
01969 j = P + 1 - bbbits;
01970 #endif
01971 #endif
01972 bb2 += j;
01973 bd2 += j;
01974 #ifdef Avoid_Underflow
01975 bd2 += scale;
01976 #endif
01977 i = bb2 < bd2 ? bb2 : bd2;
01978 if (i > bs2)
01979 i = bs2;
01980 if (i > 0) {
01981 bb2 -= i;
01982 bd2 -= i;
01983 bs2 -= i;
01984 }
01985 if (bb5 > 0) {
01986 bs = pow5mult(bs, bb5);
01987 bb1 = mult(bs, bb);
01988 Bfree(bb);
01989 bb = bb1;
01990 }
01991 if (bb2 > 0)
01992 bb = lshift(bb, bb2);
01993 if (bd5 > 0)
01994 bd = pow5mult(bd, bd5);
01995 if (bd2 > 0)
01996 bd = lshift(bd, bd2);
01997 if (bs2 > 0)
01998 bs = lshift(bs, bs2);
01999 delta = diff(bb, bd);
02000 dsign = delta->sign;
02001 delta->sign = 0;
02002 i = cmp(delta, bs);
02003 #ifdef Honor_FLT_ROUNDS
02004 if (rounding != 1) {
02005 if (i < 0) {
02006
02007 if (!delta->x[0] && delta->wds <= 1) {
02008
02009 #ifdef SET_INEXACT
02010 inexact = 0;
02011 #endif
02012 break;
02013 }
02014 if (rounding) {
02015 if (dsign) {
02016 adj = 1.;
02017 goto apply_adj;
02018 }
02019 }
02020 else if (!dsign) {
02021 adj = -1.;
02022 if (!word1(rv)
02023 && !(word0(rv) & Frac_mask)) {
02024 y = word0(rv) & Exp_mask;
02025 #ifdef Avoid_Underflow
02026 if (!scale || y > 2*P*Exp_msk1)
02027 #else
02028 if (y)
02029 #endif
02030 {
02031 delta = lshift(delta,Log2P);
02032 if (cmp(delta, bs) <= 0)
02033 adj = -0.5;
02034 }
02035 }
02036 apply_adj:
02037 #ifdef Avoid_Underflow
02038 if (scale && (y = word0(rv) & Exp_mask)
02039 <= 2*P*Exp_msk1)
02040 word0(adj) += (2*P+1)*Exp_msk1 - y;
02041 #else
02042 #ifdef Sudden_Underflow
02043 if ((word0(rv) & Exp_mask) <=
02044 P*Exp_msk1) {
02045 word0(rv) += P*Exp_msk1;
02046 dval(rv) += adj*ulp(dval(rv));
02047 word0(rv) -= P*Exp_msk1;
02048 }
02049 else
02050 #endif
02051 #endif
02052 dval(rv) += adj*ulp(dval(rv));
02053 }
02054 break;
02055 }
02056 adj = ratio(delta, bs);
02057 if (adj < 1.)
02058 adj = 1.;
02059 if (adj <= 0x7ffffffe) {
02060
02061 y = adj;
02062 if (y != adj) {
02063 if (!((rounding>>1) ^ dsign))
02064 y++;
02065 adj = y;
02066 }
02067 }
02068 #ifdef Avoid_Underflow
02069 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02070 word0(adj) += (2*P+1)*Exp_msk1 - y;
02071 #else
02072 #ifdef Sudden_Underflow
02073 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02074 word0(rv) += P*Exp_msk1;
02075 adj *= ulp(dval(rv));
02076 if (dsign)
02077 dval(rv) += adj;
02078 else
02079 dval(rv) -= adj;
02080 word0(rv) -= P*Exp_msk1;
02081 goto cont;
02082 }
02083 #endif
02084 #endif
02085 adj *= ulp(dval(rv));
02086 if (dsign)
02087 dval(rv) += adj;
02088 else
02089 dval(rv) -= adj;
02090 goto cont;
02091 }
02092 #endif
02093
02094 if (i < 0) {
02095
02096
02097
02098 if (dsign || word1(rv) || word0(rv) & Bndry_mask
02099 #ifdef IEEE_Arith
02100 #ifdef Avoid_Underflow
02101 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
02102 #else
02103 || (word0(rv) & Exp_mask) <= Exp_msk1
02104 #endif
02105 #endif
02106 ) {
02107 #ifdef SET_INEXACT
02108 if (!delta->x[0] && delta->wds <= 1)
02109 inexact = 0;
02110 #endif
02111 break;
02112 }
02113 if (!delta->x[0] && delta->wds <= 1) {
02114
02115 #ifdef SET_INEXACT
02116 inexact = 0;
02117 #endif
02118 break;
02119 }
02120 delta = lshift(delta,Log2P);
02121 if (cmp(delta, bs) > 0)
02122 goto drop_down;
02123 break;
02124 }
02125 if (i == 0) {
02126
02127 if (dsign) {
02128 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
02129 && word1(rv) == (
02130 #ifdef Avoid_Underflow
02131 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02132 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
02133 #endif
02134 0xffffffff)) {
02135
02136 word0(rv) = (word0(rv) & Exp_mask)
02137 + Exp_msk1
02138 #ifdef IBM
02139 | Exp_msk1 >> 4
02140 #endif
02141 ;
02142 word1(rv) = 0;
02143 #ifdef Avoid_Underflow
02144 dsign = 0;
02145 #endif
02146 break;
02147 }
02148 }
02149 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
02150 drop_down:
02151
02152 #ifdef Sudden_Underflow
02153 L = word0(rv) & Exp_mask;
02154 #ifdef IBM
02155 if (L < Exp_msk1)
02156 #else
02157 #ifdef Avoid_Underflow
02158 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
02159 #else
02160 if (L <= Exp_msk1)
02161 #endif
02162 #endif
02163 goto undfl;
02164 L -= Exp_msk1;
02165 #else
02166 #ifdef Avoid_Underflow
02167 if (scale) {
02168 L = word0(rv) & Exp_mask;
02169 if (L <= (2*P+1)*Exp_msk1) {
02170 if (L > (P+2)*Exp_msk1)
02171
02172
02173 break;
02174
02175 goto undfl;
02176 }
02177 }
02178 #endif
02179 L = (word0(rv) & Exp_mask) - Exp_msk1;
02180 #endif
02181 word0(rv) = L | Bndry_mask1;
02182 word1(rv) = 0xffffffff;
02183 #ifdef IBM
02184 goto cont;
02185 #else
02186 break;
02187 #endif
02188 }
02189 #ifndef ROUND_BIASED
02190 if (!(word1(rv) & LSB))
02191 break;
02192 #endif
02193 if (dsign)
02194 dval(rv) += ulp(dval(rv));
02195 #ifndef ROUND_BIASED
02196 else {
02197 dval(rv) -= ulp(dval(rv));
02198 #ifndef Sudden_Underflow
02199 if (!dval(rv))
02200 goto undfl;
02201 #endif
02202 }
02203 #ifdef Avoid_Underflow
02204 dsign = 1 - dsign;
02205 #endif
02206 #endif
02207 break;
02208 }
02209 if ((aadj = ratio(delta, bs)) <= 2.) {
02210 if (dsign)
02211 aadj = aadj1 = 1.;
02212 else if (word1(rv) || word0(rv) & Bndry_mask) {
02213 #ifndef Sudden_Underflow
02214 if (word1(rv) == Tiny1 && !word0(rv))
02215 goto undfl;
02216 #endif
02217 aadj = 1.;
02218 aadj1 = -1.;
02219 }
02220 else {
02221
02222
02223
02224 if (aadj < 2./FLT_RADIX)
02225 aadj = 1./FLT_RADIX;
02226 else
02227 aadj *= 0.5;
02228 aadj1 = -aadj;
02229 }
02230 }
02231 else {
02232 aadj *= 0.5;
02233 aadj1 = dsign ? aadj : -aadj;
02234 #ifdef Check_FLT_ROUNDS
02235 switch(Rounding) {
02236 case 2:
02237 aadj1 -= 0.5;
02238 break;
02239 case 0:
02240 case 3:
02241 aadj1 += 0.5;
02242 }
02243 #else
02244 if (Flt_Rounds == 0)
02245 aadj1 += 0.5;
02246 #endif
02247 }
02248 y = word0(rv) & Exp_mask;
02249
02250
02251
02252 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
02253 dval(rv0) = dval(rv);
02254 word0(rv) -= P*Exp_msk1;
02255 adj = aadj1 * ulp(dval(rv));
02256 dval(rv) += adj;
02257 if ((word0(rv) & Exp_mask) >=
02258 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
02259 if (word0(rv0) == Big0 && word1(rv0) == Big1)
02260 goto ovfl;
02261 word0(rv) = Big0;
02262 word1(rv) = Big1;
02263 goto cont;
02264 }
02265 else
02266 word0(rv) += P*Exp_msk1;
02267 }
02268 else {
02269 #ifdef Avoid_Underflow
02270 if (scale && y <= 2*P*Exp_msk1) {
02271 if (aadj <= 0x7fffffff) {
02272 if ((z = (uint32)aadj) <= 0)
02273 z = 1;
02274 aadj = z;
02275 aadj1 = dsign ? aadj : -aadj;
02276 }
02277 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
02278 }
02279 adj = aadj1 * ulp(dval(rv));
02280 dval(rv) += adj;
02281 #else
02282 #ifdef Sudden_Underflow
02283 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02284 dval(rv0) = dval(rv);
02285 word0(rv) += P*Exp_msk1;
02286 adj = aadj1 * ulp(dval(rv));
02287 dval(rv) += adj;
02288 #ifdef IBM
02289 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
02290 #else
02291 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
02292 #endif
02293 {
02294 if (word0(rv0) == Tiny0
02295 && word1(rv0) == Tiny1)
02296 goto undfl;
02297 word0(rv) = Tiny0;
02298 word1(rv) = Tiny1;
02299 goto cont;
02300 }
02301 else
02302 word0(rv) -= P*Exp_msk1;
02303 }
02304 else {
02305 adj = aadj1 * ulp(dval(rv));
02306 dval(rv) += adj;
02307 }
02308 #else
02309
02310
02311
02312
02313
02314
02315
02316 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
02317 aadj1 = (double)(int)(aadj + 0.5);
02318 if (!dsign)
02319 aadj1 = -aadj1;
02320 }
02321 adj = aadj1 * ulp(dval(rv));
02322 dval(rv) += adj;
02323 #endif
02324 #endif
02325 }
02326 z = word0(rv) & Exp_mask;
02327 #ifndef SET_INEXACT
02328 #ifdef Avoid_Underflow
02329 if (!scale)
02330 #endif
02331 if (y == z) {
02332
02333 L = (Long)aadj;
02334 aadj -= L;
02335
02336 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
02337 if (aadj < .4999999 || aadj > .5000001)
02338 break;
02339 }
02340 else if (aadj < .4999999/FLT_RADIX)
02341 break;
02342 }
02343 #endif
02344 cont:
02345 Bfree(bb);
02346 Bfree(bd);
02347 Bfree(bs);
02348 Bfree(delta);
02349 }
02350 #ifdef SET_INEXACT
02351 if (inexact) {
02352 if (!oldinexact) {
02353 word0(rv0) = Exp_1 + (70 << Exp_shift);
02354 word1(rv0) = 0;
02355 dval(rv0) += 1.;
02356 }
02357 }
02358 else if (!oldinexact)
02359 clear_inexact();
02360 #endif
02361 #ifdef Avoid_Underflow
02362 if (scale) {
02363 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
02364 word1(rv0) = 0;
02365 dval(rv) *= dval(rv0);
02366 #ifndef NO_ERRNO
02367
02368 if (word0(rv) == 0 && word1(rv) == 0)
02369 errno = ERANGE;
02370 #endif
02371 }
02372 #endif
02373 #ifdef SET_INEXACT
02374 if (inexact && !(word0(rv) & Exp_mask)) {
02375
02376 dval(rv0) = 1e-300;
02377 dval(rv0) *= dval(rv0);
02378 }
02379 #endif
02380 retfree:
02381 Bfree(bb);
02382 Bfree(bd);
02383 Bfree(bs);
02384 Bfree(bd0);
02385 Bfree(delta);
02386 ret:
02387 if (se)
02388 *se = (char *)s;
02389 return sign ? -dval(rv) : dval(rv);
02390 }