2009-06-28 Ozkan Sezer <sezero@users.sourceforge.net>

* mingwex/gdtoa/README.mingw mingwex/gdtoa/gdtoa_fltrnds.h: New files.
        * mingwex/gdtoa/README mingwex/gdtoa/dmisc.c mingwex/gdtoa/dtoa.c
        mingwex/gdtoa/g__fmt.c mingwex/gdtoa/g_dfmt.c mingwex/gdtoa/g_ffmt.c
        mingwex/gdtoa/g_xfmt.c mingwex/gdtoa/gd_arith.h mingwex/gdtoa/gd_qnan.h
        mingwex/gdtoa/gdtoa.c mingwex/gdtoa/gdtoa.h mingwex/gdtoa/gdtoaimp.h
        mingwex/gdtoa/gethex.c mingwex/gdtoa/gmisc.c mingwex/gdtoa/hd_init.c
        mingwex/gdtoa/hexnan.c mingwex/gdtoa/misc.c mingwex/gdtoa/qnan.c
        mingwex/gdtoa/smisc.c mingwex/gdtoa/strtodg.c mingwex/gdtoa/strtodnrp.c
        mingwex/gdtoa/strtof.c mingwex/gdtoa/strtopx.c mingwex/gdtoa/sum.c
        mingwex/gdtoa/ulp.c:  Update the gdtoa library to match the netlib.org
        sources as of Apr. 20, 2009.  Update further to match the sources in
        the mingw-w64 tree as of June 28, 2009, by removing IBM, CRAY and VAX
        code, removing KR_headers, ANSI, Void and Char ifdefs, renaming the
        double/ulong union from U to dbl_union for better grepping and white-
        space tidy-ups.
This commit is contained in:
Chris Sutcliffe 2009-07-12 22:44:37 +00:00
parent 62fb43a722
commit 15e114f5c5
28 changed files with 1345 additions and 1437 deletions

View File

@ -1,3 +1,21 @@
2009-06-28 Ozkan Sezer <sezero@users.sourceforge.net>
* mingwex/gdtoa/README.mingw mingwex/gdtoa/gdtoa_fltrnds.h: New files.
* mingwex/gdtoa/README mingwex/gdtoa/dmisc.c mingwex/gdtoa/dtoa.c
mingwex/gdtoa/g__fmt.c mingwex/gdtoa/g_dfmt.c mingwex/gdtoa/g_ffmt.c
mingwex/gdtoa/g_xfmt.c mingwex/gdtoa/gd_arith.h mingwex/gdtoa/gd_qnan.h
mingwex/gdtoa/gdtoa.c mingwex/gdtoa/gdtoa.h mingwex/gdtoa/gdtoaimp.h
mingwex/gdtoa/gethex.c mingwex/gdtoa/gmisc.c mingwex/gdtoa/hd_init.c
mingwex/gdtoa/hexnan.c mingwex/gdtoa/misc.c mingwex/gdtoa/qnan.c
mingwex/gdtoa/smisc.c mingwex/gdtoa/strtodg.c mingwex/gdtoa/strtodnrp.c
mingwex/gdtoa/strtof.c mingwex/gdtoa/strtopx.c mingwex/gdtoa/sum.c
mingwex/gdtoa/ulp.c: Update the gdtoa library to match the netlib.org
sources as of Apr. 20, 2009. Update further to match the sources in
the mingw-w64 tree as of June 28, 2009, by removing IBM, CRAY and VAX
code, removing KR_headers, ANSI, Void and Char ifdefs, renaming the
double/ulong union from U to dbl_union for better grepping and white-
space tidy-ups.
2009-06-16 Chris Sutcliffe <ir0nh34d@users.sourceforge.net>
* include/stdlib.h (_wtof): Define.

View File

@ -56,7 +56,9 @@ two letters:
whose sum is the desired value
For decimal -> binary conversions, there are three families of
helper routines: one for round-nearest:
helper routines: one for round-nearest (or the current rounding
mode on IEEE-arithmetic systems that provide the C99 fegetround()
function, if compiled with -DHonor_FLT_ROUNDS):
strtof
strtod
@ -191,6 +193,9 @@ in the buffer, if the buffer was long enough, or 0. Other forms of
conversion are easily done with the help of gdtoa(), such as %e or %f
style and conversions with direction of rounding specified (so that, if
desired, the decimal value is either >= or <= the binary value).
On IEEE-arithmetic systems that provide the C99 fegetround() function,
if compiled with -DHonor_FLT_ROUNDS, these routines honor the current
rounding mode.
For an example of more general conversions based on dtoa(), see
netlib's "printf.c from ampl/solvers".
@ -332,5 +337,21 @@ Compiling g__fmt.c, strtod.c, and strtodg.c with -DUSE_LOCALE causes
the decimal-point character to be taken from the current locale; otherwise
it is '.'.
Source files dtoa.c and strtod.c in this directory are derived from
netlib's "dtoa.c from fp" and are meant to function equivalently.
When compiled with Honor_FLT_ROUNDS #defined (on systems that provide
FLT_ROUNDS and fegetround() as specified in the C99 standard), they
honor the current rounding mode. Because FLT_ROUNDS is buggy on some
(Linux) systems -- not reflecting calls on fesetround(), as the C99
standard says it should -- when Honor_FLT_ROUNDS is #defined, the
current rounding mode is obtained from fegetround() rather than from
FLT_ROUNDS, unless Trust_FLT_ROUNDS is also #defined.
Compile with -DUSE_LOCALE to use the current locale; otherwise
decimal points are assumed to be '.'. With -DUSE_LOCALE, unless
you also compile with -DNO_LOCALE_CACHE, the details about the
current "decimal point" character string are cached and assumed not
to change during the program's execution.
Please send comments to David M. Gay (dmg at acm dot org, with " at "
changed at "@" and " dot " changed to ".").

View File

@ -0,0 +1,19 @@
The gdtoa code here is based on David M. Gay's original
gdtoa source at http://www.netlib.org/fp/ from April 20,
2009. The major changes between the original source and
the mingw port here include:
* IBM, CRAY and VAX code removed.
* KR_headers, ANSI, Void and Char ifdefs are removed.
* gdtoa symbols are prepended with "__".
* g_xfmt() uses __fpclassifyl() instead of interpreting
the flags bit-wise.
* lo0bits() and hi0bits() of misc.c replaced by wrappers
to gcc's __builtin_clz()
* The double/ulong union renamed from U to dbl_union
(grep'ped better..)
* A few compiler warning fixes here and there.
* A few other insignificant changes (if any..)
MinGW specific compile-time definitions are at the top of
gdtoaimp.h and gdtoa.h headers.

View File

@ -31,12 +31,11 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
char *
#ifdef KR_headers
rv_alloc(i) int i;
#else
rv_alloc(int i)
#ifndef MULTIPLE_THREADS
char *dtoa_result;
#endif
char *rv_alloc (int i)
{
int j, k, *r;
@ -52,14 +51,9 @@ rv_alloc(int i)
dtoa_result =
#endif
(char *)(r+1);
}
}
char *
#ifdef KR_headers
nrv_alloc(s, rve, n) char *s, **rve; int n;
#else
nrv_alloc(char *s, char **rve, int n)
#endif
char *nrv_alloc (char *s, char **rve, int n)
{
char *rv, *t;
@ -69,7 +63,7 @@ nrv_alloc(char *s, char **rve, int n)
if (rve)
*rve = t;
return rv;
}
}
/* freedtoa(s) must be used to free values s returned by dtoa
* when MULTIPLE_THREADS is #defined. It should be used in all cases,
@ -77,12 +71,7 @@ nrv_alloc(char *s, char **rve, int n)
* when MULTIPLE_THREADS is not defined.
*/
void
#ifdef KR_headers
__freedtoa(s) char *s;
#else
__freedtoa(char *s)
#endif
void __freedtoa (char *s)
{
Bigint *b = (Bigint *)((int *)s - 1);
b->maxwds = 1 << (b->k = *(int*)b);
@ -91,15 +80,9 @@ __freedtoa(char *s)
if (s == dtoa_result)
dtoa_result = 0;
#endif
}
}
int
quorem
#ifdef KR_headers
(b, S) Bigint *b, *S;
#else
(Bigint *b, Bigint *S)
#endif
int quorem (Bigint *b, Bigint *S)
{
int n;
ULong *bx, *bxe, q, *sx, *sxe;
@ -157,15 +140,16 @@ quorem
*bx++ = y & 0xffff;
#endif
#endif
}
while(sx <= sxe);
} while(sx <= sxe);
if (!*bxe) {
bx = b->x;
while(--bxe > bx && !*bxe)
--n;
b->wds = n;
}
}
}
if (cmp(b, S) >= 0) {
q++;
borrow = 0;
@ -198,15 +182,15 @@ quorem
*bx++ = y & 0xffff;
#endif
#endif
}
while(sx <= sxe);
} while(sx <= sxe);
bx = b->x;
bxe = bx + n;
if (!*bxe) {
while(--bxe > bx && !*bxe)
--n;
b->wds = n;
}
}
return q;
}
return q;
}

View File

