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