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 U 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 U 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 U d;
01195 #ifndef Sudden_Underflow
01196 int i;
01197 #endif
01198 #ifdef VAX
01199 ULong d0, d1;
01200 d0 = word0(d) >> 16 | word0(d) << 16;
01201 d1 = word1(d) >> 16 | word1(d) << 16;
01202 #else
01203 #define d0 word0(d)
01204 #define d1 word1(d)
01205 #endif
01206 dval(d) = _d;
01207
01208 #ifdef Pack_32
01209 b = Balloc(1);
01210 #else
01211 b = Balloc(2);
01212 #endif
01213 x = b->x;
01214
01215 z = d0 & Frac_mask;
01216 d0 &= 0x7fffffff;
01217 #ifdef Sudden_Underflow
01218 de = (int)(d0 >> Exp_shift);
01219 #ifndef IBM
01220 z |= Exp_msk11;
01221 #endif
01222 #else
01223 if ((de = (int)(d0 >> Exp_shift)))
01224 z |= Exp_msk1;
01225 #endif
01226 #ifdef Pack_32
01227 if ((y = d1)) {
01228 if ((k = lo0bits(&y))) {
01229 x[0] = y | z << (32 - k);
01230 z >>= k;
01231 }
01232 else
01233 x[0] = y;
01234 #ifndef Sudden_Underflow
01235 i =
01236 #endif
01237 b->wds = (x[1] = z) ? 2 : 1;
01238 }
01239 else {
01240 #ifdef DEBUG
01241 if (!z)
01242 Bug("Zero passed to d2b");
01243 #endif
01244 k = lo0bits(&z);
01245 x[0] = z;
01246 #ifndef Sudden_Underflow
01247 i =
01248 #endif
01249 b->wds = 1;
01250 k += 32;
01251 }
01252 #else
01253 if (y = d1) {
01254 if (k = lo0bits(&y))
01255 if (k >= 16) {
01256 x[0] = y | z << 32 - k & 0xffff;
01257 x[1] = z >> k - 16 & 0xffff;
01258 x[2] = z >> k;
01259 i = 2;
01260 }
01261 else {
01262 x[0] = y & 0xffff;
01263 x[1] = y >> 16 | z << 16 - k & 0xffff;
01264 x[2] = z >> k & 0xffff;
01265 x[3] = z >> k+16;
01266 i = 3;
01267 }
01268 else {
01269 x[0] = y & 0xffff;
01270 x[1] = y >> 16;
01271 x[2] = z & 0xffff;
01272 x[3] = z >> 16;
01273 i = 3;
01274 }
01275 }
01276 else {
01277 #ifdef DEBUG
01278 if (!z)
01279 Bug("Zero passed to d2b");
01280 #endif
01281 k = lo0bits(&z);
01282 if (k >= 16) {
01283 x[0] = z;
01284 i = 0;
01285 }
01286 else {
01287 x[0] = z & 0xffff;
01288 x[1] = z >> 16;
01289 i = 1;
01290 }
01291 k += 32;
01292 }
01293 while(!x[i])
01294 --i;
01295 b->wds = i + 1;
01296 #endif
01297 #ifndef Sudden_Underflow
01298 if (de) {
01299 #endif
01300 #ifdef IBM
01301 *e = (de - Bias - (P-1) << 2) + k;
01302 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
01303 #else
01304 *e = de - Bias - (P-1) + k;
01305 *bits = P - k;
01306 #endif
01307 #ifndef Sudden_Underflow
01308 }
01309 else {
01310 *e = de - Bias - (P-1) + 1 + k;
01311 #ifdef Pack_32
01312 *bits = 32*i - hi0bits(x[i-1]);
01313 #else
01314 *bits = (i+2)*16 - hi0bits(x[i]);
01315 #endif
01316 }
01317 #endif
01318 return b;
01319 }
01320 #undef d0
01321 #undef d1
01322
01323 static double
01324 ratio
01325 #ifdef KR_headers
01326 (a, b) Bigint *a, *b;
01327 #else
01328 (Bigint *a, Bigint *b)
01329 #endif
01330 {
01331 U da, db;
01332 int k, ka, kb;
01333
01334 dval(da) = b2d(a, &ka);
01335 dval(db) = b2d(b, &kb);
01336 #ifdef Pack_32
01337 k = ka - kb + 32*(a->wds - b->wds);
01338 #else
01339 k = ka - kb + 16*(a->wds - b->wds);
01340 #endif
01341 #ifdef IBM
01342 if (k > 0) {
01343 word0(da) += (k >> 2)*Exp_msk1;
01344 if (k &= 3)
01345 dval(da) *= 1 << k;
01346 }
01347 else {
01348 k = -k;
01349 word0(db) += (k >> 2)*Exp_msk1;
01350 if (k &= 3)
01351 dval(db) *= 1 << k;
01352 }
01353 #else
01354 if (k > 0)
01355 word0(da) += k*Exp_msk1;
01356 else {
01357 k = -k;
01358 word0(db) += k*Exp_msk1;
01359 }
01360 #endif
01361 return dval(da) / dval(db);
01362 }
01363
01364 static CONST double
01365 tens[] = {
01366 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
01367 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
01368 1e20, 1e21, 1e22
01369 #ifdef VAX
01370 , 1e23, 1e24
01371 #endif
01372 };
01373
01374 static CONST double
01375 #ifdef IEEE_Arith
01376 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
01377 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
01378 #ifdef Avoid_Underflow
01379 9007199254740992.*9007199254740992.e-256
01380
01381 #else
01382 1e-256
01383 #endif
01384 };
01385
01386
01387 #define Scale_Bit 0x10
01388 #define n_bigtens 5
01389 #else
01390 #ifdef IBM
01391 bigtens[] = { 1e16, 1e32, 1e64 };
01392 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
01393 #define n_bigtens 3
01394 #else
01395 bigtens[] = { 1e16, 1e32 };
01396 static CONST double tinytens[] = { 1e-16, 1e-32 };
01397 #define n_bigtens 2
01398 #endif
01399 #endif
01400
01401 #ifdef INFNAN_CHECK
01402
01403 #ifndef NAN_WORD0
01404 #define NAN_WORD0 0x7ff80000
01405 #endif
01406
01407 #ifndef NAN_WORD1
01408 #define NAN_WORD1 0
01409 #endif
01410
01411 static int
01412 match
01413 #ifdef KR_headers
01414 (sp, t) char **sp, *t;
01415 #else
01416 (CONST char **sp, char *t)
01417 #endif
01418 {
01419 int c, d;
01420 CONST char *s = *sp;
01421
01422 while((d = *t++)) {
01423 if ((c = *++s) >= 'A' && c <= 'Z')
01424 c += 'a' - 'A';
01425 if (c != d)
01426 return 0;
01427 }
01428 *sp = s + 1;
01429 return 1;
01430 }
01431
01432 #ifndef No_Hex_NaN
01433 static void
01434 hexnan
01435 #ifdef KR_headers
01436 (rvp, sp) double *rvp; CONST char **sp;
01437 #else
01438 (U *rvp, CONST char **sp)
01439 #endif
01440 {
01441 ULong c, x[2];
01442 CONST char *s;
01443 int havedig, udx0, xshift;
01444
01445 x[0] = x[1] = 0;
01446 havedig = xshift = 0;
01447 udx0 = 1;
01448 s = *sp;
01449
01450 while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
01451 ++s;
01452 if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
01453 s += 2;
01454 while((c = *(CONST unsigned char*)++s)) {
01455 if (c >= '0' && c <= '9')
01456 c -= '0';
01457 else if (c >= 'a' && c <= 'f')
01458 c += 10 - 'a';
01459 else if (c >= 'A' && c <= 'F')
01460 c += 10 - 'A';
01461 else if (c <= ' ') {
01462 if (udx0 && havedig) {
01463 udx0 = 0;
01464 xshift = 1;
01465 }
01466 continue;
01467 }
01468 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
01469 else if ( c == ')' && havedig) {
01470 *sp = s + 1;
01471 break;
01472 }
01473 else
01474 return;
01475 #else
01476 else {
01477 do {
01478 if ( c == ')') {
01479 *sp = s + 1;
01480 break;
01481 }
01482 } while((c = *++s));
01483 break;
01484 }
01485 #endif
01486 havedig = 1;
01487 if (xshift) {
01488 xshift = 0;
01489 x[0] = x[1];
01490 x[1] = 0;
01491 }
01492 if (udx0)
01493 x[0] = (x[0] << 4) | (x[1] >> 28);
01494 x[1] = (x[1] << 4) | c;
01495 }
01496 if ((x[0] &= 0xfffff) || x[1]) {
01497 word0(*rvp) = Exp_mask | x[0];
01498 word1(*rvp) = x[1];
01499 }
01500 }
01501 #endif
01502 #endif
01503
01504 double
01505 sb_strtod
01506 #ifdef KR_headers
01507 (s00, se) CONST char *s00; char **se;
01508 #else
01509 (CONST char *s00, char **se)
01510 #endif
01511 {
01512 #ifdef Avoid_Underflow
01513 int scale;
01514 #endif
01515 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
01516 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
01517 CONST char *s, *s0, *s1;
01518 double aadj, adj;
01519 U rv, rv0, aadj1;
01520 Long L;
01521 ULong y, z;
01522 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
01523 #ifdef SET_INEXACT
01524 int inexact, oldinexact;
01525 #endif
01526 #ifdef Honor_FLT_ROUNDS
01527 int rounding;
01528 #endif
01529 #ifdef USE_LOCALE
01530 CONST char *s2;
01531 #endif
01532
01533 sign = nz0 = nz = 0;
01534 dval(rv) = 0.;
01535 for(s = s00;;s++) switch(*s) {
01536 case '-':
01537 sign = 1;
01538
01539 case '+':
01540 if (*++s)
01541 goto break2;
01542
01543 case 0:
01544 goto ret0;
01545 case '\t':
01546 case '\n':
01547 case '\v':
01548 case '\f':
01549 case '\r':
01550 case ' ':
01551 continue;
01552 default:
01553 goto break2;
01554 }
01555 break2:
01556 if (*s == '0') {
01557 nz0 = 1;
01558 while(*++s == '0') ;
01559 if (!*s)
01560 goto ret;
01561 }
01562 s0 = s;
01563 y = z = 0;
01564 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
01565 if (nd < 9)
01566 y = 10*y + c - '0';
01567 else if (nd < 16)
01568 z = 10*z + c - '0';
01569 nd0 = nd;
01570 #ifdef USE_LOCALE
01571 s1 = localeconv()->decimal_point;
01572 if (c == *s1) {
01573 c = '.';
01574 if (*++s1) {
01575 s2 = s;
01576 for(;;) {
01577 if (*++s2 != *s1) {
01578 c = 0;
01579 break;
01580 }
01581 if (!*++s1) {
01582 s = s2;
01583 break;
01584 }
01585 }
01586 }
01587 }
01588 #endif
01589 if (c == '.') {
01590 c = *++s;
01591 if (!nd) {
01592 for(; c == '0'; c = *++s)
01593 nz++;
01594 if (c > '0' && c <= '9') {
01595 s0 = s;
01596 nf += nz;
01597 nz = 0;
01598 goto have_dig;
01599 }
01600 goto dig_done;
01601 }
01602 for(; c >= '0' && c <= '9'; c = *++s) {
01603 have_dig:
01604 nz++;
01605 if (c -= '0') {
01606 nf += nz;
01607 for(i = 1; i < nz; i++)
01608 if (nd++ < 9)
01609 y *= 10;
01610 else if (nd <= DBL_DIG + 1)
01611 z *= 10;
01612 if (nd++ < 9)
01613 y = 10*y + c;
01614 else if (nd <= DBL_DIG + 1)
01615 z = 10*z + c;
01616 nz = 0;
01617 }
01618 }
01619 }
01620 dig_done:
01621 e = 0;
01622 if (c == 'e' || c == 'E') {
01623 if (!nd && !nz && !nz0) {
01624 goto ret0;
01625 }
01626 s00 = s;
01627 esign = 0;
01628 switch(c = *++s) {
01629 case '-':
01630 esign = 1;
01631 case '+':
01632 c = *++s;
01633 }
01634 if (c >= '0' && c <= '9') {
01635 while(c == '0')
01636 c = *++s;
01637 if (c > '0' && c <= '9') {
01638 L = c - '0';
01639 s1 = s;
01640 while((c = *++s) >= '0' && c <= '9')
01641 L = 10*L + c - '0';
01642 if (s - s1 > 8 || L > 19999)
01643
01644
01645
01646 e = 19999;
01647 else
01648 e = (int)L;
01649 if (esign)
01650 e = -e;
01651 }
01652 else
01653 e = 0;
01654 }
01655 else
01656 s = s00;
01657 }
01658 if (!nd) {
01659 if (!nz && !nz0) {
01660 #ifdef INFNAN_CHECK
01661
01662 switch(c) {
01663 case 'i':
01664 case 'I':
01665 if (match(&s,"nf")) {
01666 --s;
01667 if (!match(&s,"inity"))
01668 ++s;
01669 word0(rv) = 0x7ff00000;
01670 word1(rv) = 0;
01671 goto ret;
01672 }
01673 break;
01674 case 'n':
01675 case 'N':
01676 if (match(&s, "an")) {
01677 word0(rv) = NAN_WORD0;
01678 word1(rv) = NAN_WORD1;
01679 #ifndef No_Hex_NaN
01680 if (*s == '(')
01681 hexnan(&rv, &s);
01682 #endif
01683 goto ret;
01684 }
01685 }
01686 #endif
01687 ret0:
01688 s = s00;
01689 sign = 0;
01690 }
01691 goto ret;
01692 }
01693 e1 = e -= nf;
01694
01695
01696
01697
01698
01699
01700 if (!nd0)
01701 nd0 = nd;
01702 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
01703 dval(rv) = y;
01704 if (k > 9) {
01705 #ifdef SET_INEXACT
01706 if (k > DBL_DIG)
01707 oldinexact = get_inexact();
01708 #endif
01709 dval(rv) = tens[k - 9] * dval(rv) + z;
01710 }
01711 bd0 = 0;
01712 if (nd <= DBL_DIG
01713 #ifndef RND_PRODQUOT
01714 #ifndef Honor_FLT_ROUNDS
01715 && Flt_Rounds == 1
01716 #endif
01717 #endif
01718 ) {
01719 if (!e)
01720 goto ret;
01721 if (e > 0) {
01722 if (e <= Ten_pmax) {
01723 #ifdef VAX
01724 goto vax_ovfl_check;
01725 #else
01726 #ifdef Honor_FLT_ROUNDS
01727
01728 if (sign) {
01729 rv = -rv;
01730 sign = 0;
01731 }
01732 #endif
01733 rounded_product(dval(rv), tens[e]);
01734 goto ret;
01735 #endif
01736 }
01737 i = DBL_DIG - nd;
01738 if (e <= Ten_pmax + i) {
01739
01740
01741
01742 #ifdef Honor_FLT_ROUNDS
01743
01744 if (sign) {
01745 rv = -rv;
01746 sign = 0;
01747 }
01748 #endif
01749 e -= i;
01750 dval(rv) *= tens[i];
01751 #ifdef VAX
01752
01753
01754
01755 vax_ovfl_check:
01756 word0(rv) -= P*Exp_msk1;
01757 rounded_product(dval(rv), tens[e]);
01758 if ((word0(rv) & Exp_mask)
01759 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
01760 goto ovfl;
01761 word0(rv) += P*Exp_msk1;
01762 #else
01763 rounded_product(dval(rv), tens[e]);
01764 #endif
01765 goto ret;
01766 }
01767 }
01768 #ifndef Inaccurate_Divide
01769 else if (e >= -Ten_pmax) {
01770 #ifdef Honor_FLT_ROUNDS
01771
01772 if (sign) {
01773 rv = -rv;
01774 sign = 0;
01775 }
01776 #endif
01777 rounded_quotient(dval(rv), tens[-e]);
01778 goto ret;
01779 }
01780 #endif
01781 }
01782 e1 += nd - k;
01783
01784 #ifdef IEEE_Arith
01785 #ifdef SET_INEXACT
01786 inexact = 1;
01787 if (k <= DBL_DIG)
01788 oldinexact = get_inexact();
01789 #endif
01790 #ifdef Avoid_Underflow
01791 scale = 0;
01792 #endif
01793 #ifdef Honor_FLT_ROUNDS
01794 if ((rounding = Flt_Rounds) >= 2) {
01795 if (sign)
01796 rounding = rounding == 2 ? 0 : 2;
01797 else
01798 if (rounding != 2)
01799 rounding = 0;
01800 }
01801 #endif
01802 #endif
01803
01804
01805
01806 if (e1 > 0) {
01807 if ((i = e1 & 15))
01808 dval(rv) *= tens[i];
01809 if (e1 &= ~15) {
01810 if (e1 > DBL_MAX_10_EXP) {
01811 ovfl:
01812 #ifndef NO_ERRNO
01813 errno = ERANGE;
01814 #endif
01815
01816 #ifdef IEEE_Arith
01817 #ifdef Honor_FLT_ROUNDS
01818 switch(rounding) {
01819 case 0:
01820 case 3:
01821 word0(rv) = Big0;
01822 word1(rv) = Big1;
01823 break;
01824 default:
01825 word0(rv) = Exp_mask;
01826 word1(rv) = 0;
01827 }
01828 #else
01829 word0(rv) = Exp_mask;
01830 word1(rv) = 0;
01831 #endif
01832 #ifdef SET_INEXACT
01833
01834 dval(rv0) = 1e300;
01835 dval(rv0) *= dval(rv0);
01836 #endif
01837 #else
01838 word0(rv) = Big0;
01839 word1(rv) = Big1;
01840 #endif
01841 if (bd0)
01842 goto retfree;
01843 goto ret;
01844 }
01845 e1 >>= 4;
01846 for(j = 0; e1 > 1; j++, e1 >>= 1)
01847 if (e1 & 1)
01848 dval(rv) *= bigtens[j];
01849
01850 word0(rv) -= P*Exp_msk1;
01851 dval(rv) *= bigtens[j];
01852 if ((z = word0(rv) & Exp_mask)
01853 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
01854 goto ovfl;
01855 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
01856
01857
01858 word0(rv) = Big0;
01859 word1(rv) = Big1;
01860 }
01861 else
01862 word0(rv) += P*Exp_msk1;
01863 }
01864 }
01865 else if (e1 < 0) {
01866 e1 = -e1;
01867 if ((i = e1 & 15))
01868 dval(rv) /= tens[i];
01869 if (e1 >>= 4) {
01870 if (e1 >= 1 << n_bigtens)
01871 goto undfl;
01872 #ifdef Avoid_Underflow
01873 if (e1 & Scale_Bit)
01874 scale = 2*P;
01875 for(j = 0; e1 > 0; j++, e1 >>= 1)
01876 if (e1 & 1)
01877 dval(rv) *= tinytens[j];
01878 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
01879 >> Exp_shift)) > 0) {
01880
01881 if (j >= 32) {
01882 word1(rv) = 0;
01883 if (j >= 53)
01884 word0(rv) = (P+2)*Exp_msk1;
01885 else
01886 word0(rv) &= 0xffffffff << (j-32);
01887 }
01888 else
01889 word1(rv) &= 0xffffffff << j;
01890 }
01891 #else
01892 for(j = 0; e1 > 1; j++, e1 >>= 1)
01893 if (e1 & 1)
01894 dval(rv) *= tinytens[j];
01895
01896 dval(rv0) = dval(rv);
01897 dval(rv) *= tinytens[j];
01898 if (!dval(rv)) {
01899 dval(rv) = 2.*dval(rv0);
01900 dval(rv) *= tinytens[j];
01901 #endif
01902 if (!dval(rv)) {
01903 undfl:
01904 dval(rv) = 0.;
01905 #ifndef NO_ERRNO
01906 errno = ERANGE;
01907 #endif
01908 if (bd0)
01909 goto retfree;
01910 goto ret;
01911 }
01912 #ifndef Avoid_Underflow
01913 word0(rv) = Tiny0;
01914 word1(rv) = Tiny1;
01915
01916
01917
01918 }
01919 #endif
01920 }
01921 }
01922
01923
01924
01925
01926
01927 bd0 = s2b(s0, nd0, nd, y);
01928
01929 for(;;) {
01930 bd = Balloc(bd0->k);
01931 Bcopy(bd, bd0);
01932 bb = d2b(dval(rv), &bbe, &bbbits);
01933 bs = i2b(1);
01934
01935 if (e >= 0) {
01936 bb2 = bb5 = 0;
01937 bd2 = bd5 = e;
01938 }
01939 else {
01940 bb2 = bb5 = -e;
01941 bd2 = bd5 = 0;
01942 }
01943 if (bbe >= 0)
01944 bb2 += bbe;
01945 else
01946 bd2 -= bbe;
01947 bs2 = bb2;
01948 #ifdef Honor_FLT_ROUNDS
01949 if (rounding != 1)
01950 bs2++;
01951 #endif
01952 #ifdef Avoid_Underflow
01953 j = bbe - scale;
01954 i = j + bbbits - 1;
01955 if (i < Emin)
01956 j += P - Emin;
01957 else
01958 j = P + 1 - bbbits;
01959 #else
01960 #ifdef Sudden_Underflow
01961 #ifdef IBM
01962 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
01963 #else
01964 j = P + 1 - bbbits;
01965 #endif
01966 #else
01967 j = bbe;
01968 i = j + bbbits - 1;
01969 if (i < Emin)
01970 j += P - Emin;
01971 else
01972 j = P + 1 - bbbits;
01973 #endif
01974 #endif
01975 bb2 += j;
01976 bd2 += j;
01977 #ifdef Avoid_Underflow
01978 bd2 += scale;
01979 #endif
01980 i = bb2 < bd2 ? bb2 : bd2;
01981 if (i > bs2)
01982 i = bs2;
01983 if (i > 0) {
01984 bb2 -= i;
01985 bd2 -= i;
01986 bs2 -= i;
01987 }
01988 if (bb5 > 0) {
01989 bs = pow5mult(bs, bb5);
01990 bb1 = mult(bs, bb);
01991 Bfree(bb);
01992 bb = bb1;
01993 }
01994 if (bb2 > 0)
01995 bb = lshift(bb, bb2);
01996 if (bd5 > 0)
01997 bd = pow5mult(bd, bd5);
01998 if (bd2 > 0)
01999 bd = lshift(bd, bd2);
02000 if (bs2 > 0)
02001 bs = lshift(bs, bs2);
02002 delta = diff(bb, bd);
02003 dsign = delta->sign;
02004 delta->sign = 0;
02005 i = cmp(delta, bs);
02006 #ifdef Honor_FLT_ROUNDS
02007 if (rounding != 1) {
02008 if (i < 0) {
02009
02010 if (!delta->x[0] && delta->wds <= 1) {
02011
02012 #ifdef SET_INEXACT
02013 inexact = 0;
02014 #endif
02015 break;
02016 }
02017 if (rounding) {
02018 if (dsign) {
02019 adj = 1.;
02020 goto apply_adj;
02021 }
02022 }
02023 else if (!dsign) {
02024 adj = -1.;
02025 if (!word1(rv)
02026 && !(word0(rv) & Frac_mask)) {
02027 y = word0(rv) & Exp_mask;
02028 #ifdef Avoid_Underflow
02029 if (!scale || y > 2*P*Exp_msk1)
02030 #else
02031 if (y)
02032 #endif
02033 {
02034 delta = lshift(delta,Log2P);
02035 if (cmp(delta, bs) <= 0)
02036 adj = -0.5;
02037 }
02038 }
02039 apply_adj:
02040 #ifdef Avoid_Underflow
02041 if (scale && (y = word0(rv) & Exp_mask)
02042 <= 2*P*Exp_msk1)
02043 word0(adj) += (2*P+1)*Exp_msk1 - y;
02044 #else
02045 #ifdef Sudden_Underflow
02046 if ((word0(rv) & Exp_mask) <=
02047 P*Exp_msk1) {
02048 word0(rv) += P*Exp_msk1;
02049 dval(rv) += adj*ulp(dval(rv));
02050 word0(rv) -= P*Exp_msk1;
02051 }
02052 else
02053 #endif
02054 #endif
02055 dval(rv) += adj*ulp(dval(rv));
02056 }
02057 break;
02058 }
02059 adj = ratio(delta, bs);
02060 if (adj < 1.)
02061 adj = 1.;
02062 if (adj <= 0x7ffffffe) {
02063
02064 y = adj;
02065 if (y != adj) {
02066 if (!((rounding>>1) ^ dsign))
02067 y++;
02068 adj = y;
02069 }
02070 }
02071 #ifdef Avoid_Underflow
02072 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02073 word0(adj) += (2*P+1)*Exp_msk1 - y;
02074 #else
02075 #ifdef Sudden_Underflow
02076 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02077 word0(rv) += P*Exp_msk1;
02078 adj *= ulp(dval(rv));
02079 if (dsign)
02080 dval(rv) += adj;
02081 else
02082 dval(rv) -= adj;
02083 word0(rv) -= P*Exp_msk1;
02084 goto cont;
02085 }
02086 #endif
02087 #endif
02088 adj *= ulp(dval(rv));
02089 if (dsign)
02090 dval(rv) += adj;
02091 else
02092 dval(rv) -= adj;
02093 goto cont;
02094 }
02095 #endif
02096
02097 if (i < 0) {
02098
02099
02100
02101 if (dsign || word1(rv) || word0(rv) & Bndry_mask
02102 #ifdef IEEE_Arith
02103 #ifdef Avoid_Underflow
02104 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
02105 #else
02106 || (word0(rv) & Exp_mask) <= Exp_msk1
02107 #endif
02108 #endif
02109 ) {
02110 #ifdef SET_INEXACT
02111 if (!delta->x[0] && delta->wds <= 1)
02112 inexact = 0;
02113 #endif
02114 break;
02115 }
02116 if (!delta->x[0] && delta->wds <= 1) {
02117
02118 #ifdef SET_INEXACT
02119 inexact = 0;
02120 #endif
02121 break;
02122 }
02123 delta = lshift(delta,Log2P);
02124 if (cmp(delta, bs) > 0)
02125 goto drop_down;
02126 break;
02127 }
02128 if (i == 0) {
02129
02130 if (dsign) {
02131 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
02132 && word1(rv) == (
02133 #ifdef Avoid_Underflow
02134 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02135 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
02136 #endif
02137 0xffffffff)) {
02138
02139 word0(rv) = (word0(rv) & Exp_mask)
02140 + Exp_msk1
02141 #ifdef IBM
02142 | Exp_msk1 >> 4
02143 #endif
02144 ;
02145 word1(rv) = 0;
02146 #ifdef Avoid_Underflow
02147 dsign = 0;
02148 #endif
02149 break;
02150 }
02151 }
02152 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
02153 drop_down:
02154
02155 #ifdef Sudden_Underflow
02156 L = word0(rv) & Exp_mask;
02157 #ifdef IBM
02158 if (L < Exp_msk1)
02159 #else
02160 #ifdef Avoid_Underflow
02161 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
02162 #else
02163 if (L <= Exp_msk1)
02164 #endif
02165 #endif
02166 goto undfl;
02167 L -= Exp_msk1;
02168 #else
02169 #ifdef Avoid_Underflow
02170 if (scale) {
02171 L = word0(rv) & Exp_mask;
02172 if (L <= (2*P+1)*Exp_msk1) {
02173 if (L > (P+2)*Exp_msk1)
02174
02175
02176 break;
02177
02178 goto undfl;
02179 }
02180 }
02181 #endif
02182 L = (word0(rv) & Exp_mask) - Exp_msk1;
02183 #endif
02184 word0(rv) = L | Bndry_mask1;
02185 word1(rv) = 0xffffffff;
02186 #ifdef IBM
02187 goto cont;
02188 #else
02189 break;
02190 #endif
02191 }
02192 #ifndef ROUND_BIASED
02193 if (!(word1(rv) & LSB))
02194 break;
02195 #endif
02196 if (dsign)
02197 dval(rv) += ulp(dval(rv));
02198 #ifndef ROUND_BIASED
02199 else {
02200 dval(rv) -= ulp(dval(rv));
02201 #ifndef Sudden_Underflow
02202 if (!dval(rv))
02203 goto undfl;
02204 #endif
02205 }
02206 #ifdef Avoid_Underflow
02207 dsign = 1 - dsign;
02208 #endif
02209 #endif
02210 break;
02211 }
02212 if ((aadj = ratio(delta, bs)) <= 2.) {
02213 if (dsign)
02214 aadj = dval(aadj1) = 1.;
02215 else if (word1(rv) || word0(rv) & Bndry_mask) {
02216 #ifndef Sudden_Underflow
02217 if (word1(rv) == Tiny1 && !word0(rv))
02218 goto undfl;
02219 #endif
02220 aadj = 1.;
02221 dval(aadj1) = -1.;
02222 }
02223 else {
02224
02225
02226
02227 if (aadj < 2./FLT_RADIX)
02228 aadj = 1./FLT_RADIX;
02229 else
02230 aadj *= 0.5;
02231 dval(aadj1) = -aadj;
02232 }
02233 }
02234 else {
02235 aadj *= 0.5;
02236 dval(aadj1) = dsign ? aadj : -aadj;
02237 #ifdef Check_FLT_ROUNDS
02238 switch(Rounding) {
02239 case 2:
02240 dval(aadj1) -= 0.5;
02241 break;
02242 case 0:
02243 case 3:
02244 dval(aadj1) += 0.5;
02245 }
02246 #else
02247 if (Flt_Rounds == 0)
02248 dval(aadj1) += 0.5;
02249 #endif
02250 }
02251 y = word0(rv) & Exp_mask;
02252
02253
02254
02255 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
02256 dval(rv0) = dval(rv);
02257 word0(rv) -= P*Exp_msk1;
02258 adj = dval(aadj1) * ulp(dval(rv));
02259 dval(rv) += adj;
02260 if ((word0(rv) & Exp_mask) >=
02261 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
02262 if (word0(rv0) == Big0 && word1(rv0) == Big1)
02263 goto ovfl;
02264 word0(rv) = Big0;
02265 word1(rv) = Big1;
02266 goto cont;
02267 }
02268 else
02269 word0(rv) += P*Exp_msk1;
02270 }
02271 else {
02272 #ifdef Avoid_Underflow
02273 if (scale && y <= 2*P*Exp_msk1) {
02274 if (aadj <= 0x7fffffff) {
02275 if ((z = (uint32)aadj) <= 0)
02276 z = 1;
02277 aadj = z;
02278 dval(aadj1) = dsign ? aadj : -aadj;
02279 }
02280 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
02281 }
02282 adj = dval(aadj1) * ulp(dval(rv));
02283 dval(rv) += adj;
02284 #else
02285 #ifdef Sudden_Underflow
02286 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02287 dval(rv0) = dval(rv);
02288 word0(rv) += P*Exp_msk1;
02289 adj = aadj1 * ulp(dval(rv));
02290 dval(rv) += adj;
02291 #ifdef IBM
02292 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
02293 #else
02294 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
02295 #endif
02296 {
02297 if (word0(rv0) == Tiny0
02298 && word1(rv0) == Tiny1)
02299 goto undfl;
02300 word0(rv) = Tiny0;
02301 word1(rv) = Tiny1;
02302 goto cont;
02303 }
02304 else
02305 word0(rv) -= P*Exp_msk1;
02306 }
02307 else {
02308 adj = aadj1 * ulp(dval(rv));
02309 dval(rv) += adj;
02310 }
02311 #else
02312
02313
02314
02315
02316
02317
02318
02319 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
02320 aadj1 = (double)(int)(aadj + 0.5);
02321 if (!dsign)
02322 aadj1 = -aadj1;
02323 }
02324 adj = aadj1 * ulp(dval(rv));
02325 dval(rv) += adj;
02326 #endif
02327 #endif
02328 }
02329 z = word0(rv) & Exp_mask;
02330 #ifndef SET_INEXACT
02331 #ifdef Avoid_Underflow
02332 if (!scale)
02333 #endif
02334 if (y == z) {
02335
02336 L = (Long)aadj;
02337 aadj -= L;
02338
02339 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
02340 if (aadj < .4999999 || aadj > .5000001)
02341 break;
02342 }
02343 else if (aadj < .4999999/FLT_RADIX)
02344 break;
02345 }
02346 #endif
02347 cont:
02348 Bfree(bb);
02349 Bfree(bd);
02350 Bfree(bs);
02351 Bfree(delta);
02352 }
02353 #ifdef SET_INEXACT
02354 if (inexact) {
02355 if (!oldinexact) {
02356 word0(rv0) = Exp_1 + (70 << Exp_shift);
02357 word1(rv0) = 0;
02358 dval(rv0) += 1.;
02359 }
02360 }
02361 else if (!oldinexact)
02362 clear_inexact();
02363 #endif
02364 #ifdef Avoid_Underflow
02365 if (scale) {
02366 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
02367 word1(rv0) = 0;
02368 dval(rv) *= dval(rv0);
02369 #ifndef NO_ERRNO
02370
02371 if (word0(rv) == 0 && word1(rv) == 0)
02372 errno = ERANGE;
02373 #endif
02374 }
02375 #endif
02376 #ifdef SET_INEXACT
02377 if (inexact && !(word0(rv) & Exp_mask)) {
02378
02379 dval(rv0) = 1e-300;
02380 dval(rv0) *= dval(rv0);
02381 }
02382 #endif
02383 retfree:
02384 Bfree(bb);
02385 Bfree(bd);
02386 Bfree(bs);
02387 Bfree(bd0);
02388 Bfree(delta);
02389 ret:
02390 if (se)
02391 *se = (char *)s;
02392 return sign ? -dval(rv) : dval(rv);
02393 }