@ -66,21 +66,13 @@ THIS SOFTWARE.
*/
#ifdef Honor_FLT_ROUNDS
#define Rounding rounding
#undef Check_FLT_ROUNDS
#define Check_FLT_ROUNDS
#else
#define Rounding Flt_Rounds
#endif
char *
__dtoa
#ifdef KR_headers
(d, mode, ndigits, decpt, sign, rve)
double d; int mode, ndigits, *decpt, *sign; char **rve;
#else
(double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
#endif
char *__dtoa (double d0, int mode, int ndigits, int *decpt, int *sign, char **rve)
{
/* Arguments ndigits, decpt, sign are similar to those
of ecvt and fcvt; trailing zeros are suppressed from
@ -116,7 +108,7 @@ __dtoa
to hold the suppressed trailing zeros.
*/
int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0,
int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
spec_case, try_quick;
Long L;
@ -125,88 +117,84 @@ __dtoa
ULong x;
#endif
Bigint *b, *b1, *delta, *mlo, *mhi, *S;
double d2, ds, eps;
union _dbl_union d, d2, eps;
double ds;
char *s, *s0;
#ifdef Honor_FLT_ROUNDS
int rounding;
#endif
#ifdef SET_INEXACT
int inexact, oldinexact;
#endif
#ifdef Honor_FLT_ROUNDS /*{*/
int Rounding;
#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
Rounding = Flt_Rounds;
#else /*}{*/
Rounding = 1;
switch(fegetround()) {
case FE_TOWARDZERO: Rounding = 0; break;
case FE_UPWARD: Rounding = 2; break;
case FE_DOWNWARD: Rounding = 3;
}
#endif /*}}*/
#endif /*}*/
#ifndef MULTIPLE_THREADS
if (dtoa_result) {
__freedtoa(dtoa_result);
dtoa_result = 0;
}
}
#endif
if (word0(d) & Sign_bit) {
d.d = d0;
if (word0(&d) & Sign_bit) {
/* set sign for everything, including 0's and NaNs */
*sign = 1;
word0(d) &= ~Sign_bit; /* clear sign bit */
}
word0(&d) &= ~Sign_bit; /* clear sign bit */
}
else
*sign = 0;
#if defined(IEEE_Arith) + defined(VAX)
#ifdef IEEE_Arith
if ((word0(d) & Exp_mask) == Exp_mask)
#else
if (word0(d) == 0x8000)
#endif
{
if ((word0(&d) & Exp_mask) == Exp_mask)
{
/* Infinity or NaN */
*decpt = 9999;
#ifdef IEEE_Arith
if (!word1(d) && !(word0(d) & 0xfffff))
if (!word1(&d) && !(word0(&d) & 0xfffff))
return nrv_alloc("Infinity", rve, 8);
#endif
return nrv_alloc("NaN", rve, 3);
}
#endif
#ifdef IBM
dval(d) += 0; /* normalize */
#endif
if (!dval(d)) {
}
if (!dval(&d)) {
*decpt = 1;
return nrv_alloc("0", rve, 1);
}
}
#ifdef SET_INEXACT
try_quick = oldinexact = get_inexact();
inexact = 1;
#endif
#ifdef Honor_FLT_ROUNDS
if ((rounding = Flt_Rounds) >= 2) {
if (Rounding >= 2) {
if (*sign)
rounding = rounding == 2 ? 0 : 2;
Rounding = Rounding == 2 ? 0 : 2;
else
if (rounding != 2)
rounding = 0;
}
if (Rounding != 2)
Rounding = 0;
}
#endif
b = d2b(dval(d), &be, &bbits);
b = d2b(dval(&d), &be, &bbits);
#ifdef Sudden_Underflow
i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
#else
if (( i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
#endif
dval(d2) = dval(d);
word0(d2) &= Frac_mask1;
word0(d2) |= Exp_11;
#ifdef IBM
if (( j = 11 - hi0bits(word0(d2) & Frac_mask) )!=0)
dval(d2) /= 1 << j;
if (( i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
#endif
dval(&d2) = dval(&d);
word0(&d2) &= Frac_mask1;
word0(&d2) |= Exp_11;
/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
* log10(x) = log(x) / log(10)
* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
* log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
* log10(&d) = (i-Bias)*log(2)/log(10) + log10(&d2)
*
* This suggests computing an approximation k to log10(d) by
* This suggests computing an approximation k to log10(&d) by
*
* k = (i - Bias)*0.301029995663981
* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
@ -224,54 +212,50 @@ __dtoa
*/
i -= Bias;
#ifdef IBM
i <<= 2;
i += j;
#endif
#ifndef Sudden_Underflow
denorm = 0;
}
}
else {
/* d is denormalized */
i = bbits + be + (Bias + (P-1) - 1);
x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
: word1(d) << (32 - i);
dval(d2) = x;
word0(d2) -= 31*Exp_msk1; /* adjust exponent */
x = i > 32 ? word0(&d) << (64 - i) | word1(&d) >> (i - 32)
: word1(&d) << (32 - i);
dval(&d2) = x;
word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
i -= (Bias + (P-1) - 1) + 1;
denorm = 1;
}
}
#endif
ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
k = (int)ds;
if (ds < 0. && ds != k)
k--; /* want k = floor(ds) */
k_check = 1;
if (k >= 0 && k <= Ten_pmax) {
if (dval(d) < tens[k])
if (dval(&d) < tens[k])
k--;
k_check = 0;
}
}
j = bbits - i - 1;
if (j >= 0) {
b2 = 0;
s2 = j;
}
}
else {
b2 = -j;
s2 = 0;
}
}
if (k >= 0) {
b5 = 0;
s5 = k;
s2 += k;
}
}
else {
b2 -= k;
b5 = -k;
s5 = 0;
}
}
if (mode < 0 || mode > 9)
mode = 0;
@ -286,12 +270,13 @@ __dtoa
if (mode > 5) {
mode -= 4;
try_quick = 0;
}
}
leftright = 1;
ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
/* silence erroneous "gcc -Wall" warning. */
switch(mode) {
case 0:
case 1:
ilim = ilim1 = -1;
i = 18;
ndigits = 0;
break;
@ -312,11 +297,11 @@ __dtoa
ilim1 = i - 1;
if (i <= 0)
i = 1;
}
}
s = s0 = rv_alloc(i);
#ifdef Honor_FLT_ROUNDS
if (mode > 1 && rounding != 1)
if (mode > 1 && Rounding != 1)
leftright = 0;
#endif
@ -325,7 +310,7 @@ __dtoa
/* Try to get by with floating-point arithmetic. */
i = 0;
dval(d2) = dval(d);
dval(&d2) = dval(&d);
k0 = k;
ilim0 = ilim;
ieps = 2; /* conservative */
@ -335,92 +320,92 @@ __dtoa
if (j & Bletch) {
/* prevent overflows */
j &= Bletch - 1;
dval(d) /= bigtens[n_bigtens-1];
dval(&d) /= bigtens[n_bigtens-1];
ieps++;
}
}
for(; j; j >>= 1, i++)
if (j & 1) {
ieps++;
ds *= bigtens[i];
}
dval(d) /= ds;
}
}
dval(&d) /= ds;
}
else if (( j1 = -k )!=0) {
dval(d) *= tens[j1 & 0xf];
dval(&d) *= tens[j1 & 0xf];
for(j = j1 >> 4; j; j >>= 1, i++)
if (j & 1) {
ieps++;
dval(d) *= bigtens[i];
}
}
if (k_check && dval(d) < 1. && ilim > 0) {
dval(&d) *= bigtens[i];
}
}
if (k_check && dval(&d) < 1. && ilim > 0) {
if (ilim1 <= 0)
goto fast_failed;
ilim = ilim1;
k--;
dval(d) *= 10.;
dval(&d) *= 10.;
ieps++;
}
dval(eps) = ieps*dval(d) + 7.;
word0(eps) -= (P-1)*Exp_msk1;
}
dval(&eps) = ieps*dval(&d) + 7.;
word0(&eps) -= (P-1)*Exp_msk1;
if (ilim == 0) {
S = mhi = 0;
dval(d) -= 5.;
if (dval(d) > dval(eps))
dval(&d) -= 5.;
if (dval(&d) > dval(&eps))
goto one_digit;
if (dval(d) < -dval(eps))
if (dval(&d) < -dval(&eps))
goto no_digits;
goto fast_failed;
}
}
#ifndef No_leftright
if (leftright) {
/* Use Steele & White method of only
* generating digits needed.
*/
dval(eps) = 0.5/tens[ilim-1] - dval(eps);
dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
for(i = 0;;) {
L = dval(d);
dval(d) -= L;
L = dval(&d);
dval(&d) -= L;
*s++ = '0' + (int)L;
if (dval(d) < dval(eps))
if (dval(&d) < dval(&eps))
goto ret1;
if (1. - dval(d) < dval(eps))
if (1. - dval(&d) < dval(&eps))
goto bump_up;
if (++i >= ilim)
break;
dval(eps) *= 10.;
dval(d) *= 10.;
}
dval(&eps) *= 10.;
dval(&d) *= 10.;
}
}
else {
#endif
/* Generate ilim digits, then fix them up. */
dval(eps) *= tens[ilim-1];
for(i = 1;; i++, dval(d) *= 10.) {
L = (Long)(dval(d));
if (!(dval(d) -= L))
dval(&eps) *= tens[ilim-1];
for(i = 1;; i++, dval(&d) *= 10.) {
L = (Long)(dval(&d));
if (!(dval(&d) -= L))
ilim = i;
*s++ = '0' + (int)L;
if (i == ilim) {
if (dval(d) > 0.5 + dval(eps))
if (dval(&d) > 0.5 + dval(&eps))
goto bump_up;
else if (dval(d) < 0.5 - dval(eps)) {
else if (dval(&d) < 0.5 - dval(&eps)) {
while(*--s == '0');
s++;
goto ret1;
}
break;
}
break;
}
#ifndef No_leftright
}
#ifndef No_leftright
}
#endif
fast_failed:
s = s0;
dval(d) = dval(d2);
dval(&d) = dval(&d2);
k = k0;
ilim = ilim0;
}
}
/* Do we have a "small" integer? */
@ -429,51 +414,51 @@ __dtoa
ds = tens[k];
if (ndigits < 0 && ilim <= 0) {
S = mhi = 0;
if (ilim < 0 || dval(d) <= 5*ds)
if (ilim < 0 || dval(&d) <= 5*ds)
goto no_digits;
goto one_digit;
}
for(i = 1;; i++, dval(d) *= 10.) {
L = (Long)(dval(d) / ds);
dval(d) -= L*ds;
}
for(i = 1;; i++, dval(&d) *= 10.) {
L = (Long)(dval(&d) / ds);
dval(&d) -= L*ds;
#ifdef Check_FLT_ROUNDS
/* If FLT_ROUNDS == 2, L will usually be high by 1 */
if (dval(d) < 0) {
if (dval(&d) < 0) {
L--;
dval(d) += ds;
}
dval(&d) += ds;
}
#endif
*s++ = '0' + (int)L;
if (!dval(d)) {
if (!dval(&d)) {
#ifdef SET_INEXACT
inexact = 0;
#endif
break;
}
}
if (i == ilim) {
#ifdef Honor_FLT_ROUNDS
if (mode > 1)
switch(rounding) {
switch(Rounding) {
case 0: goto ret1;
case 2: goto bump_up;
}
}
#endif
dval(d) += dval(d);
if (dval(d) > ds || (dval(d) == ds && L & 1)) {
dval(&d) += dval(&d);
if (dval(&d) > ds || (dval(&d) == ds && L & 1)) {
bump_up:
while(*--s == '9')
if (s == s0) {
k++;
*s = '0';
break;
}
}
++*s++;
}
break;
}
break;
}
goto ret1;
}
goto ret1;
}
m2 = b2;
m5 = b5;
@ -483,21 +468,17 @@ __dtoa
#ifndef Sudden_Underflow
denorm ? be + (Bias + (P-1) - 1 + 1) :
#endif
#ifdef IBM
1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
#else
1 + P - bbits;
#endif
b2 += i;
s2 += i;
mhi = i2b(1);
}
}
if (m2 > 0 && s2 > 0) {
i = m2 < s2 ? m2 : s2;
b2 -= i;
m2 -= i;
s2 -= i;
}
}
if (b5 > 0) {
if (leftright) {
if (m5 > 0) {
@ -505,13 +486,13 @@ __dtoa
b1 = mult(mhi, b);
Bfree(b);
b = b1;
}
}
if (( j = b5 - m5 )!=0)
b = pow5mult(b, j);
}
}
else
b = pow5mult(b, b5);
}
}
S = i2b(1);
if (s5 > 0)
S = pow5mult(S, s5);
@ -521,20 +502,20 @@ __dtoa
spec_case = 0;
if ((mode < 2 || leftright)
#ifdef Honor_FLT_ROUNDS
&& rounding == 1
&& Rounding == 1
#endif
) {
if (!word1(d) && !(word0(d) & Bndry_mask)
if (!word1(&d) && !(word0(&d) & Bndry_mask)
#ifndef Sudden_Underflow
&& word0(d) & (Exp_mask & ~Exp_msk1)
&& word0(&d) & (Exp_mask & ~Exp_msk1)
#endif
) {
/* The special case */
b2 += Log2P;
s2 += Log2P;
spec_case = 1;
}
}
}
/* Arrange for convenient computation of quotients:
* shift left if necessary so divisor has 4 leading 0 bits.
@ -555,13 +536,13 @@ __dtoa
b2 += i;
m2 += i;
s2 += i;
}
}
else if (i < 4) {
i += 28;
b2 += i;
m2 += i;
s2 += i;
}
}
if (b2 > 0)
b = lshift(b, b2);
if (s2 > 0)
@ -573,20 +554,20 @@ __dtoa
if (leftright)
mhi = multadd(mhi, 10, 0);
ilim = ilim1;
}
}
}
if (ilim <= 0 && (mode == 3 || mode == 5)) {
if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
/* no digits, fcvt style */
no_digits:
k = -1 - ndigits;
goto ret;
}
}
one_digit:
*s++ = '1';
k++;
goto ret;
}
}
if (leftright) {
if (m2 > 0)
mhi = lshift(mhi, m2);
@ -600,7 +581,7 @@ __dtoa
mhi = Balloc(mhi->k);
Bcopy(mhi, mlo);
mhi = lshift(mhi, Log2P);
}
}
for(i = 1;;i++) {
dig = quorem(b,S) + '0';
@ -612,9 +593,9 @@ __dtoa
j1 = delta->sign ? 1 : cmp(b, delta);
Bfree(delta);
#ifndef ROUND_BIASED
if (j1 == 0 && mode != 1 && !(word1(d) & 1)
if (j1 == 0 && mode != 1 && !(word1(&d) & 1)
#ifdef Honor_FLT_ROUNDS
&& rounding >= 1
&& Rounding >= 1
#endif
) {
if (dig == '9')
@ -627,11 +608,11 @@ __dtoa
#endif
*s++ = dig;
goto ret;
}
}
#endif
if (j < 0 || (j == 0 && mode != 1
#ifndef ROUND_BIASED
&& !(word1(d) & 1)
&& !(word1(&d) & 1)
#endif
)) {
if (!b->x[0] && b->wds <= 1) {
@ -639,13 +620,13 @@ __dtoa
inexact = 0;
#endif
goto accept_dig;
}
}
#ifdef Honor_FLT_ROUNDS
if (mode > 1)
switch(rounding) {
switch(Rounding) {
case 0: goto accept_dig;
case 2: goto keep_dig;
}
}
#endif /*Honor_FLT_ROUNDS*/
if (j1 > 0) {
b = lshift(b, 1);
@ -653,24 +634,24 @@ __dtoa
if ((j1 > 0 || (j1 == 0 && dig & 1))
&& dig++ == '9')
goto round_9_up;
}
}
accept_dig:
*s++ = dig;
goto ret;
}
}
if (j1 > 0) {
#ifdef Honor_FLT_ROUNDS
if (!rounding)
if (!Rounding)
goto accept_dig;
#endif
if (dig == '9') { /* possible if i == 1 */
round_9_up:
*s++ = '9';
goto roundoff;
}
}
*s++ = dig + 1;
goto ret;
}
}
#ifdef Honor_FLT_ROUNDS
keep_dig:
#endif
@ -683,9 +664,9 @@ __dtoa
else {
mlo = multadd(mlo, 10, 0);
mhi = multadd(mhi, 10, 0);
}
}
}
}
else
for(i = 1;; i++) {
*s++ = dig = quorem(b,S) + '0';
@ -694,19 +675,19 @@ __dtoa
inexact = 0;
#endif
goto ret;
}
}
if (i >= ilim)
break;
b = multadd(b, 10, 0);
}
}
/* Round off last digit */
#ifdef Honor_FLT_ROUNDS
switch(rounding) {
switch(Rounding) {
case 0: goto trimzeros;
case 2: goto roundoff;
}
}
#endif
b = lshift(b, 1);
j = cmp(b, S);
@ -717,30 +698,32 @@ __dtoa
k++;
*s++ = '1';
goto ret;
}
}
++*s++;
}
}
else {
// trimzeros:
#ifdef Honor_FLT_ROUNDS
trimzeros:
#endif
while(*--s == '0');
s++;
}
}
ret:
Bfree(S);
if (mhi) {
if (mlo && mlo != mhi)
Bfree(mlo);
Bfree(mhi);
}
}
ret1:
#ifdef SET_INEXACT
if (inexact) {
if (!oldinexact) {
word0(d) = Exp_1 + (70 << Exp_shift);
word1(d) = 0;
dval(d) += 1.;
}
word0(&d) = Exp_1 + (70 << Exp_shift);
word1(&d) = 0;
dval(&d) += 1.;
}
}
else if (!oldinexact)
clear_inexact();
#endif
@ -750,4 +733,4 @@ __dtoa
if (rve)
*rve = s;
return s0;
}
}

View File

@ -35,65 +35,108 @@ THIS SOFTWARE.
#include "locale.h"
#endif
char *
#ifdef KR_headers
__g__fmt(b, s, se, decpt, sign) char *b; char *s; char *se; int decpt; ULong sign;
#else
__g__fmt(char *b, char *s, char *se, int decpt, ULong sign)
#endif
char *__g__fmt (char *b, char *s, char *se, int decpt, ULong sign, size_t blen)
{
int i, j, k;
char *s0 = s;
char *be, *s0;
size_t len;
#ifdef USE_LOCALE
char decimalpoint = *localeconv()->decimal_point;
#ifdef NO_LOCALE_CACHE
char *decimalpoint = localeconv()->decimal_point;
size_t dlen = strlen(decimalpoint);
#else
#define decimalpoint '.'
char *decimalpoint;
static char *decimalpoint_cache;
static size_t dlen;
if (!(s0 = decimalpoint_cache)) {
s0 = localeconv()->decimal_point;
dlen = strlen(s0);
if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
strcpy(decimalpoint_cache, s0);
s0 = decimalpoint_cache;
}
}
decimalpoint = s0;
#endif
#else
#define dlen 0
#endif
s0 = s;
len = (se-s) + dlen + 6; /* 6 = sign + e+dd + trailing null */
if (blen < len)
goto ret0;
be = b + blen - 1;
if (sign)
*b++ = '-';
if (decpt <= -4 || decpt > se - s + 5) {
*b++ = *s++;
if (*s) {
*b++ = decimalpoint;
#ifdef USE_LOCALE
while((*b = *decimalpoint++))
++b;
#else
*b++ = '.';
#endif
while((*b = *s++) !=0)
b++;
}
}
*b++ = 'e';
/* sprintf(b, "%+.2d", decpt - 1); */
if (--decpt < 0) {
*b++ = '-';
decpt = -decpt;
}
}
else
*b++ = '+';
for(j = 2, k = 10; 10*k <= decpt; j++, k *= 10){}
for(;;) {
i = decpt / k;
if (b >= be)
goto ret0;
*b++ = i + '0';
if (--j <= 0)
break;
decpt -= i*k;
decpt *= 10;
}
*b = 0;
}
*b = 0;
}
else if (decpt <= 0) {
*b++ = decimalpoint;
#ifdef USE_LOCALE
while((*b = *decimalpoint++))
++b;
#else
*b++ = '.';
#endif
if (be < b - decpt + (se - s))
goto ret0;
for(; decpt < 0; decpt++)
*b++ = '0';
while((*b = *s++) !=0)
while((*b = *s++) != 0)
b++;
}
}
else {
while((*b = *s++) !=0) {
while((*b = *s++) != 0) {
b++;
if (--decpt == 0 && *s)
*b++ = decimalpoint;
if (--decpt == 0 && *s) {
#ifdef USE_LOCALE
while((*b = *decimalpoint++))
++b;
#else
*b++ = '.';
#endif
}
}
if (b + decpt > be) {
ret0:
b = 0;
goto ret;
}
for(; decpt > 0; decpt--)
*b++ = '0';
*b = 0;
}
}
ret:
__freedtoa(s0);
return b;
}
}

View File

@ -31,17 +31,17 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
char*
#ifdef KR_headers
__g_dfmt(buf, d, ndig, bufsize) char *buf; double *d; int ndig; unsigned bufsize;
#else
__g_dfmt(char *buf, double *d, int ndig, unsigned bufsize)
#endif
char *__g_dfmt (char *buf, double *d, int ndig, size_t bufsize)
{
static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, 0 };
static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, 0 };
char *b, *s, *se;
ULong bits[2], *L, sign;
int decpt, ex, i, mode;
#ifdef Honor_FLT_ROUNDS
#include "gdtoa_fltrnds.h"
#else
#define fpi &fpi0
#endif
if (ndig < 0)
ndig = 0;
@ -52,14 +52,16 @@ __g_dfmt(char *buf, double *d, int ndig, unsigned bufsize)
sign = L[_0] & 0x80000000L;
if ((L[_0] & 0x7ff00000) == 0x7ff00000) {
/* Infinity or NaN */
if (bufsize < 10)
return 0;
if (L[_0] & 0xfffff || L[_1]) {
return strcp(buf, "NaN");
}
}
b = buf;
if (sign)
*b++ = '-';
return strcp(b, "Infinity");
}
}
if (L[_1] == 0 && (L[_0] ^ sign) == 0 /*d == 0.*/) {
b = buf;
#ifndef IGNORE_ZERO_SIGN
@ -69,7 +71,7 @@ __g_dfmt(char *buf, double *d, int ndig, unsigned bufsize)
*b++ = '0';
*b = 0;
return b;
}
}
bits[0] = L[_1];
bits[1] = L[_0] & 0xfffff;
if ( (ex = (L[_0] >> 20) & 0x7ff) !=0)
@ -78,12 +80,9 @@ __g_dfmt(char *buf, double *d, int ndig, unsigned bufsize)
ex = 1;
ex -= 0x3ff + 52;
mode = 2;
if (ndig <= 0) {
if (bufsize < 25)
return 0;
if (ndig <= 0)
mode = 0;
}
i = STRTOG_Normal;
s = __gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
return __g__fmt(buf, s, se, decpt, sign);
}
s = __gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
return __g__fmt(buf, s, se, decpt, sign, bufsize);
}

View File

@ -31,17 +31,17 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
char*
#ifdef KR_headers
__g_ffmt(buf, f, ndig, bufsize) char *buf; float *f; int ndig; unsigned bufsize;
#else
__g_ffmt(char *buf, float *f, int ndig, unsigned bufsize)
#endif
char *__g_ffmt (char *buf, float *f, int ndig, size_t bufsize)
{
static FPI fpi = { 24, 1-127-24+1, 254-127-24+1, 1, 0 };
static FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, 0 };
char *b, *s, *se;
ULong bits[1], *L, sign;
int decpt, ex, i, mode;
#ifdef Honor_FLT_ROUNDS
#include "gdtoa_fltrnds.h"
#else
#define fpi &fpi0
#endif
if (ndig < 0)
ndig = 0;
@ -54,12 +54,12 @@ __g_ffmt(char *buf, float *f, int ndig, unsigned bufsize)
/* Infinity or NaN */
if (L[0] & 0x7fffff) {
return strcp(buf, "NaN");
}
}
b = buf;
if (sign)
*b++ = '-';
return strcp(b, "Infinity");
}
}
if (*f == 0.) {
b = buf;
#ifndef IGNORE_ZERO_SIGN
@ -69,7 +69,7 @@ __g_ffmt(char *buf, float *f, int ndig, unsigned bufsize)
*b++ = '0';
*b = 0;
return b;
}
}
bits[0] = L[0] & 0x7fffff;
if ( (ex = (L[0] >> 23) & 0xff) !=0)
bits[0] |= 0x800000;
@ -81,8 +81,8 @@ __g_ffmt(char *buf, float *f, int ndig, unsigned bufsize)
if (bufsize < 16)
return 0;
mode = 0;
}
i = STRTOG_Normal;
s = __gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
return __g__fmt(buf, s, se, decpt, sign);
}
i = STRTOG_Normal;
s = __gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
return __g__fmt(buf, s, se, decpt, sign, bufsize);
}

View File

