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