Clean up gdtoa code

This commit is contained in:
Justine Tunney 2021-08-19 06:07:37 -07:00
parent 7341336b1a
commit da45c7c80b
48 changed files with 2808 additions and 4721 deletions

View file

@ -1,75 +1,67 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
The author of this software is David M. Gay.
Please send bug reports to David M. Gay <dmg@acm.org>
or Justine Tunney <jtunney@gmail.com>
Copyright (C) 1998, 1999 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
*/
#include "third_party/gdtoa/gdtoa.internal.h"
/* clang-format off */
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998, 1999 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
static Bigint *
bitstob(ULong *bits, int nbits, int *bbits MTd)
bitstob(ULong *bits, int nbits, int *bbits)
{
int i, k;
Bigint *b;
ULong *be, *x, *x0;
i = ULbits;
k = 0;
while(i < nbits) {
i <<= 1;
k++;
}
#ifndef Pack_32
if (!k)
k = 1;
#endif
b = Balloc(k MTa);
}
b = Balloc(k);
be = bits + ((nbits - 1) >> kshift);
x = x0 = b->x;
do {
*x++ = *bits & ALL_ON;
#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:
ret:
return b;
}
}
/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
*
@ -106,7 +98,7 @@ bitstob(ULong *bits, int nbits, int *bbits MTd)
*/
char *
gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve)
gdtoa(const 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
@ -143,9 +135,6 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
to hold the suppressed trailing zeros.
*/
#ifdef MULTIPLE_THREADS
ThInfo *TI = 0;
#endif
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;
@ -154,53 +143,40 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
double d2, ds;
char *s, *s0;
U d, eps;
#ifndef MULTIPLE_THREADS
if (dtoa_result) {
freedtoa(dtoa_result);
dtoa_result = 0;
}
#endif
inex = 0;
kind = *kindp &= ~STRTOG_Inexact;
switch(kind & STRTOG_Retmask) {
case STRTOG_Zero:
case STRTOG_Zero:
goto ret_zero;
case STRTOG_Normal:
case STRTOG_Denormal:
case STRTOG_Normal:
case STRTOG_Denormal:
break;
case STRTOG_Infinite:
case STRTOG_Infinite:
*decpt = -32768;
return nrv_alloc("Infinity", rve, 8 MTb);
case STRTOG_NaN:
return nrv_alloc("Infinity", rve, 8);
case STRTOG_NaN:
*decpt = -32768;
return nrv_alloc("NaN", rve, 3 MTb);
default:
return nrv_alloc("NaN", rve, 3);
default:
return 0;
}
b = bitstob(bits, nbits = fpi->nbits, &bbits MTb);
}
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 MTb);
ret_zero:
Bfree(b);
ret_zero:
*decpt = 1;
return nrv_alloc("0", rve, 1 MTb);
}
return nrv_alloc("0", rve, 1);
}
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
/* 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))
@ -222,25 +198,16 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
* (We could get a more accurate k by invoking log10,
* but this is probably not worthwhile.)
*/
ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
/* correct assumption about exponent range */
if ((j = i) < 0)
j = -j;
if ((j -= 1077) > 0)
ds += j * 7e-17;
k = (int)ds;
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
// TODO: word0(&d) += (be + bbits - 1) << Exp_shift;
// error: third_party/gdtoa/gdtoa.c:244: left shift of negative value -6 'int' 20 'int'
// 4161d8: __die at libc/log/die.c:33
@ -258,69 +225,67 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
// 401d30: cosmo at libc/runtime/cosmo.S:65
// 401173: _start at libc/crt/crt.S:67
word0(&d) += (unsigned)(be + bbits - 1) << Exp_shift;
#endif
if (k >= 0 && k <= Ten_pmax) {
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;
}
}
else if (i >= -4 - Emin || i < Emin)
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:
i = (int)(nbits * .30103) + 3;
ndigits = 0;
break;
case 2:
leftright = 0;
/* no break */
case 4:
if (ndigits <= 0)
ndigits = 1;
ilim = ilim1 = i = ndigits;
break;
case 3:
leftright = 0;
/* no break */
case 5:
i = ndigits + k + 1;
ilim = i;
ilim1 = i - 1;
if (i <= 0)
i = 1;
}
s = s0 = rv_alloc(i MTb);
case 0:
case 1:
i = (int)(nbits * .30103) + 3;
ndigits = 0;
break;
case 2:
leftright = 0;
/* no break */
case 4:
if (ndigits <= 0)
ndigits = 1;
ilim = ilim1 = i = ndigits;
break;
case 3:
leftright = 0;
/* no break */
case 5:
i = ndigits + k + 1;
ilim = i;
ilim1 = i - 1;
if (i <= 0)
i = 1;
}
s = s0 = rv_alloc(i);
if (mode <= 1)
rdir = 0;
else if ( (rdir = fpi->rounding - 1) !=0) {
@ -328,18 +293,10 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
rdir = 2;
if (kind & STRTOG_Neg)
rdir = 3 - rdir;
}
}
/* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir
#ifndef IMPRECISE_INEXACT
&& k == 0
#endif
) {
if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir && k == 0) {
/* Try to get by with floating-point arithmetic. */
i = 0;
d2 = dval(&d);
k0 = k;
@ -353,13 +310,13 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
j &= Bletch - 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) {
@ -368,9 +325,9 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
if (j & 1) {
ieps++;
dval(&d) *= bigtens[i];
}
}
}
}
}
if (k_check && dval(&d) < 1. && ilim > 0) {
if (ilim1 <= 0)
goto fast_failed;
@ -378,7 +335,7 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
k--;
dval(&d) *= 10.;
ieps++;
}
}
dval(&eps) = ieps*dval(&d) + 7.;
word0(&eps) -= (P-1)*Exp_msk1;
if (ilim == 0) {
@ -389,8 +346,7 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
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.
@ -404,17 +360,16 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
if (dval(&d))
inex = STRTOG_Inexlo;
goto ret1;
}
}
if (ds - dval(&d) < dval(&eps))
goto bump_up;
if (++i >= ilim)
break;
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.) {
@ -429,22 +384,18 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
if (dval(&d))
inex = STRTOG_Inexlo;
goto ret1;
}
break;
}
break;
}
#ifndef No_leftright
}
#endif
fast_failed:
}
fast_failed:
s = s0;
dval(&d) = d2;
k = k0;
ilim = ilim0;
}
}
/* Do we have a "small" integer? */
if (be >= 0 && k <= fpi->int_max) {
/* Yes. */
ds = tens[k];
@ -453,17 +404,15 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
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;
#ifdef Check_FLT_ROUNDS
/* If FLT_ROUNDS == 2, L will usually be high by 1 */
if (dval(&d) < 0) {
L--;
dval(&d) += ds;
}
#endif
}
*s++ = '0' + (int)L;
if (dval(&d) == 0.)
break;
@ -473,32 +422,27 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
goto bump_up;
inex = STRTOG_Inexlo;
goto ret1;
}
}
dval(&d) += dval(&d);
#ifdef ROUND_BIASED
if (dval(&d) >= ds)
#else
if (dval(&d) > ds || (dval(&d) == ds && L & 1))
#endif
{
bump_up:
{
bump_up:
inex = STRTOG_Inexhi;
while(*--s == '9')
if (s == s0) {
k++;
*s = '0';
break;
}
}
++*s++;
}
}
else
inex = STRTOG_Inexlo;
break;
}
}
goto ret1;
}
goto ret1;
}
m2 = b2;
m5 = b5;
mhi = mlo = 0;
@ -509,9 +453,9 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
i = be - fpi->emin + 1;
if (mode >= 2 && ilim > 0 && ilim < i)
goto small_ilim;
}
}
else if (mode >= 2) {
small_ilim:
small_ilim:
j = ilim - 1;
if (m5 >= j)
m5 -= j;
@ -519,42 +463,40 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
s5 += j -= m5;
b5 += j;
m5 = 0;
}
}
if ((i = ilim) < 0) {
m2 -= i;
i = 0;
}
}
}
b2 += i;
s2 += i;
mhi = i2b(1 MTb);
}
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) {
mhi = pow5mult(mhi, m5 MTb);
b1 = mult(mhi, b MTb);
Bfree(b MTb);
mhi = pow5mult(mhi, m5);
b1 = mult(mhi, b);
Bfree(b);
b = b1;
}
if ( (j = b5 - m5) !=0)
b = pow5mult(b, j MTb);
}
else
b = pow5mult(b, b5 MTb);
if ( (j = b5 - m5) !=0)
b = pow5mult(b, j);
}
S = i2b(1 MTb);
else
b = pow5mult(b, b5);
}
S = i2b(1);
if (s5 > 0)
S = pow5mult(S, s5 MTb);
S = pow5mult(S, s5);
/* Check for special case that d is a normalized power of 2. */
spec_case = 0;
if (mode < 2) {
if (bbits == 1 && be0 > fpi->emin + 1) {
@ -562,9 +504,8 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
b2++;
s2++;
spec_case = 1;
}
}
}
/* Arrange for convenient computation of quotients:
* shift left if necessary so divisor has 4 leading 0 bits.
*
@ -575,190 +516,170 @@ gdtoa(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
m2 += i;
if ((b2 += i) > 0)
b = lshift(b, b2 MTb);
b = lshift(b, b2);
if ((s2 += i) > 0)
S = lshift(S, s2 MTb);
S = lshift(S, s2);
if (k_check) {
if (cmp(b,S) < 0) {
k--;
b = multadd(b, 10, 0 MTb); /* we botched the k estimate */
b = multadd(b, 10, 0); /* we botched the k estimate */
if (leftright)
mhi = multadd(mhi, 10, 0 MTb);
mhi = multadd(mhi, 10, 0);
ilim = ilim1;
}
}
}
if (ilim <= 0 && mode > 2) {
if (ilim < 0 || cmp(b,S = multadd(S,5,0 MTb)) <= 0) {
if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
/* no digits, fcvt style */
no_digits:
no_digits:
k = -1 - ndigits;
inex = STRTOG_Inexlo;
goto ret;
}
one_digit:
}
one_digit:
inex = STRTOG_Inexhi;
*s++ = '1';
k++;
goto ret;
}
}
if (leftright) {
if (m2 > 0)
mhi = lshift(mhi, m2 MTb);
mhi = lshift(mhi, m2);
/* Compute mlo -- check for special case
* that d is a normalized power of 2.
*/
mlo = mhi;
if (spec_case) {
mhi = Balloc(mhi->k MTb);
mhi = Balloc(mhi->k);
Bcopy(mhi, mlo);
mhi = lshift(mhi, 1 MTb);
}
mhi = lshift(mhi, 1);
}
for(i = 1;;i++) {
dig = quorem(b,S) + '0';
/* Do we yet have the shortest decimal string
* that will round to d?
*/
j = cmp(b, mlo);
delta = diff(S, mhi MTb);
delta = diff(S, mhi);
j1 = delta->sign ? 1 : cmp(b, delta);
Bfree(delta MTb);
#ifndef ROUND_BIASED
Bfree(delta);
if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
if (dig == '9')
goto round_9_up;
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
&& !(bits[0] & 1)
#endif
)) {
}
if (j < 0 || (j == 0 && !mode && !(bits[0] & 1))) {
if (rdir && (b->wds > 1 || b->x[0])) {
if (rdir == 2) {
inex = STRTOG_Inexlo;
goto accept;
}
}
while (cmp(S,mhi) > 0) {
*s++ = dig;
mhi1 = multadd(mhi, 10, 0 MTb);
mhi1 = multadd(mhi, 10, 0);
if (mlo == mhi)
mlo = mhi1;
mhi = mhi1;
b = multadd(b, 10, 0 MTb);
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 MTb);
b = lshift(b, 1);
j1 = cmp(b, S);
#ifdef ROUND_BIASED
if (j1 >= 0 /*)*/
#else
if ((j1 > 0 || (j1 == 0 && dig & 1))
#endif
&& dig++ == '9')
if ((j1 > 0 || (j1 == 0 && dig & 1)) && dig++ == '9')
goto round_9_up;
inex = STRTOG_Inexhi;
}
}
if (b->wds > 1 || b->x[0])
inex = STRTOG_Inexlo;
accept:
accept:
*s++ = dig;
goto ret;
}
}
if (j1 > 0 && rdir != 2) {
if (dig == '9') { /* possible if i == 1 */
round_9_up:
round_9_up:
*s++ = '9';
inex = STRTOG_Inexhi;
goto roundoff;
}
}
inex = STRTOG_Inexhi;
*s++ = dig + 1;
goto ret;
}
}
*s++ = dig;
if (i == ilim)
break;
b = multadd(b, 10, 0 MTb);
b = multadd(b, 10, 0);
if (mlo == mhi)
mlo = mhi = multadd(mhi, 10, 0 MTb);
mlo = mhi = multadd(mhi, 10, 0);
else {
mlo = multadd(mlo, 10, 0 MTb);
mhi = multadd(mhi, 10, 0 MTb);
}
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 MTb);
}
b = multadd(b, 10, 0);
}
/* Round off last digit */
if (rdir) {
if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
goto chopzeros;
goto roundoff;
}
b = lshift(b, 1 MTb);
}
b = lshift(b, 1);
j = cmp(b, S);
#ifdef ROUND_BIASED
if (j >= 0)
#else
if (j > 0 || (j == 0 && dig & 1))
#endif
{
roundoff:
{
roundoff:
inex = STRTOG_Inexhi;
while(*--s == '9')
if (s == s0) {
k++;
*s++ = '1';
goto ret;
}
}
++*s++;
}
}
else {
chopzeros:
chopzeros:
if (b->wds > 1 || b->x[0])
inex = STRTOG_Inexlo;
}
ret:
Bfree(S MTb);
}
ret:
Bfree(S);
if (mhi) {
if (mlo && mlo != mhi)
Bfree(mlo MTb);
Bfree(mhi MTb);
}
ret1:
Bfree(mlo);
Bfree(mhi);
}
ret1:
while(s > s0 && s[-1] == '0')
--s;
Bfree(b MTb);
Bfree(b);
*s = 0;
*decpt = k + 1;
if (rve)
*rve = s;
*kindp |= inex;
return s0;
}
}