@ -51,21 +51,22 @@ THIS SOFTWARE.
#define _4 0
#endif
char*
#ifdef KR_headers
__g_xfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; unsigned bufsize;
#else
__g_xfmt(char *buf, void *V, int ndig, unsigned bufsize)
#endif
char *__g_xfmt (char *buf, void *V, int ndig, size_t bufsize)
{
static FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0 };
static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0 };
char *b, *s, *se;
ULong bits[2], sign;
UShort *L;
int decpt, ex, i, mode;
#if defined(__MINGW32__) || defined(__MINGW64__)
int fptype = __fpclassifyl (*(long double*) V);
#endif /* MinGW */
#ifdef Honor_FLT_ROUNDS
#include "gdtoa_fltrnds.h"
#else
#define fpi &fpi0
#endif
if (ndig < 0)
ndig = 0;
if (bufsize < ndig + 10)
@ -73,36 +74,54 @@ __g_xfmt(char *buf, void *V, int ndig, unsigned bufsize)
L = (UShort *)V;
sign = L[_0] & 0x8000;
ex = L[_0] & 0x7fff;
ex = L[_0] & 0x7fff;
bits[1] = (L[_1] << 16) | L[_2];
bits[0] = (L[_3] << 16) | L[_4];
if (fptype & FP_NAN) /* NaN or Inf */
{
if (fptype & FP_NORMAL)
{
b = buf;
*b++ = sign ? '-': '+';
strncpy (b, "Infinity", ndig ? ndig : 8);
#if defined(__MINGW32__) || defined(__MINGW64__)
if (fptype & FP_NAN) {
/* NaN or Inf */
if (fptype & FP_NORMAL) {
b = buf;
*b++ = sign ? '-': '+';
strncpy (b, "Infinity", ndig ? ndig : 8);
return (buf + strlen (buf));
}
strncpy (buf, "NaN", ndig ? ndig : 3);
return (buf + strlen (buf));
}
strncpy (buf, "NaN", ndig ? ndig : 3);
return (buf + strlen (buf));
}
else if (fptype & FP_NORMAL) /* Normal or subnormal */
{
if (fptype & FP_ZERO)
{
}
else if (fptype & FP_NORMAL) {
/* Normal or subnormal */
if (fptype & FP_ZERO) {
i = STRTOG_Denormal;
ex = 1;
}
else
i = STRTOG_Normal;
}
#else
if (ex != 0) {
if (ex == 0x7fff) {
/* Infinity or NaN */
if (bits[0] | bits[1])
b = strcp(buf, "NaN");
else {
b = buf;
if (sign)
*b++ = '-';
b = strcp(b, "Infinity");
}
return b;
}
i = STRTOG_Normal;
}
else if (bits[0] | bits[1]) {
i = STRTOG_Denormal;
ex = 1;
}
else
i = STRTOG_Normal;
}
}
#endif
else {
i = STRTOG_Zero;
/* i = STRTOG_Zero; */
b = buf;
#ifndef IGNORE_ZERO_SIGN
if (sign)
@ -111,15 +130,14 @@ __g_xfmt(char *buf, void *V, int ndig, unsigned bufsize)
*b++ = '0';
*b = 0;
return b;
}
}
ex -= 0x3fff + 63;
mode = 2;
if (ndig <= 0) {
if (bufsize < 32)
return 0;
mode = 0;
}
s = __gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
return __g__fmt(buf, s, se, decpt, sign);
}
s = __gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
return __g__fmt(buf, s, se, decpt, sign, bufsize);
}

View File

@ -1,3 +1,6 @@
#define IEEE_8087
#define Arith_Kind_ASL 1
#define Double_Align
#define IEEE_8087
#define Arith_Kind_ASL 1
#define Double_Align
#ifdef _WIN64
#define X64_bit_pointers
#endif /* w64 */

View File

@ -1,12 +1,12 @@
#define f_QNAN 0xffc00000
#define d_QNAN0 0x0
#define d_QNAN1 0xfff80000
#define ld_QNAN0 0x0
#define ld_QNAN1 0xc0000000
#define ld_QNAN2 0xffff
#define ld_QNAN3 0x0
#define ldus_QNAN0 0x0
#define ldus_QNAN1 0x0
#define ldus_QNAN2 0x0
#define ldus_QNAN3 0xc000
#define ldus_QNAN4 0xffff
#define f_QNAN 0xffc00000
#define d_QNAN0 0x0
#define d_QNAN1 0xfff80000
#define ld_QNAN0 0x0
#define ld_QNAN1 0xc0000000
#define ld_QNAN2 0xffff
#define ld_QNAN3 0x0
#define ldus_QNAN0 0x0
#define ldus_QNAN1 0x0
#define ldus_QNAN2 0x0
#define ldus_QNAN3 0xc000
#define ldus_QNAN4 0xffff

View File

@ -31,12 +31,7 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
static Bigint *
#ifdef KR_headers
bitstob(bits, nbits, bbits) ULong *bits; int nbits; int *bbits;
#else
bitstob(ULong *bits, int nbits, int *bbits)
#endif
static Bigint *bitstob (ULong *bits, int nbits, int *bbits)
{
int i, k;
Bigint *b;
@ -47,7 +42,7 @@ bitstob(ULong *bits, int nbits, int *bbits)
while(i < nbits) {
i <<= 1;
k++;
}
}
#ifndef Pack_32
if (!k)
k = 1;
@ -60,19 +55,19 @@ bitstob(ULong *bits, int nbits, int *bbits)
#ifdef Pack_16
*x++ = (*bits >> 16) & ALL_ON;
#endif
} while(++bits <= be);
} while(++bits <= be);
i = x - x0;
while(!x0[--i])
if (!i) {
b->wds = 0;
*bbits = 0;
goto ret;
}
}
b->wds = i + 1;
*bbits = i*ULbits + 32 - hi0bits(b->x[i]);
ret:
return b;
}
}
/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
*
@ -108,15 +103,8 @@ bitstob(ULong *bits, int nbits, int *bbits)
* calculation.
*/
char *
__gdtoa
#ifdef KR_headers
(fpi, be, bits, kindp, mode, ndigits, decpt, rve)
FPI *fpi; int be; ULong *bits;
int *kindp, mode, ndigits, *decpt; char **rve;
#else
(FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve)
#endif
char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
int *decpt, char **rve)
{
/* Arguments ndigits and decpt are similar to the second and third
arguments of ecvt and fcvt; trailing zeros are suppressed from
@ -152,19 +140,20 @@ __gdtoa
to hold the suppressed trailing zeros.
*/
int bbits, b2, b5, be0, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0, inex;
int bbits, b2, b5, be0, dig, i, ieps, ilim, ilim0, ilim1, inex;
int j, j1, k, k0, k_check, kind, leftright, m2, m5, nbits;
int rdir, s2, s5, spec_case, try_quick;
Long L;
Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
double d, d2, ds, eps;
double d2, ds;
char *s, *s0;
union _dbl_union d, eps;
#ifndef MULTIPLE_THREADS
if (dtoa_result) {
__freedtoa(dtoa_result);
dtoa_result = 0;
}
}
#endif
inex = 0;
kind = *kindp &= ~STRTOG_Inexact;
@ -182,36 +171,32 @@ __gdtoa
return nrv_alloc("NaN", rve, 3);
default:
return 0;
}
}
b = bitstob(bits, nbits = fpi->nbits, &bbits);
be0 = be;
if ( (i = trailz(b)) !=0) {
rshift(b, i);
be += i;
bbits -= i;
}
}
if (!b->wds) {
Bfree(b);
ret_zero:
*decpt = 1;
return nrv_alloc("0", rve, 1);
}
}
dval(d) = b2d(b, &i);
dval(&d) = b2d(b, &i);
i = be + bbits - 1;
word0(d) &= Frac_mask1;
word0(d) |= Exp_11;
#ifdef IBM
if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
dval(d) /= 1 << j;
#endif
word0(&d) &= Frac_mask1;
word0(&d) |= Exp_11;
/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
* log10(x) = log(x) / log(10)
* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
* log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
* log10(&d) = (i-Bias)*log(2)/log(10) + log10(d2)
*
* This suggests computing an approximation k to log10(d) by
* This suggests computing an approximation k to log10(&d) by
*
* k = (i - Bias)*0.301029995663981
* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
@ -227,11 +212,7 @@ __gdtoa
* (We could get a more accurate k by invoking log10,
* but this is probably not worthwhile.)
*/
#ifdef IBM
i <<= 2;
i += j;
#endif
ds = (dval(d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
/* correct assumption about exponent range */
if ((j = i) < 0)
@ -243,50 +224,44 @@ __gdtoa
if (ds < 0. && ds != k)
k--; /* want k = floor(ds) */
k_check = 1;
#ifdef IBM
j = be + bbits - 1;
if ( (j1 = j & 3) !=0)
dval(d) *= 1 << j1;
word0(d) += j << Exp_shift - 2 & Exp_mask;
#else
word0(d) += (be + bbits - 1) << Exp_shift;
#endif
word0(&d) += (be + bbits - 1) << Exp_shift;
if (k >= 0 && k <= Ten_pmax) {
if (dval(d) < tens[k])
if (dval(&d) < tens[k])
k--;
k_check = 0;
}
}
j = bbits - i - 1;
if (j >= 0) {
b2 = 0;
s2 = j;
}
}
else {
b2 = -j;
s2 = 0;
}
}
if (k >= 0) {
b5 = 0;
s5 = k;
s2 += k;
}
}
else {
b2 -= k;
b5 = -k;
s5 = 0;
}
}
if (mode < 0 || mode > 9)
mode = 0;
try_quick = 1;
if (mode > 5) {
mode -= 4;
try_quick = 0;
}
}
leftright = 1;
ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
/* silence erroneous "gcc -Wall" warning. */
switch(mode) {
case 0:
case 1:
ilim = ilim1 = -1;
i = (int)(nbits * .30103) + 3;
ndigits = 0;
break;
@ -307,7 +282,7 @@ __gdtoa
ilim1 = i - 1;
if (i <= 0)
i = 1;
}
}
s = s0 = rv_alloc(i);
if ( (rdir = fpi->rounding - 1) !=0) {
@ -315,7 +290,7 @@ __gdtoa
rdir = 2;
if (kind & STRTOG_Neg)
rdir = 3 - rdir;
}
}
/* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
@ -328,11 +303,7 @@ __gdtoa
/* Try to get by with floating-point arithmetic. */
i = 0;
d2 = dval(d);
#ifdef IBM
if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
dval(d) /= 1 << j;
#endif
d2 = dval(&d);
k0 = k;
ilim0 = ilim;
ieps = 2; /* conservative */
@ -342,99 +313,97 @@ __gdtoa
if (j & Bletch) {
/* prevent overflows */
j &= Bletch - 1;
dval(d) /= bigtens[n_bigtens-1];
dval(&d) /= bigtens[n_bigtens-1];
ieps++;
}
}
for(; j; j >>= 1, i++)
if (j & 1) {
ieps++;
ds *= bigtens[i];
}
}
}
}
else {
ds = 1.;
if ( (j1 = -k) !=0) {
dval(d) *= tens[j1 & 0xf];
dval(&d) *= tens[j1 & 0xf];
for(j = j1 >> 4; j; j >>= 1, i++)
if (j & 1) {
ieps++;
dval(d) *= bigtens[i];
}
}
dval(&d) *= bigtens[i];
}
}
if (k_check && dval(d) < 1. && ilim > 0) {
}
if (k_check && dval(&d) < 1. && ilim > 0) {
if (ilim1 <= 0)
goto fast_failed;
ilim = ilim1;
k--;
dval(d) *= 10.;
dval(&d) *= 10.;
ieps++;
}
dval(eps) = ieps*dval(d) + 7.;
word0(eps) -= (P-1)*Exp_msk1;
}
dval(&eps) = ieps*dval(&d) + 7.;
word0(&eps) -= (P-1)*Exp_msk1;
if (ilim == 0) {
S = mhi = 0;
dval(d) -= 5.;
if (dval(d) > dval(eps))
dval(&d) -= 5.;
if (dval(&d) > dval(&eps))
goto one_digit;
if (dval(d) < -dval(eps))
if (dval(&d) < -dval(&eps))
goto no_digits;
goto fast_failed;
}
}
#ifndef No_leftright
if (leftright) {
/* Use Steele & White method of only
* generating digits needed.
*/
dval(eps) = ds*0.5/tens[ilim-1] - dval(eps);
dval(&eps) = ds*0.5/tens[ilim-1] - dval(&eps);
for(i = 0;;) {
L = (Long)(dval(d)/ds);
dval(d) -= L*ds;
L = (Long)(dval(&d)/ds);
dval(&d) -= L*ds;
*s++ = '0' + (int)L;
if (dval(d) < dval(eps)) {
if (dval(d))
if (dval(&d) < dval(&eps)) {
if (dval(&d))
inex = STRTOG_Inexlo;
goto ret1;
}
if (ds - dval(d) < dval(eps))
}
if (ds - dval(&d) < dval(&eps))
goto bump_up;
if (++i >= ilim)
break;
dval(eps) *= 10.;
dval(d) *= 10.;
}
dval(&eps) *= 10.;
dval(&d) *= 10.;
}
}
else {
#endif
/* Generate ilim digits, then fix them up. */
dval(eps) *= tens[ilim-1];
for(i = 1;; i++, dval(d) *= 10.) {
if ( (L = (Long)(dval(d)/ds)) !=0)
dval(d) -= L*ds;
dval(&eps) *= tens[ilim-1];
for(i = 1;; i++, dval(&d) *= 10.) {
if ( (L = (Long)(dval(&d)/ds)) !=0)
dval(&d) -= L*ds;
*s++ = '0' + (int)L;
if (i == ilim) {
ds *= 0.5;
if (dval(d) > ds + dval(eps))
if (dval(&d) > ds + dval(&eps))
goto bump_up;
else if (dval(d) < ds - dval(eps)) {
while(*--s == '0'){}
s++;
if (dval(d))
else if (dval(&d) < ds - dval(&eps)) {
if (dval(&d))
inex = STRTOG_Inexlo;
goto ret1;
}
break;
goto clear_trailing0;
}
break;
}
#ifndef No_leftright
}
#ifndef No_leftright
}
#endif
fast_failed:
s = s0;
dval(d) = d2;
dval(&d) = d2;
k = k0;
ilim = ilim0;
}
}
/* Do we have a "small" integer? */
@ -443,22 +412,22 @@ __gdtoa
ds = tens[k];
if (ndigits < 0 && ilim <= 0) {
S = mhi = 0;
if (ilim < 0 || dval(d) <= 5*ds)
if (ilim < 0 || dval(&d) <= 5*ds)
goto no_digits;
goto one_digit;
}
for(i = 1;; i++, dval(d) *= 10.) {
L = dval(d) / ds;
dval(d) -= L*ds;
}
for(i = 1;; i++, dval(&d) *= 10.) {
L = dval(&d) / ds;
dval(&d) -= L*ds;
#ifdef Check_FLT_ROUNDS
/* If FLT_ROUNDS == 2, L will usually be high by 1 */
if (dval(d) < 0) {
if (dval(&d) < 0) {
L--;
dval(d) += ds;
}
dval(&d) += ds;
}
#endif
*s++ = '0' + (int)L;
if (dval(d) == 0.)
if (dval(&d) == 0.)
break;
if (i == ilim) {
if (rdir) {
@ -466,9 +435,9 @@ __gdtoa
goto bump_up;
inex = STRTOG_Inexlo;
goto ret1;
}
dval(d) += dval(d);
if (dval(d) > ds || (dval(d) == ds && L & 1)) {
}
dval(&d) += dval(&d);
if (dval(&d) > ds || (dval(&d) == ds && L & 1)) {
bump_up:
inex = STRTOG_Inexhi;
while(*--s == '9')
@ -476,16 +445,20 @@ __gdtoa
k++;
*s = '0';
break;
}
}
++*s++;
}
else
inex = STRTOG_Inexlo;
break;
}
else {
inex = STRTOG_Inexlo;
clear_trailing0:
while(*--s == '0'){}
++s;
}
break;
}
goto ret1;
}
goto ret1;
}
m2 = b2;
m5 = b5;
@ -496,7 +469,7 @@ __gdtoa
if (be - i++ < fpi->emin)
/* denormal */
i = be - fpi->emin + 1;
}
}
else {
j = ilim - 1;
if (m5 >= j)
@ -505,22 +478,22 @@ __gdtoa
s5 += j -= m5;
b5 += j;
m5 = 0;
}
}
if ((i = ilim) < 0) {
m2 -= i;
i = 0;
}
}
}
b2 += i;
s2 += i;
mhi = i2b(1);
}
}
if (m2 > 0 && s2 > 0) {
i = m2 < s2 ? m2 : s2;
b2 -= i;
m2 -= i;
s2 -= i;
}
}
if (b5 > 0) {
if (leftright) {
if (m5 > 0) {
@ -528,13 +501,13 @@ __gdtoa
b1 = mult(mhi, b);
Bfree(b);
b = b1;
}
}
if ( (j = b5 - m5) !=0)
b = pow5mult(b, j);
}
}
else
b = pow5mult(b, b5);
}
}
S = i2b(1);
if (s5 > 0)
S = pow5mult(S, s5);
@ -548,8 +521,8 @@ __gdtoa
b2++;
s2++;
spec_case = 1;
}
}
}
/* Arrange for convenient computation of quotients:
* shift left if necessary so divisor has 4 leading 0 bits.
@ -558,28 +531,11 @@ __gdtoa
* and for all and pass them and a shift to quorem, so it
* can do shifts and ors to compute the numerator for q.
*/
#ifdef Pack_32
if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) !=0)
i = 32 - i;
#else
if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) !=0)
i = 16 - i;
#endif
if (i > 4) {
i -= 4;
b2 += i;
m2 += i;
s2 += i;
}
else if (i < 4) {
i += 28;
b2 += i;
m2 += i;
s2 += i;
}
if (b2 > 0)
i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
m2 += i;
if ((b2 += i) > 0)
b = lshift(b, b2);
if (s2 > 0)
if ((s2 += i) > 0)
S = lshift(S, s2);
if (k_check) {
if (cmp(b,S) < 0) {
@ -588,8 +544,8 @@ __gdtoa
if (leftright)
mhi = multadd(mhi, 10, 0);
ilim = ilim1;
}
}
}
if (ilim <= 0 && mode > 2) {
if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
/* no digits, fcvt style */
@ -597,13 +553,13 @@ __gdtoa
k = -1 - ndigits;
inex = STRTOG_Inexlo;
goto ret;
}
}
one_digit:
inex = STRTOG_Inexhi;
*s++ = '1';
k++;
goto ret;
}
}
if (leftright) {
if (m2 > 0)
mhi = lshift(mhi, m2);
@ -617,7 +573,7 @@ __gdtoa
mhi = Balloc(mhi->k);
Bcopy(mhi, mlo);
mhi = lshift(mhi, 1);
}
}
for(i = 1;;i++) {
dig = quorem(b,S) + '0';
@ -635,14 +591,14 @@ __gdtoa
if (j <= 0) {
if (b->wds > 1 || b->x[0])
inex = STRTOG_Inexlo;
}
}
else {
dig++;
inex = STRTOG_Inexhi;
}
}
*s++ = dig;
goto ret;
}
}
#endif
if (j < 0 || (j == 0 && !mode
#ifndef ROUND_BIASED
@ -653,7 +609,7 @@ __gdtoa
if (rdir == 2) {
inex = STRTOG_Inexlo;
goto accept;
}
}
while (cmp(S,mhi) > 0) {
*s++ = dig;
mhi1 = multadd(mhi, 10, 0);
@ -662,12 +618,12 @@ __gdtoa
mhi = mhi1;
b = multadd(b, 10, 0);
dig = quorem(b,S) + '0';
}
}
if (dig++ == '9')
goto round_9_up;
inex = STRTOG_Inexhi;
goto accept;
}
}
if (j1 > 0) {
b = lshift(b, 1);
j1 = cmp(b, S);
@ -675,24 +631,24 @@ __gdtoa
&& dig++ == '9')
goto round_9_up;
inex = STRTOG_Inexhi;
}
}
if (b->wds > 1 || b->x[0])
inex = STRTOG_Inexlo;
accept:
*s++ = dig;
goto ret;
}
}
if (j1 > 0 && rdir != 2) {
if (dig == '9') { /* possible if i == 1 */
round_9_up:
*s++ = '9';
inex = STRTOG_Inexhi;
goto roundoff;
}
}
inex = STRTOG_Inexhi;
*s++ = dig + 1;
goto ret;
}
}
*s++ = dig;
if (i == ilim)
break;
@ -702,16 +658,16 @@ __gdtoa
else {
mlo = multadd(mlo, 10, 0);
mhi = multadd(mhi, 10, 0);
}
}
}
}
else
for(i = 1;; i++) {
*s++ = dig = quorem(b,S) + '0';
if (i >= ilim)
break;
b = multadd(b, 10, 0);
}
}
/* Round off last digit */
@ -719,7 +675,7 @@ __gdtoa
if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
goto chopzeros;
goto roundoff;
}
}
b = lshift(b, 1);
j = cmp(b, S);
if (j > 0 || (j == 0 && dig & 1)) {
@ -730,23 +686,23 @@ __gdtoa
k++;
*s++ = '1';
goto ret;
}
}
++*s++;
}
}
else {
chopzeros:
if (b->wds > 1 || b->x[0])
inex = STRTOG_Inexlo;
while(*--s == '0'){}
s++;
}
++s;
}
ret:
Bfree(S);
if (mhi) {
if (mlo && mlo != mhi)
Bfree(mlo);
Bfree(mhi);
}
}
ret1:
Bfree(b);
*s = 0;
@ -755,4 +711,4 @@ __gdtoa
*rve = s;
*kindp |= inex;
return s0;
}
}

