mirrors/cosmopolitan
1
0
Fork 0
mirror of https://github.com/jart/cosmopolitan.git synced 2025-01-31 19:43:32 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions
cosmopolitan/third_party/gdtoa/strtod.c
Justine Tunney 226aaf3547 Improve memory safety 
This commit makes numerous refinements to cosmopolitan memory handling.

The default stack size has been reduced from 2mb to 128kb. A new macro
is now provided so you can easily reconfigure the stack size to be any
value you want. Work around the breaking change by adding to your main:

    STATIC_STACK_SIZE(0x00200000);  // 2mb stack

If you're not sure how much stack you need, then you can use:

    STATIC_YOINK("stack_usage_logging");

After which you can `sort -nr o/$MODE/stack.log`. Based on the unit test
suite, nothing in the Cosmopolitan repository (except for Python) needs
a stack size greater than 30kb. There are also new macros for detecting
the size and address of the stack at runtime, e.g. GetStackAddr(). We
also now support sigaltstack() so if you want to see nice looking crash
reports whenever a stack overflow happens, you can put this in main():

    ShowCrashReports();

Under `make MODE=dbg` and `make MODE=asan` the unit testing framework
will now automatically print backtraces of memory allocations when
things like memory leaks happen. Bugs are now fixed in ASAN global
variable overrun detection. The memtrack and asan runtimes also handle
edge cases now. The new tools helped to identify a few memory leaks,
which are fixed by this change.

This change should fix an issue reported in #288 with ARG_MAX limits.
Fixing this doubled the performance of MKDEPS.COM and AR.COM yet again.
2021-10-13 17:27:13 -07:00

686 lines
17 KiB
C
Raw Blame History