View File

@ -36,6 +36,15 @@ THIS SOFTWARE.
#define GDTOA_H_INCLUDED
#include "gd_arith.h"
#include <stddef.h> /* for size_t */
#if defined(__MINGW32__) || defined(__MINGW64__)
/* keep the 'Long' definition as 'long' for compatibility
* with older/other software. long in w64 is 32 bits anyway..
*/
#define Long long /* int */
#undef NO_LONG_LONG /* we have long long type */
#endif /* MinGW */
#ifndef Long
#define Long long
@ -47,25 +56,7 @@ typedef unsigned Long ULong;
typedef unsigned short UShort;
#endif
#ifndef ANSI
#ifdef KR_headers
#define ANSI(x) ()
#define Void /*nothing*/
#else
#define ANSI(x) x
#define Void void
#endif
#endif /* ANSI */
#ifndef CONST
#ifdef KR_headers
#define CONST /* blank */
#else
#define CONST const
#endif
#endif /* CONST */
enum { /* return values from strtodg */
enum { /* return values from strtodg */
STRTOG_Zero = 0,
STRTOG_Normal = 1,
STRTOG_Denormal = 2,
@ -77,49 +68,49 @@ typedef unsigned short UShort;
/* The following may be or-ed into one of the above values. */
STRTOG_Neg = 0x08,
STRTOG_Inexlo = 0x10,
STRTOG_Inexhi = 0x20,
STRTOG_Neg = 0x08, /* does not affect STRTOG_Inexlo or STRTOG_Inexhi */
STRTOG_Inexlo = 0x10, /* returned result rounded toward zero */
STRTOG_Inexhi = 0x20, /* returned result rounded away from zero */
STRTOG_Inexact = 0x30,
STRTOG_Underflow= 0x40,
STRTOG_Overflow = 0x80
};
};
typedef struct
typedef struct
FPI {
int nbits;
int emin;
int emax;
int rounding;
int sudden_underflow;
} FPI;
} FPI;
enum { /* FPI.rounding values: same as FLT_ROUNDS */
FPI_Round_zero = 0,
FPI_Round_near = 1,
FPI_Round_up = 2,
FPI_Round_down = 3
};
};
#ifdef __cplusplus
extern "C" {
#endif
extern char* __dtoa ANSI((double d, int mode, int ndigits, int *decpt,
int *sign, char **rve));
extern char* __gdtoa ANSI((FPI *fpi, int be, ULong *bits, int *kindp,
int mode, int ndigits, int *decpt, char **rve));
extern void __freedtoa ANSI((char*));
extern char* __dtoa (double d, int mode, int ndigits, int *decpt,
int *sign, char **rve);
extern char* __gdtoa (FPI *fpi, int be, ULong *bits, int *kindp,
int mode, int ndigits, int *decpt, char **rve);
extern void __freedtoa (char *);
extern int __strtodg ANSI((CONST char*, char**, FPI*, Long*, ULong*));
extern float __strtof ANSI((CONST char *, char **));
extern double __strtod ANSI((CONST char *, char **));
extern long double strtold ANSI((CONST char *, char **));
extern float __strtof (const char *, char **);
extern double __strtod (const char *, char **);
extern long double __strtold (const char *, char **);
extern int __strtodg (const char *, char **, FPI *, Long *, ULong *);
extern char* __g__fmt ANSI((char *, char *, char *e, int, ULong));
extern char* __g_dfmt ANSI((char*, double*, int, unsigned));
extern char* __g_ffmt ANSI((char*, float*, int, unsigned));
extern char* __g_xfmt ANSI((char*, void*, int, unsigned));
extern char* __g__fmt (char*, char*, char*, int, ULong, size_t);
extern char* __g_dfmt (char*, double*, int, size_t);
extern char* __g_ffmt (char*, float*, int, size_t);
extern char* __g_xfmt (char*, void*, int, size_t);
#ifdef __cplusplus
}

View File

@ -0,0 +1,18 @@
FPI *fpi, fpi1;
int Rounding;
#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
Rounding = Flt_Rounds;
#else /*}{*/
Rounding = 1;
switch(fegetround()) {
case FE_TOWARDZERO: Rounding = 0; break;
case FE_UPWARD: Rounding = 2; break;
case FE_DOWNWARD: Rounding = 3;
}
#endif /*}}*/
fpi = &fpi0;
if (Rounding != 1) {
fpi1 = fpi0;
fpi = &fpi1;
fpi1.rounding = Rounding;
}

View File

@ -2,7 +2,7 @@
The author of this software is David M. Gay.
Copyright (C) 1998-2008 by Lucent Technologies
Copyright (C) 1998-2000 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
@ -126,17 +126,22 @@ THIS SOFTWARE.
* conversions of IEEE doubles in single-threaded executions with
* 8-byte pointers, PRIVATE_MEM >= 7400 appears to suffice; with
* 4-byte pointers, PRIVATE_MEM >= 7112 appears adequate.
* #define INFNAN_CHECK on IEEE systems to cause strtod to check for
* Infinity and NaN (case insensitively).
* #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
* #defined automatically on IEEE systems. On such systems,
* when INFNAN_CHECK is #defined, strtod checks
* for Infinity and NaN (case insensitively).
* When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
* strtodg also accepts (case insensitively) strings of the form
* NaN(x), where x is a string of hexadecimal digits and spaces;
* if there is only one string of hexadecimal digits, it is taken
* for the fraction bits of the resulting NaN; if there are two or
* more strings of hexadecimal digits, each string is assigned
* to the next available sequence of 32-bit words of fractions
* bits (starting with the most significant), right-aligned in
* each sequence.
* NaN(x), where x is a string of hexadecimal digits (optionally
* preceded by 0x or 0X) and spaces; if there is only one string
* of hexadecimal digits, it is taken for the fraction bits of the
* resulting NaN; if there are two or more strings of hexadecimal
* digits, each string is assigned to the next available sequence
* of 32-bit words of fractions bits (starting with the most
* significant), right-aligned in each sequence.
* Unless GDTOA_NON_PEDANTIC_NANCHECK is #defined, input "NaN(...)"
* is consumed even when ... has the wrong form (in which case the
* "(...)" is consumed but ignored).
* #define MULTIPLE_THREADS if the system offers preemptively scheduled
* multiple threads. In this case, you must provide (or suitably
* #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
@ -157,11 +162,6 @@ THIS SOFTWARE.
* #define NO_STRING_H to use private versions of memcpy.
* On some K&R systems, it may also be necessary to
* #define DECLARE_SIZE_T in this case.
* #define YES_ALIAS to permit aliasing certain double values with
* arrays of ULongs. This leads to slightly better code with
* some compilers and was always used prior to 19990916, but it
* is not strictly legal and can cause trouble with aggressively
* optimizing compilers (e.g., gcc 2.95.1 under -O2).
* #define USE_LOCALE to use the current locale's decimal_point value.
*/
@ -170,9 +170,15 @@ THIS SOFTWARE.
#include "gdtoa.h"
#include "gd_qnan.h"
#define INFNAN_CHECK 1
#if defined(__MINGW32__) || defined(__MINGW64__)
#define MULTIPLE_THREADS 1
#define USE_LOCALE 1
#define NO_LOCALE_CACHE 1
#endif /* MinGW */
#ifdef Honor_FLT_ROUNDS
#include <fenv.h>
#endif
#ifdef DEBUG
#include <stdio.h>
@ -182,14 +188,8 @@ THIS SOFTWARE.
#include <stdlib.h>
#include <string.h>
#ifdef KR_headers
#define Char char
#else
#define Char void
#endif
#ifdef MALLOC
extern Char *MALLOC ANSI((size_t));
extern void *MALLOC (size_t);
#else
#define MALLOC malloc
#endif
@ -204,6 +204,14 @@ extern Char *MALLOC ANSI((size_t));
#endif
#include <errno.h>
#ifdef NO_ERRNO
#define SET_ERRNO(x)
#else
#define SET_ERRNO(x) \
errno = (x)
#endif
#ifdef Bad_float_h
#ifdef IEEE_Arith
@ -264,27 +272,16 @@ extern "C" {
Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
#endif
typedef union { double d; ULong L[2]; } U;
typedef union _dbl_union { double d; ULong L[2]; } dbl_union;
#ifdef YES_ALIAS
#define dval(x) x
#ifdef IEEE_8087
#define word0(x) ((ULong *)&x)[1]
#define word1(x) ((ULong *)&x)[0]
#define word0(x) (x)->L[1]
#define word1(x) (x)->L[0]
#else
#define word0(x) ((ULong *)&x)[0]
#define word1(x) ((ULong *)&x)[1]
#define word0(x) (x)->L[0]
#define word1(x) (x)->L[1]
#endif
#else /* !YES_ALIAS */
#ifdef IEEE_8087
#define word0(x) ((U*)&x)->L[1]
#define word1(x) ((U*)&x)->L[0]
#else
#define word0(x) ((U*)&x)->L[0]
#define word1(x) ((U*)&x)->L[1]
#endif
#define dval(x) ((U*)&x)->d
#endif /* YES_ALIAS */
#define dval(x) (x)->d
/* The following definition of Storeinc is appropriate for MIPS processors.
* An alternative that might be better on some machines is
@ -403,11 +400,7 @@ typedef union { double d; ULong L[2]; } U;
#ifdef RND_PRODQUOT
#define rounded_product(a,b) a = rnd_prod(a, b)
#define rounded_quotient(a,b) a = rnd_quot(a, b)
#ifdef KR_headers
extern double rnd_prod(), rnd_quot();
#else
extern double rnd_prod(double, double), rnd_quot(double, double);
#endif
#else
#define rounded_product(a,b) a *= b
#define rounded_quotient(a,b) a /= b
@ -458,23 +451,22 @@ extern double rnd_prod(double, double), rnd_quot(double, double);
#define FREE_DTOA_LOCK(n) /*nothing*/
#endif
#define Kmax 15
#define Kmax 9
#define Bigint __Bigint
struct
struct
Bigint {
struct Bigint *next;
int k, maxwds, sign, wds;
ULong x[1];
};
typedef struct Bigint Bigint;
};
typedef struct Bigint Bigint;
#ifdef NO_STRING_H
#ifdef DECLARE_SIZE_T
typedef unsigned int size_t;
#endif
extern void memcpy_D2A ANSI((void*, const void*, size_t));
extern void memcpy_D2A (void*, const void*, size_t);
#define Bcopy(x,y) memcpy_D2A(&x->sign,&y->sign,y->wds*sizeof(ULong) + 2*sizeof(int))
#else /* !NO_STRING_H */
#define Bcopy(x,y) memcpy(&x->sign,&y->sign,y->wds*sizeof(ULong) + 2*sizeof(int))
@ -484,15 +476,15 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t));
static inline int
__lo0bits_D2A (ULong *y)
{
int ret = __builtin_ctz(*y);
*y = *y >> ret;
return ret;
int ret = __builtin_ctz(*y);
*y = *y >> ret;
return ret;
}
static inline int
__hi0bits_D2A (ULong y)
{
return __builtin_clz(y);
return __builtin_clz(y);
}
#endif
@ -543,48 +535,48 @@ __hi0bits_D2A (ULong y)
#define trailz __trailz_D2A
#define ulp __ulp_D2A
extern char *dtoa_result;
extern CONST double bigtens[], tens[], tinytens[];
extern unsigned char hexdig[];
extern char *dtoa_result;
extern const double bigtens[], tens[], tinytens[];
extern unsigned char hexdig[];
extern Bigint *Balloc ANSI((int));
extern void Bfree ANSI((Bigint*));
extern void ULtof ANSI((ULong*, ULong*, Long, int));
extern void ULtod ANSI((ULong*, ULong*, Long, int));
extern void ULtodd ANSI((ULong*, ULong*, Long, int));
extern void ULtoQ ANSI((ULong*, ULong*, Long, int));
extern void ULtox ANSI((UShort*, ULong*, Long, int));
extern void ULtoxL ANSI((ULong*, ULong*, Long, int));
extern ULong any_on ANSI((Bigint*, int));
extern double b2d ANSI((Bigint*, int*));
extern int cmp ANSI((Bigint*, Bigint*));
extern void copybits ANSI((ULong*, int, Bigint*));
extern Bigint *d2b ANSI((double, int*, int*));
extern int decrement ANSI((Bigint*));
extern Bigint *diff ANSI((Bigint*, Bigint*));
extern int gethex ANSI((CONST char**, FPI*, Long*, Bigint**, int));
extern void hexdig_init_D2A(Void);
extern int hexnan ANSI((CONST char**, FPI*, ULong*));
extern int hi0bits_D2A ANSI((ULong));
extern Bigint *i2b ANSI((int));
extern Bigint *increment ANSI((Bigint*));
extern int lo0bits ANSI((ULong*));
extern Bigint *lshift ANSI((Bigint*, int));
extern int match ANSI((CONST char**, char*));
extern Bigint *mult ANSI((Bigint*, Bigint*));
extern Bigint *multadd ANSI((Bigint*, int, int));
extern char *nrv_alloc ANSI((char*, char **, int));
extern Bigint *pow5mult ANSI((Bigint*, int));
extern int quorem ANSI((Bigint*, Bigint*));
extern double ratio ANSI((Bigint*, Bigint*));
extern void rshift ANSI((Bigint*, int));
extern char *rv_alloc ANSI((int));
extern Bigint *s2b ANSI((CONST char*, int, int, ULong));
extern Bigint *set_ones ANSI((Bigint*, int));
extern char *strcp ANSI((char*, const char*));
extern Bigint *sum ANSI((Bigint*, Bigint*));
extern int trailz ANSI((Bigint*));
extern double ulp ANSI((double));
extern Bigint *Balloc (int);
extern void Bfree (Bigint*);
extern void ULtof (ULong*, ULong*, Long, int);
extern void ULtod (ULong*, ULong*, Long, int);
extern void ULtodd (ULong*, ULong*, Long, int);
extern void ULtoQ (ULong*, ULong*, Long, int);
extern void ULtox (UShort*, ULong*, Long, int);
extern void ULtoxL (ULong*, ULong*, Long, int);
extern ULong any_on (Bigint*, int);
extern double b2d (Bigint*, int*);
extern int cmp (Bigint*, Bigint*);
extern void copybits (ULong*, int, Bigint*);
extern Bigint *d2b (double, int*, int*);
extern void decrement (Bigint*);
extern Bigint *diff (Bigint*, Bigint*);
extern int gethex (const char**, FPI*, Long*, Bigint**, int);
extern void hexdig_init_D2A(void);
extern int hexnan (const char**, FPI*, ULong*);
extern int hi0bits_D2A (ULong);
extern Bigint *i2b (int);
extern Bigint *increment (Bigint*);
extern int lo0bits (ULong*);
extern Bigint *lshift (Bigint*, int);
extern int match (const char**, char*);
extern Bigint *mult (Bigint*, Bigint*);
extern Bigint *multadd (Bigint*, int, int);
extern char *nrv_alloc (char*, char **, int);
extern Bigint *pow5mult (Bigint*, int);
extern int quorem (Bigint*, Bigint*);
extern double ratio (Bigint*, Bigint*);
extern void rshift (Bigint*, int);
extern char *rv_alloc (int);
extern Bigint *s2b (const char*, int, int, ULong, int);
extern Bigint *set_ones (Bigint*, int);
extern char *strcp (char*, const char*);
extern Bigint *sum (Bigint*, Bigint*);
extern int trailz (Bigint*);
extern double ulp (dbl_union *);
#ifdef __cplusplus
}
@ -599,6 +591,10 @@ __hi0bits_D2A (ULong y)
* (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
*/
#ifdef IEEE_Arith
#ifndef NO_INFNAN_CHECK
#undef INFNAN_CHECK
#define INFNAN_CHECK
#endif
#ifdef IEEE_MC68k
#define _0 0
#define _1 1

View File

@ -35,29 +35,38 @@ THIS SOFTWARE.
#include "locale.h"
#endif
int
#ifdef KR_headers
gethex(sp, fpi, exp, bp, sign)
CONST char **sp; FPI *fpi; Long *exp; Bigint **bp; int sign;
#else
gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
#endif
int gethex (const char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
{
Bigint *b;
CONST unsigned char *decpt, *s0, *s, *s1;
int esign, havedig, irv, k, n, nbits, up, zret;
const unsigned char *decpt, *s0, *s, *s1;
int big, esign, havedig, irv, j, k, n, n0, nbits, up, zret;
ULong L, lostbits, *x;
Long e, e1;
#ifdef USE_LOCALE
unsigned char decimalpoint = *localeconv()->decimal_point;
int i;
const unsigned char *decimalpoint;
#ifdef NO_LOCALE_CACHE
decimalpoint = (unsigned char *)localeconv()->decimal_point;
#else
#define decimalpoint '.'
static unsigned char *decimalpoint_cache;
if (!(s0 = decimalpoint_cache)) {
s0 = (unsigned char *)localeconv()->decimal_point;
decimalpoint_cache = (unsigned char *)
MALLOC(strlen((char *)s0) + 1);
if (decimalpoint_cache) {
strcpy((char *)decimalpoint_cache, (char *)s0);
s0 = decimalpoint_cache;
}
}
decimalpoint = s0;
#endif
#endif
if (!hexdig['0'])
hexdig_init_D2A();
*bp = 0;
havedig = 0;
s0 = *(CONST unsigned char **)sp + 2;
s0 = *(const unsigned char **)sp + 2;
while(s0[havedig] == '0')
havedig++;
s0 += havedig;
@ -65,11 +74,21 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
decpt = 0;
zret = 0;
e = 0;
if (!hexdig[*s]) {
if (hexdig[*s])
havedig++;
else {
zret = 1;
if (*s != decimalpoint)
#ifdef USE_LOCALE
for(i = 0; decimalpoint[i]; ++i) {
if (s[i] != decimalpoint[i])
goto pcheck;
}
decpt = s += i;
#else
if (*s != '.')
goto pcheck;
decpt = ++s;
#endif
if (!hexdig[*s])
goto pcheck;
while(*s == '0')
@ -78,64 +97,134 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
zret = 0;
havedig = 1;
s0 = s;
}
}
while(hexdig[*s])
s++;
if (*s == decimalpoint && !decpt) {
#ifdef USE_LOCALE
if (*s == *decimalpoint && !decpt) {
for(i = 1; decimalpoint[i]; ++i) {
if (s[i] != decimalpoint[i])
goto pcheck;
}
decpt = s += i;
#else
if (*s == '.' && !decpt) {
decpt = ++s;
#endif
while(hexdig[*s])
s++;
}
}/*}*/
if (decpt)
e = -(((Long)(s-decpt)) << 2);
pcheck:
s1 = s;
big = esign = 0;
switch(*s) {
case 'p':
case 'P':
esign = 0;
switch(*++s) {
case '-':
esign = 1;
/* no break */
case '+':
s++;
}
}
if ((n = hexdig[*s]) == 0 || n > 0x19) {
s = s1;
break;
}
}
e1 = n - 0x10;
while((n = hexdig[*++s]) !=0 && n <= 0x19)
while((n = hexdig[*++s]) !=0 && n <= 0x19) {
if (e1 & 0xf8000000)
big = 1;
e1 = 10*e1 + n - 0x10;
}
if (esign)
e1 = -e1;
e += e1;
}
}
*sp = (char*)s;
if (!havedig)
*sp = (char*)s0 - 1;
if (zret)
return havedig ? STRTOG_Zero : STRTOG_NoNumber;
return STRTOG_Zero;
if (big) {
if (esign) {
switch(fpi->rounding) {
case FPI_Round_up:
if (sign)
break;
goto ret_tiny;
case FPI_Round_down:
if (!sign)
break;
goto ret_tiny;
}
goto retz;
ret_tiny:
b = Balloc(0);
b->wds = 1;
b->x[0] = 1;
goto dret;
}
switch(fpi->rounding) {
case FPI_Round_near:
goto ovfl1;
case FPI_Round_up:
if (!sign)
goto ovfl1;
goto ret_big;
case FPI_Round_down:
if (sign)
goto ovfl1;
goto ret_big;
}
ret_big:
nbits = fpi->nbits;
n0 = n = nbits >> kshift;
if (nbits & kmask)
++n;
for(j = n, k = 0; j >>= 1; ++k);
*bp = b = Balloc(k);
b->wds = n;
for(j = 0; j < n0; ++j)
b->x[j] = ALL_ON;
if (n > n0)
b->x[j] = ULbits >> (ULbits - (nbits & kmask));
*exp = fpi->emin;
return STRTOG_Normal | STRTOG_Inexlo;
}
n = s1 - s0 - 1;
for(k = 0; n > 7; n >>= 1)
for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
k++;
b = Balloc(k);
x = b->x;
n = 0;
L = 0;
#ifdef USE_LOCALE
for(i = 0; decimalpoint[i+1]; ++i);
#endif
while(s1 > s0) {
if (*--s1 == decimalpoint)
#ifdef USE_LOCALE
if (*--s1 == decimalpoint[i]) {
s1 -= i;
continue;
if (n == 32) {
}
#else
if (*--s1 == '.')
continue;
#endif
if (n == ULbits) {
*x++ = L;
L = 0;
n = 0;
}
}
L |= (hexdig[*s1] & 0x0f) << n;
n += 4;
}
}
*x++ = L;
b->wds = n = x - b->x;
n = 32*n - hi0bits(L);
n = ULbits*n - hi0bits(L);
nbits = fpi->nbits;
lostbits = 0;
x = b->x;
@ -146,25 +235,26 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
k = n - 1;
if (x[k>>kshift] & 1 << (k & kmask)) {
lostbits = 2;
if (k > 1 && any_on(b,k-1))
if (k > 0 && any_on(b,k))
lostbits = 3;
}
}
}
rshift(b, n);
e += n;
}
}
else if (n < nbits) {
n = nbits - n;
b = lshift(b, n);
e -= n;
x = b->x;
}
}
if (e > fpi->emax) {
ovfl:
Bfree(b);
*bp = 0;
ovfl1:
SET_ERRNO(ERANGE);
return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
}
}
irv = STRTOG_Normal;
if (e < fpi->emin) {
irv = STRTOG_Denormal;
@ -182,17 +272,20 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
case FPI_Round_down:
if (sign) {
one_bit:
*exp = fpi->emin;
x[0] = b->wds = 1;
dret:
*bp = b;
*exp = fpi->emin;
SET_ERRNO(ERANGE);
return STRTOG_Denormal | STRTOG_Inexhi
| STRTOG_Underflow;
}
}
Bfree(b);
*bp = 0;
return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow;
}
}
Bfree(b);
retz:
SET_ERRNO(ERANGE);
return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow;
}
k = n - 1;
if (lostbits)
lostbits = 1;
@ -203,7 +296,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
nbits -= n;
rshift(b,n);
e = fpi->emin;
}
}
if (lostbits) {
up = 0;
switch(fpi->rounding) {
@ -211,7 +304,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
break;
case FPI_Round_near:
if (lostbits & 2
&& (lostbits & 1) | (x[0] & 1))
&& (lostbits | x[0]) & 1)
up = 1;
break;
case FPI_Round_up:
@ -219,7 +312,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
break;
case FPI_Round_down:
up = sign;
}
}
if (up) {
k = b->wds;
b = increment(b);
@ -228,20 +321,20 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
if (nbits == fpi->nbits - 1
&& x[nbits >> kshift] & 1 << (nbits & kmask))
irv = STRTOG_Normal;
}
}
else if (b->wds > k
|| ((n = nbits & kmask) !=0
&& hi0bits(x[k-1]) < 32-n)) {
&& hi0bits(x[k-1]) < 32-n)) {
rshift(b,1);
if (++e > fpi->emax)
goto ovfl;
}
irv |= STRTOG_Inexhi;
}
irv |= STRTOG_Inexhi;
}
else
irv |= STRTOG_Inexlo;
}
}
*bp = b;
*exp = e;
return irv;
}
}

View File

@ -31,12 +31,7 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
void
#ifdef KR_headers
rshift(b, k) Bigint *b; int k;
#else
rshift(Bigint *b, int k)
#endif
void rshift (Bigint *b, int k)
{
ULong *x, *x1, *xe, y;
int n;
@ -52,24 +47,19 @@ rshift(Bigint *b, int k)
while(x < xe) {
*x1++ = (y | (*x << n)) & ALL_ON;
y = *x++ >> k;
}
}
if ((*x1 = y) !=0)
x1++;
}
}
else
while(x < xe)
*x1++ = *x++;
}
}
if ((b->wds = x1 - b->x) == 0)
b->x[0] = 0;
}
}
int
#ifdef KR_headers
trailz(b) Bigint *b;
#else
trailz(Bigint *b)
#endif
int trailz (Bigint *b)
{
ULong L, *x, *xe;
int n = 0;
@ -81,6 +71,6 @@ trailz(Bigint *b)
if (x < xe) {
L = *x;
n += lo0bits(&L);
}
return n;
}
return n;
}

View File

@ -31,25 +31,19 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
unsigned char hexdig[256];
unsigned char hexdig[256];
static void
#ifdef KR_headers
htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc;
#else
htinit(unsigned char *h, unsigned char *s, int inc)
#endif
static void htinit (unsigned char *h, unsigned char *s, int inc)
{
int i, j;
for(i = 0; (j = s[i]) !=0; i++)
h[j] = i + inc;
}
}
void
hexdig_init_D2A(Void)
void hexdig_init_D2A (void)
{
#define USC (unsigned char *)
htinit(hexdig, USC "0123456789", 0x10);
htinit(hexdig, USC "abcdef", 0x10 + 10);
htinit(hexdig, USC "ABCDEF", 0x10 + 10);
}
}

View File

@ -31,12 +31,7 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
static void
#ifdef KR_headers
L_shift(x, x1, i) ULong *x; ULong *x1; int i;
#else
L_shift(ULong *x, ULong *x1, int i)
#endif
static void L_shift (ULong *x, ULong *x1, int i)
{
int j;
@ -46,19 +41,13 @@ L_shift(ULong *x, ULong *x1, int i)
do {
*x |= x[1] << j;
x[1] >>= i;
} while(++x < x1);
}
} while(++x < x1);
}
int
#ifdef KR_headers
hexnan(sp, fpi, x0)
CONST char **sp; FPI *fpi; ULong *x0;
#else
hexnan( CONST char **sp, FPI *fpi, ULong *x0)
#endif
int hexnan (const char **sp, FPI *fpi, ULong *x0)
{
ULong c, h, *x, *x1, *xe;
CONST char *s;
const char *s;
int havedig, hd0, i, nbits;
if (!hexdig['0'])
@ -71,7 +60,13 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
x1 = xe = x;
havedig = hd0 = i = 0;
s = *sp;
while((c = *(CONST unsigned char*)++s)) {
/* allow optional initial 0x or 0X */
while((c = *(const unsigned char*)(s+1)) && c <= ' ')
++s;
if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
&& *(const unsigned char*)(s+3) > ' ')
s += 2;
while((c = *(const unsigned char*)++s)) {
if (!(h = hexdig[c])) {
if (c <= ' ') {
if (hd0 < havedig) {
@ -80,29 +75,42 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
if (x <= x0) {
i = 8;
continue;
}
}
hd0 = havedig;
*--x = 0;
x1 = x;
i = 0;
}
continue;
}
while(*(const unsigned char*)(s+1) <= ' ')
++s;
if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
&& *(const unsigned char*)(s+3) > ' ')
s += 2;
continue;
}
if (/*(*/ c == ')' && havedig) {
*sp = s + 1;
break;
}
return STRTOG_NaN;
}
#ifndef GDTOA_NON_PEDANTIC_NANCHECK
do {
if (/*(*/ c == ')') {
*sp = s + 1;
break;
}
} while((c = *++s));
#endif
return STRTOG_NaN;
}
havedig++;
if (++i > 8) {
if (x <= x0)
continue;
i = 1;
*--x = 0;
}
*x = (*x << 4) | (h & 0xf);
}
*x = (*x << 4) | (h & 0xf);
}
if (!havedig)
return STRTOG_NaN;
if (x < x1 && i < 8)
@ -113,19 +121,19 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
while(x <= xe);
do *x1++ = 0;
while(x1 <= xe);
}
}
else {
/* truncate high-order word if necessary */
if ( (i = nbits & (ULbits-1)) !=0)
*xe &= ((ULong)0xffffffff) >> (ULbits - i);
}
}
for(x1 = xe;; --x1) {
if (*x1 != 0)
break;
if (x1 == x0) {
*x1 = 1;
break;
}
}
return STRTOG_NaNbits;
}
return STRTOG_NaNbits;
}

View File