/*-*- 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 "libc/errno.h"
#include "libc/runtime/fenv.h"
#include "third_party/gdtoa/gdtoa.internal.h"
/* clang-format off */
#define Avoid_Underflow
#define dplen 1
/* The factor of 2^106 in tiny__gdtoa_tens[4] helps us avoid setting the underflow */
/* flag unnecessarily. It leads to a song and dance at the end of strtod. */
static const double tiny__gdtoa_tens[] = {
1e-16, 1e-32, 1e-64, 1e-128,
9007199254740992.*9007199254740992.e-256
};
static double
s__gdtoa_ulp(U *x, int scale)
{
U u;
int i;
double rv;
rv = __gdtoa_ulp(x);
if (!scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
return rv; /* Is there an example where i <= 0 ? */
word0(&u) = Exp_1 + (i << Exp_shift);
word1(&u) = 0;
return rv * u.d;
}
double
strtod(const char *s00, char **se)
{
int scale;
int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
const char *s, *s0, *s1;
double aadj;
Long L;
U adj, aadj1, rv, rv0;
ULong y, z;
Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
ULong Lsb, Lsb1;
int Rounding;
Rounding = FLT_ROUNDS;
sign = nz0 = nz = decpt = 0;
dval(&rv) = 0.;
for(s = s00;;s++) switch(*s) {
case '-':
sign = 1;
/* no break */
case '+':
if (*++s)
goto break2;
/* no break */
case 0:
goto ret0;
case '\t':
case '\n':
case '\v':
case '\f':
case '\r':
case ' ':
continue;
default:
goto break2;
}
break2:
if (*s == '0') {
{
static const FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI, 0 /*unused*/ };
Long exp;
ULong bits[2];
switch(s[1]) {
case 'x':
case 'X':
{
FPI fpi1 = fpi;
fpi1.rounding = Rounding;
switch((i = __gdtoa_gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
case STRTOG_NoNumber:
s = s00;
sign = 0;
case STRTOG_Zero:
break;
default:
if (bb) {
__gdtoa_copybits(bits, fpi.nbits, bb);
__gdtoa_Bfree(bb);
}
__gdtoa_ULtod(((U*)&rv)->L, bits, exp, i);
}}
goto ret;
}
}
nz0 = 1;
while(*++s == '0') ;
if (!*s)
goto ret;
}
s0 = s;
y = z = 0;
for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
if (nd < 9)
y = 10*y + c - '0';
else if (nd < DBL_DIG + 2)
z = 10*z + c - '0';
nd0 = nd;
if (c == '.') {
c = *++s;
decpt = 1;
if (!nd) {
for(; c == '0'; c = *++s)
nz++;
if (c > '0' && c <= '9') {
s0 = s;
nf += nz;
nz = 0;
goto have_dig;
}
goto dig_done;
}
for(; c >= '0' && c <= '9'; c = *++s) {
have_dig:
nz++;
if (c -= '0') {
nf += nz;
for(i = 1; i < nz; i++)
if (nd++ < 9)
y *= 10;
else if (nd <= DBL_DIG + 2)
z *= 10;
if (nd++ < 9)
y = 10*y + c;
else if (nd <= DBL_DIG + 2)
z = 10*z + c;
nz = 0;
}
}
}/*}*/
dig_done:
e = 0;
if (c == 'e' || c == 'E') {
if (!nd && !nz && !nz0) {
goto ret0;
}
s00 = s;
esign = 0;
switch(c = *++s) {
case '-':
esign = 1;
case '+':
c = *++s;
}
if (c >= '0' && c <= '9') {
while(c == '0')
c = *++s;
if (c > '0' && c <= '9') {
L = c - '0';
s1 = s;
while((c = *++s) >= '0' && c <= '9')
L = 10*L + c - '0';
if (s - s1 > 8 || L > 19999)
/* Avoid confusion from exponents
* so large that e might overflow.
*/
e = 19999; /* safe for 16 bit ints */
else
e = (int)L;
if (esign)
e = -e;
}
else
e = 0;
}
else
s = s00;
}
if (!nd) {
if (!nz && !nz0) {
/* Check for Nan and Infinity */
ULong bits[2];
static FPI fpinan = /* only 52 explicit bits */
{ 52, 1-1023-53+1, 2046-1023-53+1, 1, SI, 0 /*unused*/ };
if (!decpt)
switch(c) {
case 'i':
case 'I':
if (__gdtoa_match(&s,"nf")) {
--s;
if (!__gdtoa_match(&s,"inity"))
++s;
word0(&rv) = 0x7ff00000;
word1(&rv) = 0;
goto ret;
}
break;
case 'n':
case 'N':
if (__gdtoa_match(&s, "an")) {
if (*s == '(' /*)*/
&& __gdtoa_hexnan(&s, &fpinan, bits)
== STRTOG_NaNbits) {
word0(&rv) = 0x7ff00000 | bits[1];
word1(&rv) = bits[0];
}
else {
word0(&rv) = NAN_WORD0;
word1(&rv) = NAN_WORD1;
}
goto ret;
}
}
ret0:
s = s00;
sign = 0;
}
goto ret;
}
e1 = e -= nf;
/* Now we have nd0 digits, starting at s0, followed by a
* decimal point, followed by nd-nd0 digits. The number we're
* after is the integer represented by those digits times
* 10**e */
if (!nd0)
nd0 = nd;
k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;
dval(&rv) = y;
if (k > 9) {
dval(&rv) = __gdtoa_tens[k - 9] * dval(&rv) + z;
}
bd0 = 0;
if (nd <= DBL_DIG) {
if (!e)
goto ret;
if (e > 0) {
if (e <= Ten_pmax) {
/* round correctly FLT_ROUNDS = 2 or 3 */
if (sign) {
rv.d = -rv.d;
sign = 0;
}
/* rv = */ rounded_product(dval(&rv), __gdtoa_tens[e]);
goto ret;
}
i = DBL_DIG - nd;
if (e <= Ten_pmax + i) {
/* A fancier test would sometimes let us do
* this for larger i values.
*/
/* round correctly FLT_ROUNDS = 2 or 3 */
if (sign) {
rv.d = -rv.d;
sign = 0;
}
e -= i;
dval(&rv) *= __gdtoa_tens[i];
/* rv = */ rounded_product(dval(&rv), __gdtoa_tens[e]);
goto ret;
}
}
else if (e >= -Ten_pmax) {
/* round correctly FLT_ROUNDS = 2 or 3 */
if (sign) {
rv.d = -rv.d;
sign = 0;
}
/* rv = */ rounded_quotient(dval(&rv), __gdtoa_tens[-e]);
goto ret;
}
}
e1 += nd - k;
scale = 0;
if (Rounding >= 2) {
if (sign)
Rounding = Rounding == 2 ? 0 : 2;
else
if (Rounding != 2)
Rounding = 0;
}
/* Get starting approximation = rv * 10**e1 */
if (e1 > 0) {
if ( (i = e1 & 15) !=0)
dval(&rv) *= __gdtoa_tens[i];
if (e1 &= ~15) {
if (e1 > DBL_MAX_10_EXP) {
ovfl:
/* Can't trust HUGE_VAL */
switch(Rounding) {
case 0: /* toward 0 */
case 3: /* toward -infinity */
word0(&rv) = Big0;
word1(&rv) = Big1;
break;
default:
word0(&rv) = Exp_mask;
word1(&rv) = 0;
}
range_err:
if (bd0) {
__gdtoa_Bfree(bb);
__gdtoa_Bfree(bd);
__gdtoa_Bfree(bs);
__gdtoa_Bfree(bd0);
__gdtoa_Bfree(delta);
}
errno = ERANGE;
goto ret;
}
e1 >>= 4;
for(j = 0; e1 > 1; j++, e1 >>= 1)
if (e1 & 1)
dval(&rv) *= __gdtoa_bigtens[j];
/* The last __gdtoa_multiplication could overflow. */
word0(&rv) -= P*Exp_msk1;
dval(&rv) *= __gdtoa_bigtens[j];
if ((z = word0(&rv) & Exp_mask)
> Exp_msk1*(DBL_MAX_EXP+Bias-P))
goto ovfl;
if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
/* set to largest number */
/* (Can't trust DBL_MAX) */
word0(&rv) = Big0;
word1(&rv) = Big1;
}
else
word0(&rv) += P*Exp_msk1;
}
}
else if (e1 < 0) {
e1 = -e1;
if ( (i = e1 & 15) !=0)
dval(&rv) /= __gdtoa_tens[i];
if (e1 >>= 4) {
if (e1 >= 1 << n___gdtoa_bigtens)
goto undfl;
if (e1 & Scale_Bit)
scale = 2*P;
for(j = 0; e1 > 0; j++, e1 >>= 1)
if (e1 & 1)
dval(&rv) *= tiny__gdtoa_tens[j];
if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
>> Exp_shift)) > 0) {
/* scaled rv is denormal; zap j low bits */
if (j >= 32) {
word1(&rv) = 0;
if (j >= 53)
word0(&rv) = (P+2)*Exp_msk1;
else
word0(&rv) &= 0xffffffff << (j-32);
}
else
word1(&rv) &= 0xffffffff << j;
}
if (!dval(&rv)) {
undfl:
dval(&rv) = 0.;
if (Rounding == 2)
word1(&rv) = 1;
goto range_err;
}
}
}
/* Now the hard part -- adjusting rv to the correct value.*/
/* Put digits into bd: true value = bd * 10^e */
bd0 = __gdtoa_s2b(s0, nd0, nd, y, dplen);
for(;;) {
bd = __gdtoa_Balloc(bd0->k);
Bcopy(bd, bd0);
bb = __gdtoa_d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
bs = __gdtoa_i2b(1);
if (e >= 0) {
bb2 = bb5 = 0;
bd2 = bd5 = e;
}
else {
bb2 = bb5 = -e;
bd2 = bd5 = 0;
}
if (bbe >= 0)
bb2 += bbe;
else
bd2 -= bbe;
bs2 = bb2;
if (Rounding != 1)
bs2++;
Lsb = LSB;
Lsb1 = 0;
j = bbe - scale;
i = j + bbbits - 1; /* logb(rv) */
j = P + 1 - bbbits;
if (i < Emin) { /* denormal */
i = Emin - i;
j -= i;
if (i < 32)
Lsb <<= i;
else
Lsb1 = Lsb << (i-32);
}
bb2 += j;
bd2 += j;
bd2 += scale;
i = bb2 < bd2 ? bb2 : bd2;
if (i > bs2)
i = bs2;
if (i > 0) {
bb2 -= i;
bd2 -= i;
bs2 -= i;
}
if (bb5 > 0) {
bs = __gdtoa_pow5mult(bs, bb5);
bb1 = __gdtoa_mult(bs, bb);
__gdtoa_Bfree(bb);
bb = bb1;
}
if (bb2 > 0)
bb = __gdtoa_lshift(bb, bb2);
if (bd5 > 0)
bd = __gdtoa_pow5mult(bd, bd5);
if (bd2 > 0)
bd = __gdtoa_lshift(bd, bd2);
if (bs2 > 0)
bs = __gdtoa_lshift(bs, bs2);
delta = __gdtoa_diff(bb, bd);
dsign = delta->sign;
delta->sign = 0;
i = __gdtoa_cmp(delta, bs);
if (Rounding != 1) {
if (i < 0) {
/* Error is less than an __gdtoa_ulp */
if (!delta->x[0] && delta->wds <= 1) {
/* exact */
break;
}
if (Rounding) {
if (dsign) {
dval(&adj) = 1.;
goto apply_adj;
}
}
else if (!dsign) {
dval(&adj) = -1.;
if (!word1(&rv)
&& !(word0(&rv) & Frac_mask)) {
y = word0(&rv) & Exp_mask;
if (!scale || y > 2*P*Exp_msk1)
{
delta = __gdtoa_lshift(delta,Log2P);
if (__gdtoa_cmp(delta, bs) <= 0)
dval(&adj) = -0.5;
}
}
apply_adj:
if (scale && (y = word0(&rv) & Exp_mask)
<= 2*P*Exp_msk1)
word0(&adj) += (2*P+1)*Exp_msk1 - y;
dval(&rv) += adj.d*__gdtoa_ulp(&rv);
}
break;
}
dval(&adj) = __gdtoa_ratio(delta, bs);
if (adj.d < 1.)
dval(&adj) = 1.;
if (adj.d <= 0x7ffffffe) {
/* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */
y = adj.d;
if (y != adj.d) {
if (!((Rounding>>1) ^ dsign))
y++;
dval(&adj) = y;
}
}
if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
word0(&adj) += (2*P+1)*Exp_msk1 - y;
dval(&adj) *= __gdtoa_ulp(&rv); /* XXX */
if (dsign) {
if (word0(&rv) == Big0 && word1(&rv) == Big1)
goto ovfl;
dval(&rv) += adj.d;
}
else
dval(&rv) -= adj.d;
goto cont;
}
if (i < 0) {
/* Error is less than half an __gdtoa_ulp -- check for
* special case of mantissa a power of two.
*/
if (dsign || word1(&rv) || word0(&rv) & Bndry_mask ||
(word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1) {
break;
}
if (!delta->x[0] && delta->wds <= 1) {
/* exact result */
break;
}
delta = __gdtoa_lshift(delta,Log2P);
if (__gdtoa_cmp(delta, bs) > 0)
goto drop_down;
break;
}
if (i == 0) {
/* exactly half-way between */
if (dsign) {
if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
&& word1(&rv) == (
(scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
0xffffffff)) {
/*boundary case -- __gdtoa_increment exponent*/
if (word0(&rv) == Big0 && word1(&rv) == Big1)
goto ovfl;
word0(&rv) = (word0(&rv) & Exp_mask)
+ Exp_msk1
;
word1(&rv) = 0;
dsign = 0;
break;
}
}
else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
drop_down:
/* boundary case -- __gdtoa_decrement exponent */
if (scale) {
L = word0(&rv) & Exp_mask;
if (L <= (2*P+1)*Exp_msk1) {
if (L > (P+2)*Exp_msk1)
/* round even ==> */
/* accept rv */
break;
/* rv = smallest denormal */
goto undfl;
}
}
L = (word0(&rv) & Exp_mask) - Exp_msk1;
word0(&rv) = L | Bndry_mask1;
word1(&rv) = 0xffffffff;
break;
}
if (Lsb1) {
if (!(word0(&rv) & Lsb1))
break;
}
else if (!(word1(&rv) & Lsb))
break;
if (dsign)
dval(&rv) += s__gdtoa_ulp(&rv, scale);
else {
dval(&rv) -= s__gdtoa_ulp(&rv, scale);
if (!dval(&rv))
goto undfl;
}
dsign = 1 - dsign;
break;
}
if ((aadj = __gdtoa_ratio(delta, bs)) <= 2.) {
if (dsign)
aadj = dval(&aadj1) = 1.;
else if (word1(&rv) || word0(&rv) & Bndry_mask) {
if (word1(&rv) == Tiny1 && !word0(&rv))
goto undfl;
aadj = 1.;
dval(&aadj1) = -1.;
}
else {
/* special case -- power of FLT_RADIX to be */
/* rounded down... */
if (aadj < 2./FLT_RADIX)
aadj = 1./FLT_RADIX;
else
aadj *= 0.5;
dval(&aadj1) = -aadj;
}
}
else {
aadj *= 0.5;
dval(&aadj1) = dsign ? aadj : -aadj;
switch(Rounding) {
case 2: /* towards +infinity */
dval(&aadj1) -= 0.5;
break;
case 0: /* towards 0 */
case 3: /* towards -infinity */
dval(&aadj1) += 0.5;
}
}
y = word0(&rv) & Exp_mask;
/* Check for overflow */
if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
dval(&rv0) = dval(&rv);
word0(&rv) -= P*Exp_msk1;
dval(&adj) = dval(&aadj1) * __gdtoa_ulp(&rv);
dval(&rv) += dval(&adj);
if ((word0(&rv) & Exp_mask) >=
Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
goto ovfl;
word0(&rv) = Big0;
word1(&rv) = Big1;
goto cont;
}
else
word0(&rv) += P*Exp_msk1;
}
else {
if (scale && y <= 2*P*Exp_msk1) {
if (aadj <= 0x7fffffff) {
if ((z = aadj) <= 0)
z = 1;
aadj = z;
dval(&aadj1) = dsign ? aadj : -aadj;
}
word0(&aadj1) += (2*P+1)*Exp_msk1 - y;
}
dval(&adj) = dval(&aadj1) * __gdtoa_ulp(&rv);
dval(&rv) += dval(&adj);
}
z = word0(&rv) & Exp_mask;
if (!scale)
if (y == z) {
/* Can we stop now? */
L = (Long)aadj;
aadj -= L;
/* The tolerances below are conservative. */
if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
if (aadj < .4999999 || aadj > .5000001)
break;
}
else if (aadj < .4999999/FLT_RADIX)
break;
}
cont:
__gdtoa_Bfree(bb);
__gdtoa_Bfree(bd);
__gdtoa_Bfree(bs);
__gdtoa_Bfree(delta);
}
__gdtoa_Bfree(bb);
__gdtoa_Bfree(bd);
__gdtoa_Bfree(bs);
__gdtoa_Bfree(bd0);
__gdtoa_Bfree(delta);
if (scale) {
word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
word1(&rv0) = 0;
dval(&rv) *= dval(&rv0);
/* try to avoid the bug of testing an 8087 register value */
if (!(word0(&rv) & Exp_mask))
errno = ERANGE;
}
ret:
if (se)
*se = (char *)s;
return sign ? -dval(&rv) : dval(&rv);
}
Reference in a new issue View git blame Copy permalink