@ -30,16 +30,16 @@ THIS SOFTWARE.
* with " at " changed at "@" and " dot " changed to "."). */
#ifdef __MINGW32__
/* we have to include windows.h before gdtoa headers, otherwise
defines cause conflicts. */
#if defined(__MINGW32__) || defined(__MINGW64__)
/* we have to include windows.h before gdtoa
headers, otherwise defines cause conflicts. */
#define WIN32_LEAN_AND_MEAN
#include <windows.h>
#define NLOCKS 2
#ifdef USE_WIN32_SL
/* Use spin locks. */
/* Use spin locks. */
static long dtoa_sl[NLOCKS];
#define ACQUIRE_DTOA_LOCK(n) \
@ -47,7 +47,8 @@ static long dtoa_sl[NLOCKS];
Sleep (0);
#define FREE_DTOA_LOCK(n) InterlockedExchange (&dtoa_sl[n], 0);
#else
#else /* USE_WIN32_SL */
#include <stdlib.h>
static CRITICAL_SECTION dtoa_CritSec[NLOCKS];
static long dtoa_CS_init = 0;
@ -55,69 +56,59 @@ static long dtoa_CS_init = 0;
1 = initializing
2 = initialized
3 = deleted
*/
static void dtoa_lock_cleanup()
*/
static void dtoa_lock_cleanup (void)
{
long last_CS_init = InterlockedExchange (&dtoa_CS_init,3);
if (2 == last_CS_init)
{
int i;
for (i = 0; i < NLOCKS; i++)
DeleteCriticalSection (&dtoa_CritSec[i]);
}
long last_CS_init = InterlockedExchange (&dtoa_CS_init,3);
if (2 == last_CS_init) {
int i;
for (i = 0; i < NLOCKS; i++)
DeleteCriticalSection (&dtoa_CritSec[i]);
}
}
static void dtoa_lock(int n)
static void dtoa_lock (int n)
{
if (2 == dtoa_CS_init)
{
EnterCriticalSection (&dtoa_CritSec[n]);
return;
}
if (2 == dtoa_CS_init) {
EnterCriticalSection (&dtoa_CritSec[n]);
return;
}
else if (0 == dtoa_CS_init) {
long last_CS_init = InterlockedExchange (&dtoa_CS_init, 1);
if (0 == last_CS_init) {
int i;
for (i = 0; i < NLOCKS; i++)
InitializeCriticalSection (&dtoa_CritSec[i]);
atexit (dtoa_lock_cleanup);
dtoa_CS_init = 2;
}
else if (2 == last_CS_init)
dtoa_CS_init = 2;
}
/* Another thread is initializing. Wait. */
while (1 == dtoa_CS_init)
Sleep (1);
else if (0 == dtoa_CS_init)
{
long last_CS_init = InterlockedExchange (&dtoa_CS_init, 1);
if (0 == last_CS_init)
{
int i;
for (i = 0; i < NLOCKS; i++)
InitializeCriticalSection (&dtoa_CritSec[i]);
atexit (dtoa_lock_cleanup);
dtoa_CS_init = 2;
}
else if (2 == last_CS_init)
dtoa_CS_init = 2;
}
/* Another thread is initializing. Wait. */
while (1 == dtoa_CS_init)
Sleep (1);
/* It had better be initialized now. */
if (2 == dtoa_CS_init)
EnterCriticalSection(&dtoa_CritSec[n]);
/* It had better be initialized now. */
if (2 == dtoa_CS_init)
EnterCriticalSection(&dtoa_CritSec[n]);
}
static void dtoa_unlock(int n)
static void dtoa_unlock (int n)
{
if (2 == dtoa_CS_init)
LeaveCriticalSection (&dtoa_CritSec[n]);
if (2 == dtoa_CS_init)
LeaveCriticalSection (&dtoa_CritSec[n]);
}
#define ACQUIRE_DTOA_LOCK(n) dtoa_lock(n)
#define FREE_DTOA_LOCK(n) dtoa_unlock(n)
#endif
#endif /* USE_WIN32_SL */
#endif /* __MINGW32__ */
#endif /* __MINGW32__ / __MINGW64__ */
#include "gdtoaimp.h"
#ifndef MULTIPLE_THREADS
char *dtoa_result;
#endif
static Bigint *freelist[Kmax+1];
static Bigint *freelist[Kmax+1];
#ifndef Omit_Private_Memory
#ifndef PRIVATE_MEM
#define PRIVATE_MEM 2304
@ -126,13 +117,7 @@ static void dtoa_unlock(int n)
static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
#endif
Bigint *
Balloc
#ifdef KR_headers
(k) int k;
#else
(int k)
#endif
Bigint *Balloc (int k)
{
int x;
Bigint *rv;
@ -141,9 +126,11 @@ Balloc
#endif
ACQUIRE_DTOA_LOCK(0);
if ( (rv = freelist[k]) !=0) {
/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
/* but this case seems very unlikely. */
if (k <= Kmax && (rv = freelist[k]) !=0) {
freelist[k] = rv->next;
}
}
else {
x = 1 << k;
#ifdef Omit_Private_Memory
@ -151,51 +138,44 @@ Balloc
#else
len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
/sizeof(double);
if (pmem_next - private_mem + len <= PRIVATE_mem) {
if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
rv = (Bigint*)pmem_next;
pmem_next += len;
}
}
else
rv = (Bigint*)MALLOC(len*sizeof(double));
#endif
rv->k = k;
rv->maxwds = x;
}
}
FREE_DTOA_LOCK(0);
rv->sign = rv->wds = 0;
return rv;
}
}
void
Bfree
#ifdef KR_headers
(v) Bigint *v;
#else
(Bigint *v)
#endif
void Bfree (Bigint *v)
{
if (v) {
ACQUIRE_DTOA_LOCK(0);
v->next = freelist[v->k];
freelist[v->k] = v;
FREE_DTOA_LOCK(0);
if (v->k > Kmax)
free((void*)v);
else {
ACQUIRE_DTOA_LOCK(0);
v->next = freelist[v->k];
freelist[v->k] = v;
FREE_DTOA_LOCK(0);
}
}
}
// Shift y so lowest bit is 1 and return the number of bits y was
// shifted.
// With __GNUC__, we use an inline wrapper for __builtin_clz()
/* lo0bits(): Shift y so lowest bit is 1 and return the
* number of bits y was shifted.
* With GCC, we use an inline wrapper for __builtin_clz()
*/
#ifndef __GNUC__
int
lo0bits
#ifdef KR_headers
(y) ULong *y;
#else
(ULong *y)
#endif
int lo0bits (ULong *y)
{
register int k;
register ULong x = *y;
int k;
ULong x = *y;
if (x & 7) {
if (x & 1)
@ -203,46 +183,39 @@ lo0bits
if (x & 2) {
*y = x >> 1;
return 1;
}
}
*y = x >> 2;
return 2;
}
}
k = 0;
if (!(x & 0xffff)) {
k = 16;
x >>= 16;
}
}
if (!(x & 0xff)) {
k += 8;
x >>= 8;
}
}
if (!(x & 0xf)) {
k += 4;
x >>= 4;
}
}
if (!(x & 0x3)) {
k += 2;
x >>= 2;
}
}
if (!(x & 1)) {
k++;
x >>= 1;
if (!x)
return 32;
}
}
*y = x;
return k;
}
}
#endif /* __GNUC__ */
#endif /* __GNUC__ */
Bigint *
multadd
#ifdef KR_headers
(b, m, a) Bigint *b; int m, a;
#else
(Bigint *b, int m, int a) /* multiply by m and add a */
#endif
Bigint *multadd (Bigint *b, int m, int a) /* multiply by m and add a */
{
int i, wds;
#ifdef ULLong
@ -278,66 +251,54 @@ multadd
*x++ = y & 0xffff;
#endif
#endif
}
while(++i < wds);
} while(++i < wds);
if (carry) {
if (wds >= b->maxwds) {
b1 = Balloc(b->k+1);
Bcopy(b1, b);
Bfree(b);
b = b1;
}
}
b->x[wds++] = carry;
b->wds = wds;
}
return b;
}
return b;
}
// With __GNUC__, we use an inline wrapper for __builtin_clz()
/* hi0bits();
* With GCC, we use an inline wrapper for __builtin_clz()
*/
#ifndef __GNUC__
int
hi0bits_D2A
#ifdef KR_headers
(x) register ULong x;
#else
(register ULong x)
#endif
int hi0bits_D2A (ULong x)
{
register int k = 0;
int k = 0;
if (!(x & 0xffff0000)) {
k = 16;
x <<= 16;
}
}
if (!(x & 0xff000000)) {
k += 8;
x <<= 8;
}
}
if (!(x & 0xf0000000)) {
k += 4;
x <<= 4;
}
}
if (!(x & 0xc0000000)) {
k += 2;
x <<= 2;
}
}
if (!(x & 0x80000000)) {
k++;
if (!(x & 0x40000000))
return 32;
}
return k;
}
return k;
}
#endif /* __GNUC__ */
#endif
Bigint *
i2b
#ifdef KR_headers
(i) int i;
#else
(int i)
#endif
Bigint *i2b (int i)
{
Bigint *b;
@ -345,15 +306,9 @@ i2b
b->x[0] = i;
b->wds = 1;
return b;
}
}
Bigint *
mult
#ifdef KR_headers
(a, b) Bigint *a, *b;
#else
(Bigint *a, Bigint *b)
#endif
Bigint *mult (Bigint *a, Bigint *b)
{
Bigint *c;
int k, wa, wb, wc;
@ -372,7 +327,7 @@ mult
c = a;
a = b;
b = c;
}
}
k = a->k;
wa = a->wds;
wb = b->wds;
@ -397,11 +352,10 @@ mult
z = *x++ * (ULLong)y + *xc + carry;
carry = z >> 32;
*xc++ = z & 0xffffffffUL;
}
while(x < xae);
} while(x < xae);
*xc = carry;
}
}
}
#else
#ifdef Pack_32
for(; xb < xbe; xb++, xc0++) {
@ -415,10 +369,9 @@ mult
z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
carry = z2 >> 16;
Storeinc(xc, z2, z);
}
while(x < xae);
} while(x < xae);
*xc = carry;
}
}
if ( (y = *xb >> 16) !=0) {
x = xa;
xc = xc0;
@ -430,11 +383,10 @@ mult
Storeinc(xc, z, z2);
z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
carry = z2 >> 16;
}
while(x < xae);
} while(x < xae);
*xc = z2;
}
}
}
#else
for(; xb < xbe; xc0++) {
if ( (y = *xb++) !=0) {
@ -445,27 +397,20 @@ mult
z = *x++ * y + *xc + carry;
carry = z >> 16;
*xc++ = z & 0xffff;
}
while(x < xae);
} while(x < xae);
*xc = carry;
}
}
}
#endif
#endif
for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
c->wds = wc;
return c;
}
}
static Bigint *p5s;
static Bigint *p5s;
Bigint *
pow5mult
#ifdef KR_headers
(b, k) Bigint *b; int k;
#else
(Bigint *b, int k)
#endif
Bigint *pow5mult (Bigint *b, int k)
{
Bigint *b1, *p5, *p51;
int i;
@ -483,19 +428,19 @@ pow5mult
if (!(p5 = p5s)) {
p5 = p5s = i2b(625);
p5->next = 0;
}
}
FREE_DTOA_LOCK(1);
#else
p5 = p5s = i2b(625);
p5->next = 0;
#endif
}
}
for(;;) {
if (k & 1) {
b1 = mult(b, p5);
Bfree(b);
b = b1;
}
}
if (!(k >>= 1))
break;
if ((p51 = p5->next) == 0) {
@ -504,26 +449,19 @@ pow5mult
if (!(p51 = p5->next)) {
p51 = p5->next = mult(p5,p5);
p51->next = 0;
}
}
FREE_DTOA_LOCK(1);
#else
p51 = p5->next = mult(p5,p5);
p51->next = 0;
#endif
}
p5 = p51;
}
return b;
p5 = p51;
}
return b;
}
Bigint *
lshift
#ifdef KR_headers
(b, k) Bigint *b; int k;
#else
(Bigint *b, int k)
#endif
Bigint *lshift (Bigint *b, int k)
{
int i, k1, n, n1;
Bigint *b1;
@ -547,8 +485,7 @@ lshift
do {
*x1++ = *x << k | z;
z = *x++ >> k1;
}
while(x < xe);
} while(x < xe);
if ((*x1 = z) !=0)
++n1;
#else
@ -557,27 +494,20 @@ lshift
do {
*x1++ = *x << k & 0xffff | z;
z = *x++ >> k1;
}
while(x < xe);
} while(x < xe);
if (*x1 = z)
++n1;
#endif
}
}
else do
*x1++ = *x++;
while(x < xe);
b1->wds = n1 - 1;
Bfree(b);
return b1;
}
}
int
cmp
#ifdef KR_headers
(a, b) Bigint *a, *b;
#else
(Bigint *a, Bigint *b)
#endif
int cmp (Bigint *a, Bigint *b)
{
ULong *xa, *xa0, *xb, *xb0;
int i, j;
@ -601,17 +531,11 @@ cmp
return *xa < *xb ? -1 : 1;
if (xa <= xa0)
break;
}
return 0;
}
return 0;
}
Bigint *
diff
#ifdef KR_headers
(a, b) Bigint *a, *b;
#else
(Bigint *a, Bigint *b)
#endif
Bigint *diff (Bigint *a, Bigint *b)
{
Bigint *c;
int i, wa, wb;
@ -631,13 +555,13 @@ diff
c->wds = 1;
c->x[0] = 0;
return c;
}
}
if (i < 0) {
c = a;
a = b;
b = c;
i = 1;
}
}
else
i = 0;
c = Balloc(a->k);
@ -655,13 +579,12 @@ diff
y = (ULLong)*xa++ - *xb++ - borrow;
borrow = y >> 32 & 1UL;
*xc++ = y & 0xffffffffUL;
}
while(xb < xbe);
} while(xb < xbe);
while(xa < xae) {
y = *xa++ - borrow;
borrow = y >> 32 & 1UL;
*xc++ = y & 0xffffffffUL;
}
}
#else
#ifdef Pack_32
do {
@ -670,52 +593,40 @@ diff
z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
borrow = (z & 0x10000) >> 16;
Storeinc(xc, z, y);
}
while(xb < xbe);
} while(xb < xbe);
while(xa < xae) {
y = (*xa & 0xffff) - borrow;
borrow = (y & 0x10000) >> 16;
z = (*xa++ >> 16) - borrow;
borrow = (z & 0x10000) >> 16;
Storeinc(xc, z, y);
}
}
#else
do {
y = *xa++ - *xb++ - borrow;
borrow = (y & 0x10000) >> 16;
*xc++ = y & 0xffff;
}
while(xb < xbe);
} while(xb < xbe);
while(xa < xae) {
y = *xa++ - borrow;
borrow = (y & 0x10000) >> 16;
*xc++ = y & 0xffff;
}
}
#endif
#endif
while(!*--xc)
wa--;
c->wds = wa;
return c;
}
}
double
b2d
#ifdef KR_headers
(a, e) Bigint *a; int *e;
#else
(Bigint *a, int *e)
#endif
double b2d (Bigint *a, int *e)
{
ULong *xa, *xa0, w, y, z;
int k;
double d;
#ifdef VAX
ULong d0, d1;
#else
#define d0 word0(d)
#define d1 word1(d)
#endif
union _dbl_union d;
#define d0 word0(&d)
#define d1 word1(&d)
xa0 = a->x;
xa = xa0 + a->wds;
@ -731,17 +642,17 @@ b2d
w = xa > xa0 ? *--xa : 0;
d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
goto ret_d;
}
}
z = xa > xa0 ? *--xa : 0;
if (k -= Ebits) {
d0 = Exp_1 | y << k | z >> (32 - k);
y = xa > xa0 ? *--xa : 0;
d1 = z << k | y >> (32 - k);
}
}
else {
d0 = Exp_1 | y;
d1 = z;
}
}
#else
if (k < Ebits + 16) {
z = xa > xa0 ? *--xa : 0;
@ -750,7 +661,7 @@ b2d
y = xa > xa0 ? *--xa : 0;
d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
goto ret_d;
}
}
z = xa > xa0 ? *--xa : 0;
w = xa > xa0 ? *--xa : 0;
k -= Ebits + 16;
@ -759,37 +670,23 @@ b2d
d1 = w << k + 16 | y << k;
#endif
ret_d:
#ifdef VAX
word0(d) = d0 >> 16 | d0 << 16;
word1(d) = d1 >> 16 | d1 << 16;
#endif
return dval(d);
}
return dval(&d);
#undef d0
#undef d1
}
Bigint *
d2b
#ifdef KR_headers
(d, e, bits) double d; int *e, *bits;
#else
(double d, int *e, int *bits)
#endif
Bigint *d2b (double dd, int *e, int *bits)
{
Bigint *b;
union _dbl_union d;
#ifndef Sudden_Underflow
int i;
#endif
int de, k;
ULong *x, y, z;
#ifdef VAX
ULong d0, d1;
d0 = word0(d) >> 16 | word0(d) << 16;
d1 = word1(d) >> 16 | word1(d) << 16;
#else
#define d0 word0(d)
#define d1 word1(d)
#endif
#define d0 word0(&d)
#define d1 word1(&d)
d.d = dd;
#ifdef Pack_32
b = Balloc(1);
@ -802,9 +699,7 @@ d2b
d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
#ifdef Sudden_Underflow
de = (int)(d0 >> Exp_shift);
#ifndef IBM
z |= Exp_msk11;
#endif
#else
if ( (de = (int)(d0 >> Exp_shift)) !=0)
z |= Exp_msk1;
@ -814,19 +709,15 @@ d2b
if ( (k = lo0bits(&y)) !=0) {
x[0] = y | z << (32 - k);
z >>= k;
}
}
else
x[0] = y;
#ifndef Sudden_Underflow
i =
#endif
b->wds = (x[1] = z) !=0 ? 2 : 1;
}
}
else {
#ifdef DEBUG
if (!z)
Bug("Zero passed to d2b");
#endif
k = lo0bits(&z);
x[0] = z;
#ifndef Sudden_Underflow
@ -834,7 +725,7 @@ d2b
#endif
b->wds = 1;
k += 32;
}
}
#else
if ( (y = d1) !=0) {
if ( (k = lo0bits(&y)) !=0)
@ -843,22 +734,22 @@ d2b
x[1] = z >> k - 16 & 0xffff;
x[2] = z >> k;
i = 2;
}
}
else {
x[0] = y & 0xffff;
x[1] = y >> 16 | z << 16 - k & 0xffff;
x[2] = z >> k & 0xffff;
x[3] = z >> k+16;
i = 3;
}
}
else {
x[0] = y & 0xffff;
x[1] = y >> 16;
x[2] = z & 0xffff;
x[3] = z >> 16;
i = 3;
}
}
}
else {
#ifdef DEBUG
if (!z)
@ -868,14 +759,14 @@ d2b
if (k >= 16) {
x[0] = z;
i = 0;
}
}
else {
x[0] = z & 0xffff;
x[1] = z >> 16;
i = 1;
}
k += 32;
}
k += 32;
}
while(!x[i])
--i;
b->wds = i + 1;
@ -883,15 +774,10 @@ d2b
#ifndef Sudden_Underflow
if (de) {
#endif
#ifdef IBM
*e = (de - Bias - (P-1) << 2) + k;
*bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
#else
*e = de - Bias - (P-1) + k;
*bits = P - k;
#endif
#ifndef Sudden_Underflow
}
}
else {
*e = de - Bias - (P-1) + 1 + k;
#ifdef Pack_32
@ -899,64 +785,39 @@ d2b
#else
*bits = (i+2)*16 - hi0bits(x[i]);
#endif
}
}
#endif
return b;
}
#undef d0
#undef d1
}
CONST double
#ifdef IEEE_Arith
const double
bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
};
#else
#ifdef IBM
bigtens[] = { 1e16, 1e32, 1e64 };
CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
#else
bigtens[] = { 1e16, 1e32 };
CONST double tinytens[] = { 1e-16, 1e-32 };
#endif
#endif
const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
CONST double
const double
tens[] = {
1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1e20, 1e21, 1e22
#ifdef VAX
, 1e23, 1e24
#endif
};
};
char *
#ifdef KR_headers
strcp_D2A(a, b) char *a; char *b;
#else
strcp_D2A(char *a, CONST char *b)
#endif
char *strcp_D2A (char *a, const char *b)
{
while((*a = *b++))
a++;
return a;
}
}
#ifdef NO_STRING_H
Char *
#ifdef KR_headers
memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
#else
memcpy_D2A(void *a1, void *b1, size_t len)
#endif
void *memcpy_D2A (void *a1, void *b1, size_t len)
{
register char *a = (char*)a1, *ae = a + len;
register char *b = (char*)b1, *a0 = a;
char *a = (char*)a1, *ae = a + len;
char *b = (char*)b1, *a0 = a;
while(a < ae)
*a++ = *b++;
return a0;
}
}
#endif /* NO_STRING_H */

View File

@ -74,12 +74,16 @@ main(void)
double d;
Ulong L[4];
#ifndef NO_LONG_LONG
unsigned short u[5];
/* need u[8] instead of u[5] for 64 bit */
unsigned short u[8];
long double D;
#endif
} U;
U a, b, c;
int i;
a.L[0]=a.L[1]=a.L[2]=a.L[3]=0;
b.L[0]=b.L[1]=b.L[2]=b.L[3]=0;
c.L[0]=c.L[1]=c.L[2]=c.L[3]=0;
a.L[0] = b.L[0] = 0x7f800000;
c.f = a.f - b.f;

View File

@ -31,13 +31,7 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
Bigint *
s2b
#ifdef KR_headers
(s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
#else
(CONST char *s, int nd0, int nd, ULong y9)
#endif
Bigint *s2b (const char *s, int nd0, int nd, ULong y9, int dplen)
{
Bigint *b;
int i, k;
@ -60,82 +54,51 @@ s2b
s += 9;
do b = multadd(b, 10, *s++ - '0');
while(++i < nd0);
s++;
}
s += dplen;
}
else
s += 10;
s += dplen + 9;
for(; i < nd; i++)
b = multadd(b, 10, *s++ - '0');
return b;
}
}
double
ratio
#ifdef KR_headers
(a, b) Bigint *a, *b;
#else
(Bigint *a, Bigint *b)
#endif
double ratio (Bigint *a, Bigint *b)
{
double da, db;
union _dbl_union da, db;
int k, ka, kb;
dval(da) = b2d(a, &ka);
dval(db) = b2d(b, &kb);
dval(&da) = b2d(a, &ka);
dval(&db) = b2d(b, &kb);
k = ka - kb + ULbits*(a->wds - b->wds);
#ifdef IBM
if (k > 0) {
word0(da) += (k >> 2)*Exp_msk1;
if (k &= 3)
dval(da) *= 1 << k;
}
else {
k = -k;
word0(db) += (k >> 2)*Exp_msk1;
if (k &= 3)
dval(db) *= 1 << k;
}
#else
if (k > 0)
word0(da) += k*Exp_msk1;
word0(&da) += k*Exp_msk1;
else {
k = -k;
word0(db) += k*Exp_msk1;
}
#endif
return dval(da) / dval(db);
word0(&db) += k*Exp_msk1;
}
return dval(&da) / dval(&db);
}
#ifdef INFNAN_CHECK
int
match
#ifdef KR_headers
(sp, t) char **sp, *t;
#else
(CONST char **sp, char *t)
#endif
int match (const char **sp, char *t)
{
int c, d;
CONST char *s = *sp;
const char *s = *sp;
while( (d = *t++) !=0) {
if ((c = *++s) >= 'A' && c <= 'Z')
c += 'a' - 'A';
if (c != d)
return 0;
}
}
*sp = s + 1;
return 1;
}
}
#endif /* INFNAN_CHECK */
void
#ifdef KR_headers
copybits(c, n, b) ULong *c; int n; Bigint *b;
#else
copybits(ULong *c, int n, Bigint *b)
#endif
void copybits (ULong *c, int n, Bigint *b)
{
ULong *ce, *x, *xe;
#ifdef Pack_16
@ -158,14 +121,9 @@ copybits(ULong *c, int n, Bigint *b)
#endif
while(c < ce)
*c++ = 0;
}
}
ULong
#ifdef KR_headers
any_on(b, k) Bigint *b; int k;
#else
any_on(Bigint *b, int k)
#endif
ULong any_on (Bigint *b, int k)
{
int n, nwds;
ULong *x, *x0, x1, x2;
@ -181,11 +139,11 @@ any_on(Bigint *b, int k)
x1 <<= k;
if (x1 != x2)
return 1;
}
}
x0 = x;
x += n;
while(x > x0)
if (*--x)
return 1;
return 0;
}
}

View File

@ -35,21 +35,13 @@ THIS SOFTWARE.
#include "locale.h"
#endif
static CONST int
static const int
fivesbits[] = { 0, 3, 5, 7, 10, 12, 14, 17, 19, 21,
24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
47, 49, 52
#ifdef VAX
, 54, 56
#endif
};
};
Bigint *
#ifdef KR_headers
increment(b) Bigint *b;
#else
increment(Bigint *b)
#endif
Bigint *increment (Bigint *b)
{
ULong *x, *xe;
Bigint *b1;
@ -64,9 +56,9 @@ increment(Bigint *b)
if (*x < (ULong)0xffffffffL) {
++*x;
return b;
}
}
*x++ = 0;
} while(x < xe);
} while(x < xe);
#else
do {
y = *x + carry;
@ -74,7 +66,7 @@ increment(Bigint *b)
*x++ = y & 0xffff;
if (!carry)
return b;
} while(x < xe);
} while(x < xe);
if (carry)
#endif
{
@ -83,18 +75,13 @@ increment(Bigint *b)
Bcopy(b1,b);
Bfree(b);
b = b1;
}
b->x[b->wds++] = 1;
}
return b;
b->x[b->wds++] = 1;
}
return b;
}
int
#ifdef KR_headers
decrement(b) Bigint *b;
#else
decrement(Bigint *b)
#endif
void decrement (Bigint *b)
{
ULong *x, *xe;
#ifdef Pack_16
@ -108,26 +95,19 @@ decrement(Bigint *b)
if (*x) {
--*x;
break;
}
*x++ = 0xffffffffL;
}
while(x < xe);
*x++ = 0xffffffffL;
} while(x < xe);
#else
do {
y = *x - borrow;
borrow = (y & 0x10000) >> 16;
*x++ = y & 0xffff;
} while(borrow && x < xe);
} while(borrow && x < xe);
#endif
return STRTOG_Inexlo;
}
}
static int
#ifdef KR_headers
all_on(b, n) Bigint *b; int n;
#else
all_on(Bigint *b, int n)
#endif
static int all_on (Bigint *b, int n)
{
ULong *x, *xe;
@ -139,14 +119,9 @@ all_on(Bigint *b, int n)
if (n &= kmask)
return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON;
return 1;
}
}
Bigint *
#ifdef KR_headers
set_ones(b, n) Bigint *b; int n;
#else
set_ones(Bigint *b, int n)
#endif
Bigint *set_ones (Bigint *b, int n)
{
int k;
ULong *x, *xe;
@ -155,7 +130,7 @@ set_ones(Bigint *b, int n)
if (b->k < k) {
Bfree(b);
b = Balloc(k);
}
}
k = n >> kshift;
if (n &= kmask)
k++;
@ -167,30 +142,24 @@ set_ones(Bigint *b, int n)
if (n)
x[-1] >>= ULbits - n;
return b;
}
}
static int
rvOK
#ifdef KR_headers
(d, fpi, exp, bits, exact, rd, irv)
double d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
#else
(double d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
#endif
static int rvOK (dbl_union *d, FPI *fpi, Long *exp, ULong *bits,
int exact, int rd, int *irv)
{
Bigint *b;
ULong carry, inex, lostbits;
int bdif, e, j, k, k1, nb, rv;
carry = rv = 0;
b = d2b(d, &e, &bdif);
b = d2b(dval(d), &e, &bdif);
bdif -= nb = fpi->nbits;
e += bdif;
if (bdif <= 0) {
if (exact)
goto trunc;
goto ret;
}
}
if (P == nb) {
if (
#ifndef IMPRECISE_INEXACT
@ -204,11 +173,11 @@ rvOK
#endif
) goto trunc;
goto ret;
}
}
switch(rd) {
case 1:
case 1: /* round down (toward -Infinity) */
goto trunc;
case 2:
case 2: /* round up (toward +Infinity) */
break;
default: /* round near */
k = bdif - 1;
@ -220,11 +189,11 @@ rvOK
if (b->x[0] & 2)
break;
goto trunc;
}
}
if (b->x[k>>kshift] & ((ULong)1 << (k & kmask)))
break;
goto trunc;
}
}
/* "break" cases: round up 1 bit, then truncate; bdif > 0 */
carry = 1;
trunc:
@ -243,9 +212,9 @@ rvOK
lostbits = b->x[0] & 1;
rshift(b, 1);
e++;
}
}
}
}
else if (bdif < 0)
b = lshift(b, -bdif);
if (e < fpi->emin) {
@ -254,7 +223,7 @@ rvOK
if (k > nb || fpi->sudden_underflow) {
b->wds = inex = 0;
*irv = STRTOG_Underflow | STRTOG_Inexlo;
}
}
else {
k1 = k - 1;
if (k1 > 0 && !lostbits)
@ -268,19 +237,17 @@ rvOK
if (carry) {
b = increment(b);
inex = STRTOG_Inexhi | STRTOG_Underflow;
}
}
else if (lostbits)
inex = STRTOG_Inexlo | STRTOG_Underflow;
}
}
}
else if (e > fpi->emax) {
e = fpi->emax + 1;
*irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
#ifndef NO_ERRNO
errno = ERANGE;
#endif
SET_ERRNO(ERANGE);
b->wds = inex = 0;
}
}
*exp = e;
copybits(bits, nb, b);
*irv |= inex;
@ -288,54 +255,55 @@ rvOK
ret:
Bfree(b);
return rv;
}
}
static int
#ifdef KR_headers
mantbits(d) double d;
#else
mantbits(double d)
#endif
static int mantbits (dbl_union *d)
{
ULong L;
#ifdef VAX
L = word1(d) << 16 | word1(d) >> 16;
if (L)
#else
if ( (L = word1(d)) !=0)
#endif
return P - lo0bits(&L);
#ifdef VAX
L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11;
#else
L = word0(d) | Exp_msk1;
#endif
return P - 32 - lo0bits(&L);
}
}
int
__strtodg
#ifdef KR_headers
(s00, se, fpi, exp, bits)
CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits;
#else
(CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
#endif
int __strtodg (const char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
{
int abe, abits, asub;
int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm;
int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
int sudden_underflow;
CONST char *s, *s0, *s1;
double adj, adj0, rv, tol;
const char *s, *s0, *s1;
double adj0, tol;
Long L;
ULong y, z;
union _dbl_union adj, rv;
ULong *b, *be, y, z;
Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
#ifdef USE_LOCALE /*{{*/
#ifdef NO_LOCALE_CACHE
char *decimalpoint = localeconv()->decimal_point;
int dplen = strlen(decimalpoint);
#else
char *decimalpoint;
static char *decimalpoint_cache;
static int dplen;
if (!(s0 = decimalpoint_cache)) {
s0 = localeconv()->decimal_point;
if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
strcpy(decimalpoint_cache, s0);
s0 = decimalpoint_cache;
}
dplen = strlen(s0);
}
decimalpoint = (char*)s0;
#endif /*NO_LOCALE_CACHE*/
#else /*USE_LOCALE}{*/
#define dplen 1
#endif /*USE_LOCALE}}*/
irv = STRTOG_Zero;
denorm = sign = nz0 = nz = 0;
dval(rv) = 0.;
dval(&rv) = 0.;
rvb = 0;
nbits = fpi->nbits;
for(s = s00;;s++) switch(*s) {
@ -360,7 +328,7 @@ __strtodg
continue;
default:
goto break2;
}
}
break2:
if (*s == '0') {
#ifndef NO_HEX_FP
@ -371,15 +339,15 @@ __strtodg
if (irv == STRTOG_NoNumber) {
s = s00;
sign = 0;
}
}
goto ret;
}
}
#endif
nz0 = 1;
while(*++s == '0') ;
if (!*s)
goto ret;
}
}
sudden_underflow = fpi->sudden_underflow;
s0 = s;
y = z = 0;
@ -390,13 +358,17 @@ __strtodg
z = 10*z + c - '0';
nd0 = nd;
#ifdef USE_LOCALE
if (c == *localeconv()->decimal_point)
if (c == *decimalpoint) {
for(i = 1; decimalpoint[i]; ++i)
if (s[i] != decimalpoint[i])
goto dig_done;
s += i;
c = *s;
#else
if (c == '.')
#endif
{
decpt = 1;
if (c == '.') {
c = *++s;
#endif
decpt = 1;
if (!nd) {
for(; c == '0'; c = *++s)
nz++;
@ -405,9 +377,9 @@ __strtodg
nf += nz;
nz = 0;
goto have_dig;
}
goto dig_done;
}
goto dig_done;
}
for(; c >= '0' && c <= '9'; c = *++s) {
have_dig:
nz++;
@ -423,9 +395,9 @@ __strtodg
else if (nd <= DBL_DIG + 1)
z = 10*z + c;
nz = 0;
}
}
}
}/*}*/
dig_done:
e = 0;
if (c == 'e' || c == 'E') {
@ -433,7 +405,7 @@ __strtodg
irv = STRTOG_NoNumber;
s = s00;
goto ret;
}
}
s00 = s;
esign = 0;
switch(c = *++s) {
@ -441,7 +413,7 @@ __strtodg
esign = 1;
case '+':
c = *++s;
}
}
if (c >= '0' && c <= '9') {
while(c == '0')
c = *++s;
@ -459,13 +431,13 @@ __strtodg
e = (int)L;
if (esign)
e = -e;
}
}
else
e = 0;
}
}
else
s = s00;
}
}
if (!nd) {
if (!nz && !nz0) {
#ifdef INFNAN_CHECK
@ -480,7 +452,7 @@ __strtodg
++s;
irv = STRTOG_Infinite;
goto infnanexp;
}
}
break;
case 'n':
case 'N':
@ -492,14 +464,14 @@ __strtodg
irv = hexnan(&s, fpi, bits);
#endif
goto infnanexp;
}
}
}
}
#endif /* INFNAN_CHECK */
irv = STRTOG_NoNumber;
s = s00;
}
goto ret;
}
goto ret;
}
irv = STRTOG_Normal;
e1 = e -= nf;
@ -513,7 +485,7 @@ __strtodg
break;
case FPI_Round_down:
rd = 1 + sign;
}
}
/* Now we have nd0 digits, starting at s0, followed by a
* decimal point, followed by nd-nd0 digits. The number we're
@ -523,28 +495,24 @@ __strtodg
if (!nd0)
nd0 = nd;
k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
dval(rv) = y;
dval(&rv) = y;
if (k > 9)
dval(rv) = tens[k - 9] * dval(rv) + z;
dval(&rv) = tens[k - 9] * dval(&rv) + z;
bd0 = 0;
if (nbits <= P && nd <= DBL_DIG) {
if (!e) {
if (rvOK(dval(rv), fpi, exp, bits, 1, rd, &irv))
if (rvOK(&rv, fpi, exp, bits, 1, rd, &irv))
goto ret;
}
}
else if (e > 0) {
if (e <= Ten_pmax) {
#ifdef VAX
goto vax_ovfl_check;
#else
i = fivesbits[e] + mantbits(dval(rv)) <= P;
/* rv = */ rounded_product(dval(rv), tens[e]);
if (rvOK(dval(rv), fpi, exp, bits, i, rd, &irv))
i = fivesbits[e] + mantbits(&rv) <= P;
/* rv = */ rounded_product(dval(&rv), tens[e]);
if (rvOK(&rv, fpi, exp, bits, i, rd, &irv))
goto ret;
e1 -= e;
goto rv_notOK;
#endif
}
}
i = DBL_DIG - nd;
if (e <= Ten_pmax + i) {
/* A fancier test would sometimes let us do
@ -552,37 +520,22 @@ __strtodg
*/
e2 = e - i;
e1 -= i;
dval(rv) *= tens[i];
#ifdef VAX
/* VAX exponent range is so narrow we must
* worry about overflow here...
*/
vax_ovfl_check:
dval(adj) = dval(rv);
word0(adj) -= P*Exp_msk1;
/* adj = */ rounded_product(dval(adj), tens[e2]);
if ((word0(adj) & Exp_mask)
> Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
goto rv_notOK;
word0(adj) += P*Exp_msk1;
dval(rv) = dval(adj);
#else
/* rv = */ rounded_product(dval(rv), tens[e2]);
#endif
if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
dval(&rv) *= tens[i];
/* rv = */ rounded_product(dval(&rv), tens[e2]);
if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
goto ret;
e1 -= e2;
}
}
}
#ifndef Inaccurate_Divide
else if (e >= -Ten_pmax) {
/* rv = */ rounded_quotient(dval(rv), tens[-e]);
if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
/* rv = */ rounded_quotient(dval(&rv), tens[-e]);
if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
goto ret;
e1 -= e;
}
#endif
}
#endif
}
rv_notOK:
e1 += nd - k;
@ -591,62 +544,54 @@ __strtodg
e2 = 0;
if (e1 > 0) {
if ( (i = e1 & 15) !=0)
dval(rv) *= tens[i];
dval(&rv) *= tens[i];
if (e1 &= ~15) {
e1 >>= 4;
while(e1 >= (1 << (n_bigtens-1))) {
e2 += ((word0(rv) & Exp_mask)
e2 += ((word0(&rv) & Exp_mask)
>> Exp_shift1) - Bias;
word0(rv) &= ~Exp_mask;
word0(rv) |= Bias << Exp_shift1;
dval(rv) *= bigtens[n_bigtens-1];
word0(&rv) &= ~Exp_mask;
word0(&rv) |= Bias << Exp_shift1;
dval(&rv) *= bigtens[n_bigtens-1];
e1 -= 1 << (n_bigtens-1);
}
e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
word0(rv) &= ~Exp_mask;
word0(rv) |= Bias << Exp_shift1;
}
e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
word0(&rv) &= ~Exp_mask;
word0(&rv) |= Bias << Exp_shift1;
for(j = 0; e1 > 0; j++, e1 >>= 1)
if (e1 & 1)
dval(rv) *= bigtens[j];
}
dval(&rv) *= bigtens[j];
}
}
else if (e1 < 0) {
e1 = -e1;
if ( (i = e1 & 15) !=0)
dval(rv) /= tens[i];
dval(&rv) /= tens[i];
if (e1 &= ~15) {
e1 >>= 4;
while(e1 >= (1 << (n_bigtens-1))) {
e2 += ((word0(rv) & Exp_mask)
e2 += ((word0(&rv) & Exp_mask)
>> Exp_shift1) - Bias;
word0(rv) &= ~Exp_mask;
word0(rv) |= Bias << Exp_shift1;
dval(rv) *= tinytens[n_bigtens-1];
word0(&rv) &= ~Exp_mask;
word0(&rv) |= Bias << Exp_shift1;
dval(&rv) *= tinytens[n_bigtens-1];
e1 -= 1 << (n_bigtens-1);
}
e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
word0(rv) &= ~Exp_mask;
word0(rv) |= Bias << Exp_shift1;
}
e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
word0(&rv) &= ~Exp_mask;
word0(&rv) |= Bias << Exp_shift1;
for(j = 0; e1 > 0; j++, e1 >>= 1)
if (e1 & 1)
dval(rv) *= tinytens[j];
}
dval(&rv) *= tinytens[j];
}
#ifdef IBM
/* e2 is a correction to the (base 2) exponent of the return
* value, reflecting adjustments above to avoid overflow in the
* native arithmetic. For native IBM (base 16) arithmetic, we
* must multiply e2 by 4 to change from base 16 to 2.
*/
e2 <<= 2;
#endif
rvb = d2b(dval(rv), &rve, &rvbits); /* rv = rvb * 2^rve */
}
rvb = d2b(dval(&rv), &rve, &rvbits); /* rv = rvb * 2^rve */
rve += e2;
if ((j = rvbits - nbits) > 0) {
rshift(rvb, j);
rvbits = nbits;
rve += j;
}
}
bb0 = 0; /* trailing zero bits in rvb */
e2 = rve + rvbits - nbits;
if (e2 > fpi->emax + 1)
@ -658,7 +603,7 @@ __strtodg
if (j > 0) {
rvb = lshift(rvb, j);
rvbits += j;
}
}
else if (j < 0) {
rvbits += j;
if (rvbits <= 0) {
@ -669,22 +614,22 @@ __strtodg
*exp = emin;
irv = STRTOG_Underflow | STRTOG_Inexlo;
goto ret;
}
rvb->x[0] = rvb->wds = rvbits = 1;
}
rvb->x[0] = rvb->wds = rvbits = 1;
}
else
rshift(rvb, -j);
}
}
rve = rve1 = emin;
if (sudden_underflow && e2 + 1 < emin)
goto ufl;
}
}
/* Now the hard part -- adjusting rv to the correct value.*/
/* Put digits into bd: true value = bd * 10^e */
bd0 = s2b(s0, nd0, nd, y);
bd0 = s2b(s0, nd0, nd, y, dplen);
for(;;) {
bd = Balloc(bd0->k);
@ -698,11 +643,11 @@ __strtodg
if (e >= 0) {
bb2 = bb5 = 0;
bd2 = bd5 = e;
}
}
else {
bb2 = bb5 = -e;
bd2 = bd5 = 0;
}
}
if (bbe >= 0)
bb2 += bbe;
else
@ -721,13 +666,13 @@ __strtodg
bb2 -= i;
bd2 -= i;
bs2 -= i;
}
}
if (bb5 > 0) {
bs = pow5mult(bs, bb5);
bb1 = mult(bs, bb);
Bfree(bb);
bb = bb1;
}
}
bb2 -= bb0;
if (bb2 > 0)
bb = lshift(bb, bb2);
@ -754,7 +699,7 @@ __strtodg
if (dsign != 0) {
irv |= STRTOG_Inexhi;
goto adj1;
}
}
irv |= STRTOG_Inexlo;
if (rve1 == emin)
goto adj1;
@ -762,16 +707,16 @@ __strtodg
i++, j -= ULbits) {
if (rvb->x[i] & ALL_ON)
goto adj1;
}
}
if (j > 1 && lo0bits(rvb->x + i) < j - 1)
goto adj1;
rve = rve1 - 1;
rvb = set_ones(rvb, rvbits = nbits);
break;
}
}
irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
break;
}
}
if (i < 0) {
/* Error is less than half an ulp -- check for
* special case of mantissa a power of two.
@ -785,9 +730,9 @@ __strtodg
if (cmp(delta, bs) > 0) {
irv = STRTOG_Normal | STRTOG_Inexlo;
goto drop_down;
}
break;
}
break;
}
if (i == 0) {
/* exactly half-way between */
if (dsign) {
@ -799,9 +744,9 @@ __strtodg
irv = STRTOG_Normal | STRTOG_Inexhi;
denorm = 0;
break;
}
irv = STRTOG_Normal | STRTOG_Inexlo;
}
irv = STRTOG_Normal | STRTOG_Inexlo;
}
else if (bbbits == 1) {
irv = STRTOG_Normal;
drop_down:
@ -811,55 +756,53 @@ __strtodg
if (rvb->wds == 1 && rvb->x[0] == 1)
sudden_underflow = 1;
break;
}
}
rve -= nbits;
rvb = set_ones(rvb, rvbits = nbits);
break;
}
}
else
irv = STRTOG_Normal | STRTOG_Inexhi;
if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1))
break;
if (dsign) {
rvb = increment(rvb);
if ( (j = rvbits & kmask) !=0)
j = ULbits - j;
if (hi0bits(rvb->x[(rvb->wds - 1) >> kshift])
!= j)
j = kmask & (ULbits - (rvbits & kmask));
if (hi0bits(rvb->x[rvb->wds - 1]) != j)
rvbits++;
irv = STRTOG_Normal | STRTOG_Inexhi;
}
}
else {
if (bbbits == 1)
goto undfl;
decrement(rvb);
irv = STRTOG_Normal | STRTOG_Inexlo;
}
break;
}
if ((dval(adj) = ratio(delta, bs)) <= 2.) {
break;
}
if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
adj1:
inex = STRTOG_Inexlo;
if (dsign) {
asub = 0;
inex = STRTOG_Inexhi;
}
}
else if (denorm && bbbits <= 1) {
undfl:
rvb->wds = 0;
rve = emin;
irv = STRTOG_Underflow | STRTOG_Inexlo;
break;
}
adj0 = dval(adj) = 1.;
}
adj0 = dval(&adj) = 1.;
}
else {
adj0 = dval(adj) *= 0.5;
adj0 = dval(&adj) *= 0.5;
if (dsign) {
asub = 0;
inex = STRTOG_Inexlo;
}
if (dval(adj) < 2147483647.) {
}
if (dval(&adj) < 2147483647.) {
L = adj0;
adj0 -= L;
switch(rd) {
@ -876,22 +819,22 @@ __strtodg
inc_L:
L++;
inex = STRTOG_Inexact - inex;
}
}
dval(adj) = L;
}
}
dval(&adj) = L;
}
}
y = rve + rvbits;
/* adj *= ulp(dval(rv)); */
/* adj *= ulp(&rv); */
/* if (asub) rv -= adj; else rv += adj; */
if (!denorm && rvbits < nbits) {
rvb = lshift(rvb, j = nbits - rvbits);
rve -= j;
rvbits = nbits;
}
ab = d2b(dval(adj), &abe, &abits);
}
ab = d2b(dval(&adj), &abe, &abits);
if (abe < 0)
rshift(ab, -abe);
else if (abe > 0)
@ -911,15 +854,15 @@ __strtodg
if (rve1 == emin) {
--rvbits;
denorm = 1;
}
}
else {
rvb = lshift(rvb, 1);
--rve;
--rve1;
L = finished = 0;
}
}
}
}
else {
rvb = sum(rvb, ab);
k = rvb->wds - 1;
@ -928,15 +871,15 @@ __strtodg
if (denorm) {
if (++rvbits == nbits)
denorm = 0;
}
}
else {
rshift(rvb, 1);
rve++;
rve1++;
L = 0;
}
}
}
}
Bfree(ab);
Bfree(rvb0);
if (finished)
@ -945,32 +888,32 @@ __strtodg
z = rve + rvbits;
if (y == z && L) {
/* Can we stop now? */
tol = dval(adj) * 5e-16; /* > max rel error */
dval(adj) = adj0 - .5;
if (dval(adj) < -tol) {
tol = dval(&adj) * 5e-16; /* > max rel error */
dval(&adj) = adj0 - .5;
if (dval(&adj) < -tol) {
if (adj0 > tol) {
irv |= inex;
break;
}
}
else if (dval(adj) > tol && adj0 < 1. - tol) {
irv |= inex;
break;
}
}
else if (dval(&adj) > tol && adj0 < 1. - tol) {
irv |= inex;
break;
}
}
bb0 = denorm ? 0 : trailz(rvb);
Bfree(bb);
Bfree(bd);
Bfree(bs);
Bfree(delta);
}
}
if (!denorm && (j = nbits - rvbits)) {
if (j > 0)
rvb = lshift(rvb, j);
else
rshift(rvb, -j);
rve -= j;
}
}
*exp = rve;
Bfree(bb);
Bfree(bd);
@ -978,28 +921,52 @@ __strtodg
Bfree(bd0);
Bfree(delta);
if (rve > fpi->emax) {
switch(fpi->rounding & 3) {
case FPI_Round_near:
goto huge;
case FPI_Round_up:
if (!sign)
goto huge;
break;
case FPI_Round_down:
if (sign)
goto huge;
}
/* Round to largest representable magnitude */
Bfree(rvb);
rvb = 0;
irv = STRTOG_Normal | STRTOG_Inexlo;
*exp = fpi->emax;
b = bits;
be = b + ((fpi->nbits + 31) >> 5);
while(b < be)
*b++ = -1;
if ((j = fpi->nbits & 0x1f))
*--be >>= (32 - j);
goto ret;
huge:
rvb->wds = 0;
irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
#ifndef NO_ERRNO
errno = ERANGE;
#endif
SET_ERRNO(ERANGE);
infnanexp:
*exp = fpi->emax + 1;
}
}
ret:
if (denorm) {
if (sudden_underflow) {
rvb->wds = 0;
irv = STRTOG_Underflow | STRTOG_Inexlo;
}
SET_ERRNO(ERANGE);
}
else {
irv = (irv & ~STRTOG_Retmask) |
(rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
if (irv & STRTOG_Inexact)
if (irv & STRTOG_Inexact) {
irv |= STRTOG_Underflow;
SET_ERRNO(ERANGE);
}
}
}
if (se)
*se = (char *)s;
if (sign)
@ -1007,6 +974,6 @@ __strtodg
if (rvb) {
copybits(bits, nbits, rvb);
Bfree(rvb);
}
return irv;
}
return irv;
}

View File

@ -37,12 +37,7 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
double
#ifdef KR_headers
__strtod(s, sp) CONST char *s; char **sp;
#else
__strtod(CONST char *s, char **sp)
#endif
double __strtod (const char *s, char **sp)
{
static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
ULong bits[2];
@ -80,8 +75,8 @@ __strtod(CONST char *s, char **sp)
case STRTOG_NaNbits:
u.L[_0] = 0x7ff00000 | bits[1];
u.L[_1] = bits[0];
}
}
if (k & STRTOG_Neg)
u.L[_0] |= 0x80000000L;
return u.d;
}
}

View File

@ -31,20 +31,20 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
float
#ifdef KR_headers
__strtof(s, sp) CONST char *s; char **sp;
#else
__strtof(CONST char *s, char **sp)
#endif
float __strtof (const char *s, char **sp)
{
static FPI fpi = { 24, 1-127-24+1, 254-127-24+1, 1, SI };
static FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, SI };
ULong bits[1];
Long exp;
int k;
union { ULong L[1]; float f; } u;
union { ULong L[1]; float f; } u = { { 0 } };
#ifdef Honor_FLT_ROUNDS
#include "gdtoa_fltrnds.h"
#else
#define fpi &fpi0
#endif
k = __strtodg(s, sp, &fpi, &exp, bits);
k = __strtodg(s, sp, fpi, &exp, bits);
switch(k & STRTOG_Retmask) {
case STRTOG_NoNumber:
case STRTOG_Zero:
@ -53,7 +53,7 @@ __strtof(CONST char *s, char **sp)
case STRTOG_Normal:
case STRTOG_NaNbits:
u.L[0] = (bits[0] & 0x7fffff) | (exp + 0x7f + 23) << 23;
u.L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
break;
case STRTOG_Denormal:
@ -66,11 +66,11 @@ __strtof(CONST char *s, char **sp)
case STRTOG_NaN:
u.L[0] = f_QNAN;
}
}
if (k & STRTOG_Neg)
u.L[0] |= 0x80000000L;
return u.f;
}
}
float __cdecl
strtof (const char * __restrict__ src, char ** __restrict__ endptr)

View File

@ -51,20 +51,26 @@ THIS SOFTWARE.
#define _4 0
#endif
static int
#ifdef KR_headers
__strtopx(s, sp, V) CONST char *s; char **sp; long double *V;
#else
__strtopx(CONST char *s, char **sp, long double *V)
#endif
typedef union lD {
UShort L[5];
long double D;
} lD;
static int __strtopx (const char *s, char **sp, lD *V)
{
static FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
ULong bits[2];
Long exp;
int k;
UShort *L = (UShort*)V;
UShort *L = & (V->L[0]);
#ifdef Honor_FLT_ROUNDS
#include "gdtoa_fltrnds.h"
#else
#define fpi &fpi0
#endif
V->D = 0.0L;
k = __strtodg(s, sp, &fpi, &exp, bits);
k = __strtodg(s, sp, fpi, &exp, bits);
switch(k & STRTOG_Retmask) {
case STRTOG_NoNumber:
case STRTOG_Zero:
@ -97,23 +103,22 @@ __strtopx(CONST char *s, char **sp, long double *V)
L[2] = ldus_QNAN2;
L[3] = ldus_QNAN3;
L[4] = ldus_QNAN4;
}
}
if (k & STRTOG_Neg)
L[_0] |= 0x8000;
return k;
}
long double
__cdecl
__strtold (const char * __restrict__ src, char ** __restrict__ endptr)
{
long double ret;
__strtopx(src, endptr, &ret);
return ret;
}
long double
__cdecl
long double __cdecl
__strtold (const char * __restrict__ src, char ** __restrict__ endptr)
{
lD ret;
ret.D = 0.0L;
__strtopx(src, endptr, &ret);
return ret.D;
}
long double __cdecl
strtold (const char * __restrict__ src, char ** __restrict__ endptr)
__attribute__((alias("__strtold")));

View File

@ -31,12 +31,7 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
Bigint *
#ifdef KR_headers
sum(a, b) Bigint *a; Bigint *b;
#else
sum(Bigint *a, Bigint *b)
#endif
Bigint *sum (Bigint *a, Bigint *b)
{
Bigint *c;
ULong carry, *xc, *xa, *xb, *xe, y;
@ -46,7 +41,7 @@ sum(Bigint *a, Bigint *b)
if (a->wds < b->wds) {
c = b; b = a; a = c;
}
}
c = Balloc(a->k);
c->wds = a->wds;
carry = 0;
@ -61,8 +56,7 @@ sum(Bigint *a, Bigint *b)
z = (*xa++ >> 16) + (*xb++ >> 16) + carry;
carry = (z & 0x10000) >> 16;
Storeinc(xc, z, y);
}
while(xc < xe);
} while(xc < xe);
xe += a->wds - b->wds;
while(xc < xe) {
y = (*xa & 0xffff) + carry;
@ -70,20 +64,19 @@ sum(Bigint *a, Bigint *b)
z = (*xa++ >> 16) + carry;
carry = (z & 0x10000) >> 16;
Storeinc(xc, z, y);
}
}
#else
do {
y = *xa++ + *xb++ + carry;
carry = (y & 0x10000) >> 16;
*xc++ = y & 0xffff;
}
while(xc < xe);
} while(xc < xe);
xe += a->wds - b->wds;
while(xc < xe) {
y = *xa++ + carry;
carry = (y & 0x10000) >> 16;
*xc++ = y & 0xffff;
}
}
#endif
if (carry) {
if (c->wds == c->maxwds) {
@ -91,8 +84,8 @@ sum(Bigint *a, Bigint *b)
Bcopy(b, c);
Bfree(c);
c = b;
}
c->x[c->wds++] = 1;
}
return c;
c->x[c->wds++] = 1;
}
return c;
}

View File

@ -31,40 +31,31 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
double
ulp
#ifdef KR_headers
(x) double x;
#else
(double x)
#endif
double ulp (dbl_union *x)
{
Long L;
double a;
union _dbl_union a;
L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
#ifndef Sudden_Underflow
if (L > 0) {
#endif
#ifdef IBM
L |= Exp_msk1 >> 4;
#endif
word0(a) = L;
word1(a) = 0;
word0(&a) = L;
word1(&a) = 0;
#ifndef Sudden_Underflow
}
}
else {
L = -L >> Exp_shift;
if (L < Exp_shift) {
word0(a) = 0x80000 >> L;
word1(a) = 0;
}
else {
word0(a) = 0;
L -= Exp_shift;
word1(a) = L >= 31 ? 1 : 1 << (31 - L);
}
word0(&a) = 0x80000 >> L;
word1(&a) = 0;
}
else {
word0(&a) = 0;
L -= Exp_shift;
word1(&a) = L >= 31 ? 1 : 1 << (31 - L);
}
#endif
return a;
}
#endif
return dval(&a);
}