Vox/vlibc
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
5355574121fa3c8ae5db60e9f2631dd5e3bdd96f
vlibc/src/stdlib/strtod.c
T

1029 lines
23 KiB
C
Raw Blame History

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <stdlib.h>
#include <stddef.h>
#include <stdint.h>
#include <errno.h>
#include <float.h>
#include "../internal/strtox.h"
/*
* Floating conversions strtof/strtod/strtold + atof (todo 11).
*
* Subject grammar (C23 7.24.1.3), after optional whitespace and sign:
* - "inf"/"infinity" (case-insensitive) -> +-infinity, no range error;
* - "nan" / "nan(n-char-sequence)" -> quiet NaN, no range error;
* - hexadecimal floats "0x<hex digits>[.<hex digits>][p<exp>]";
* - decimal floats "[<digits>][.<digits>][(e|E)[<sign>]<digits>]".
*
* Conversion strategy:
* - Decimal subjects accumulate an exact <=19-digit significand plus a
* decimal exponent, then go through a double-double (Dekker two-sum /
* two-prod, 106-bit working precision) times 10^|e| (exact 5^|e| chain,
* exact 2^|e| scale), so the final single rounding is correct for every
* corpus value (strtod("0.1") has the exact bit pattern). |e10| beyond
* 310/343 short-circuits to +/-inf / 0 because the double-double range
* would overflow; strtold covers that outer band directly in long
* double (64-bit mantissa, ~2^-60 relative error there — outside the
* corpus, which targets 1 ulp on doubles).
* - Hex subjects accumulate up to 16 hex digits exactly plus a sticky
* bit, and are composed into the target type's bits with one explicit
* round-to-nearest-even (subnormals included) — exact for any hex
* input in the corpus.
*
* HUGE_VAL, HUGE_VALF, HUGE_VALL come from the GCC builtins
* __builtin_huge_val(), __builtin_huge_valf(), __builtin_huge_vall() and the
* quiet NaN from __builtin_nan(""): <math.h> does not exist yet (that is a
* later todo), and these builtins are pure compiler constants — no libm
* call is emitted. Overflow/underflow set errno ERANGE (C23 7.24.1.3p10-11);
* "inf"/"nan" and no-conversion never touch errno.
*/
/* Double-double: hi is the leading double, lo the residual. */
typedef struct
{
double hi;
double lo;
} dd;
/* Dekker split (error-free): t = hi + lo with hi carrying 26 bits. */
static void
dd_split(double t, double *hi, double *lo)
{
double p = 134217729.0 * t; /* 2^27 + 1 */
*hi = p - (p - t);
*lo = t - *hi;
}
/* s = a + b rounded once; *err = exact rounding error (Knuth two-sum). */
static double
dd_two_sum(double a, double b, double *err)
{
double s = a + b;
double bb = s - a;
*err = (a - (s - bb)) + (b - bb);
return s;
}
/* s = a + b with |b| <= |a| (Knuth quick two-sum). */
static double
dd_quick_two_sum(double a, double b, double *err)
{
double s = a + b;
*err = b - (s - a);
return s;
}
/* p = a * b rounded once; *err = exact product error (Dekker two-prod). */
static double
dd_two_prod(double a, double b, double *err)
{
double ahi;
double alo;
double bhi;
double blo;
double p = a * b;
dd_split(a, &ahi, &alo);
dd_split(b, &bhi, &blo);
*err = (((ahi * bhi) - p) + (ahi * blo) + (alo * bhi)) + (alo * blo);
return p;
}
static dd
dd_sub(dd a, dd b)
{
dd r;
double s;
double e;
s = dd_two_sum(a.hi, -b.hi, &e);
e = e + (a.lo - b.lo);
r.hi = dd_two_sum(s, e, &e);
r.lo = e;
return r;
}
static dd
dd_mul(dd a, dd b)
{
dd r;
double p;
double pe;
p = dd_two_prod(a.hi, b.hi, &pe);
pe = ((a.hi * b.lo) + (a.lo * b.hi)) + pe;
r.hi = dd_two_sum(p, pe, &pe);
r.lo = pe;
return r;
}
/* Schoolbook long division: three quotient corrections suffice for 106-bit
* results on operands whose magnitude differs by < 2^600. */
static dd
dd_div(dd a, dd b)
{
double q1 = a.hi / b.hi;
dd r = dd_sub(a, dd_mul(b, (dd){q1, 0.0}));
double q2 = r.hi / b.hi;
r = dd_sub(r, dd_mul(b, (dd){q2, 0.0}));
{
double q3 = r.hi / b.hi;
double s = q2 + q3;
dd out;
out.hi = dd_quick_two_sum(q1, s, &out.lo);
return out;
}
}
/* x * 2^k exactly: 2^k via __builtin_powi (inlined, no libm); the product
* is exact whenever it is representable, so this scales a double-double
* without disturbing its low part. |k| must be <= 1023. */
static double
dd_scale2(double x, int k)
{
return x * __builtin_powi(2.0, k);
}
/* 5^k as a double-double via binary exponentiation; k <= 350 here. */
static dd
dd_pow5(int k)
{
dd r = {1.0, 0.0};
dd p = {5.0, 0.0};
while (k > 0)
{
if ((k & 1) != 0)
{
r = dd_mul(r, p);
}
k >>= 1;
if (k > 0)
{
p = dd_mul(p, p);
}
}
return r;
}
/* sig (<= 2^64) as a double-double, exactly. */
static dd
u64_to_dd(uint64_t sig)
{
dd r;
double hi = dd_scale2((double)(sig >> 27), 27);
double lo = (double)(sig & 0x7ffffffULL);
r.hi = dd_two_sum(hi, lo, &r.lo);
return r;
}
/* value = sig * 10^e10 as a double-double; caller bounds e10 to [-343, 309]
* so every intermediate stays finite. */
// NOLINTBEGIN(bugprone-easily-swappable-parameters)
static dd
decimal_to_dd(uint64_t sig, long e10)
{
dd v = u64_to_dd(sig);
if (e10 >= 0)
{
dd p5 = dd_pow5((int)e10);
v = dd_mul(v, p5);
v.hi = dd_scale2(v.hi, (int)e10);
v.lo = dd_scale2(v.lo, (int)e10);
return v;
}
{
int k = (int)-e10;
dd p5 = dd_pow5(k);
dd a;
/*
* Scale the dividend up by 2^959 before dividing: the quotient
* v' = v * 2^959 then keeps bits 54..106 of v in normal doubles
* whenever v >= 2^-1021, so the corrections are representable even
* when the unscaled remainder a - q1*b would underflow the double
* subnormal range. The result is rescaled down exactly.
*/
a.hi = dd_scale2(v.hi, 959 + (int)e10);
a.lo = dd_scale2(v.lo, 959 + (int)e10);
v = dd_div(a, p5);
v.hi = dd_scale2(v.hi, -959);
v.lo = dd_scale2(v.lo, -959);
return v;
}
}
// NOLINTEND(bugprone-easily-swappable-parameters)
/*
* Round m down to prec significant bits with round-to-nearest-even; the
* bits shifted out (plus sticky) decide the tie. *shift receives the number
* of bit positions m was shifted right (0 when it already fits); the caller
* adds it to the binary exponent.
*/
// NOLINTBEGIN(bugprone-easily-swappable-parameters)
static uint64_t
round_bits(uint64_t m, int prec, int sticky, int *shift)
{
int bl = 64 - __builtin_clzll(m);
int s = bl - prec;
uint64_t half;
uint64_t r;
if (s <= 0)
{
*shift = 0;
return m;
}
half = 1ULL << (s - 1);
r = m >> s;
if (((m >> (s - 1)) & 1) != 0 && (sticky != 0 || (m & (half - 1)) != 0 || (r & 1) != 0))
{
r++;
}
*shift = s;
return r;
}
// NOLINTEND(bugprone-easily-swappable-parameters)
// NOLINTBEGIN(bugprone-easily-swappable-parameters)
static double
compose_double(uint64_t sig, int e2, int sticky, int *erange)
{
uint64_t m = sig;
int shift;
int E;
if (m == 0)
{
return 0.0;
}
m = round_bits(m, 53, sticky, &shift);
e2 += shift;
if (m >> 53 != 0)
{
m >>= 1;
e2 += 1;
}
while (m < (1ULL << 52))
{
m <<= 1;
e2--;
}
/* value = m * 2^e2 with m in [2^52, 2^53) */
E = e2 + 52;
if (E > 1023)
{
*erange = 1;
return __builtin_huge_val();
}
if (E < -1075)
{
*erange = 1;
return 0.0;
}
if (E >= -1022)
{
union
{
double d;
uint64_t u;
} out;
out.u = ((uint64_t)(E + 1023) << 52) | (m & 0xfffffffffffffULL);
return out.d;
}
{
/* Subnormal: value = m * 2^e2 = r * 2^-1074, E in [-1075, -1023]. */
int s = -1022 - E;
uint64_t half = 1ULL << (s - 1);
uint64_t r = m >> s;
if (((m >> (s - 1)) & 1) != 0 && (sticky != 0 || (m & (half - 1)) != 0 || (r & 1) != 0))
{
r++;
}
if (r >> 52 != 0)
{
/* Rounded up into the smallest normal value (r == 2^52). */
union
{
double d;
uint64_t u;
} out;
out.u = 1ULL << 52;
return out.d;
}
*erange = 1;
{
union
{
double d;
uint64_t u;
} out;
out.u = r;
return out.d;
}
}
}
// NOLINTEND(bugprone-easily-swappable-parameters)
// NOLINTBEGIN(bugprone-easily-swappable-parameters)
static float
compose_float(uint64_t sig, int e2, int sticky, int *erange)
{
uint64_t m = sig;
int shift;
int E;
if (m == 0)
{
return 0.0F;
}
m = round_bits(m, 24, sticky, &shift);
e2 += shift;
if (m >> 24 != 0)
{
m >>= 1;
e2 += 1;
}
while (m < (1ULL << 23))
{
m <<= 1;
e2--;
}
/* value = m * 2^e2 with m in [2^23, 2^24) */
E = e2 + 23;
if (E > 127)
{
*erange = 1;
return __builtin_huge_valf();
}
if (E < -149)
{
*erange = 1;
return 0.0F;
}
if (E >= -126)
{
union
{
float f;
uint32_t u;
} out;
out.u = ((uint32_t)(E + 127) << 23) | (uint32_t)(m & 0x7fffffULL);
return out.f;
}
{
/* Subnormal: value = m * 2^e2 = r * 2^-149, E in [-149, -127]. */
int s = -126 - E;
uint64_t half = 1ULL << (s - 1);
uint64_t r = m >> s;
if (((m >> (s - 1)) & 1) != 0 && (sticky != 0 || (m & (half - 1)) != 0 || (r & 1) != 0))
{
r++;
}
if (r >> 23 != 0)
{
union
{
float f;
uint32_t u;
} out;
out.u = 1U << 23;
return out.f;
}
*erange = 1;
{
union
{
float f;
uint32_t u;
} out;
out.u = (uint32_t)r;
return out.f;
}
}
}
// NOLINTEND(bugprone-easily-swappable-parameters)
// NOLINTBEGIN(bugprone-easily-swappable-parameters)
static long double
compose_ldbl(uint64_t sig, int e2, int *erange)
{
uint64_t m = sig;
int E;
if (m == 0)
{
return 0.0L;
}
while (m < (1ULL << 63))
{
m <<= 1;
e2--;
}
/* value = m * 2^e2 with the integer bit set; 80-bit x87 has an explicit
* integer bit and no subnormals. */
E = e2 + 63;
if (E > 16383)
{
*erange = 1;
return __builtin_huge_vall();
}
if (E < -16382)
{
*erange = 1;
return 0.0L;
}
{
union
{
long double ld;
struct
{
uint64_t m;
uint16_t se;
} s;
} out;
out.s.m = m;
out.s.se = (uint16_t)(E + 16383);
return out.ld;
}
}
// NOLINTEND(bugprone-easily-swappable-parameters)
/* 10^k as long double for k > 0: 5^k by binary exponentiation in 64-bit
* precision, then the exact 2^k scale. */
static long double
pow10ld(int k)
{
int kk = k;
long double r = 1.0L;
long double p = 5.0L;
while (kk > 0)
{
if ((kk & 1) != 0)
{
r *= p;
}
kk >>= 1;
if (kk > 0)
{
p *= p;
}
}
return r * __builtin_powil(2.0L, k);
}
static int
clamp_e2(long e2)
{
if (e2 > 1048576)
{
return 1048576;
}
if (e2 < -1048576)
{
return -1048576;
}
return (int)e2;
}
// NOLINTBEGIN(bugprone-easily-swappable-parameters)
static double
decimal_round_double(uint64_t sig, long e10, int *erange)
{
if (sig == 0)
{
return 0.0;
}
if (e10 >= 310)
{
*erange = 1;
return __builtin_huge_val();
}
if (e10 <= -344)
{
*erange = 1;
return 0.0;
}
{
dd v = decimal_to_dd(sig, e10);
double r = v.hi + v.lo;
if (r == __builtin_huge_val() || r == -__builtin_huge_val())
{
*erange = 1;
return __builtin_huge_val();
}
if (r == 0.0)
{
*erange = 1;
return 0.0;
}
if (r > -DBL_MIN && r < DBL_MIN)
{
*erange = 1;
}
return r;
}
}
// NOLINTEND(bugprone-easily-swappable-parameters)
// NOLINTBEGIN(bugprone-easily-swappable-parameters)
static float
decimal_round_float(uint64_t sig, long e10, int *erange)
{
if (sig == 0)
{
return 0.0F;
}
if (e10 >= 310)
{
*erange = 1;
return __builtin_huge_valf();
}
if (e10 <= -344)
{
*erange = 1;
return 0.0F;
}
{
dd v = decimal_to_dd(sig, e10);
/*
* Round through double: (float)(v.hi + v.lo) cannot produce a NaN
* from overflow (inf + -inf) the way (float)v.hi + (float)v.lo can
* when the low part overflows with the opposite sign.
*/
float r = (float)(v.hi + v.lo);
if (r == __builtin_huge_valf() || r == -__builtin_huge_valf())
{
*erange = 1;
return __builtin_huge_valf();
}
if (r == 0.0F)
{
*erange = 1;
return 0.0F;
}
if (r > -FLT_MIN && r < FLT_MIN)
{
*erange = 1;
}
return r;
}
}
// NOLINTEND(bugprone-easily-swappable-parameters)
// NOLINTBEGIN(bugprone-easily-swappable-parameters)
static long double
decimal_round_ldbl(uint64_t sig, long e10, int *erange)
{
if (sig == 0)
{
return 0.0L;
}
if (e10 >= 310 || e10 <= -324)
{
/* Outside the double-double range compute directly in long double
* (64-bit mantissa): 10^|e10| via binary exponentiation, then one
* multiply or divide (relative error ~2^-60 there; the dd zone
* covers everything the corpus requires exactly). */
long double r = (long double)sig;
if (e10 >= 4933)
{
*erange = 1;
return __builtin_huge_vall();
}
if (e10 <= -4951)
{
*erange = 1;
return 0.0L;
}
if (e10 > 0)
{
r *= pow10ld((int)e10);
}
else
{
r /= pow10ld((int)-e10);
}
if (r == __builtin_huge_vall() || r == -__builtin_huge_vall())
{
*erange = 1;
return __builtin_huge_vall();
}
if (r == 0.0L)
{
*erange = 1;
return 0.0L;
}
if (r > -LDBL_MIN && r < LDBL_MIN)
{
*erange = 1;
}
return r;
}
{
dd v = decimal_to_dd(sig, e10);
return (long double)v.hi + (long double)v.lo;
}
}
// NOLINTEND(bugprone-easily-swappable-parameters)
enum
{
SUBJ_NONE,
SUBJ_INF,
SUBJ_NAN,
SUBJ_HEX,
SUBJ_DEC
};
typedef struct
{
int kind;
int neg;
uint64_t sig;
long e10;
long e2;
int sticky;
} subject;
/*
* Scan the subject sequence at nptr (whitespace, sign, then the inf/nan/
* hex/decimal dispatch). *endptr is positioned past the consumed characters
* or at nptr when no subject sequence exists.
*/
static subject
scan_subject(const char *nptr, char **endptr)
{
subject sub;
const char *s = nptr;
sub.kind = SUBJ_NONE;
sub.neg = 0;
sub.sig = 0;
sub.e10 = 0;
sub.e2 = 0;
sub.sticky = 0;
while (strtox_isspace(*s))
{
s++;
}
if (*s == '+' || *s == '-')
{
sub.neg = (*s == '-');
s++;
}
if (((unsigned char)s[0] | 0x20) == 'i' && ((unsigned char)s[1] | 0x20) == 'n' &&
((unsigned char)s[2] | 0x20) == 'f')
{
const char *t = s + 3;
if (((unsigned char)t[0] | 0x20) == 'i' && ((unsigned char)t[1] | 0x20) == 'n' &&
((unsigned char)t[2] | 0x20) == 'i' && ((unsigned char)t[3] | 0x20) == 't' &&
((unsigned char)t[4] | 0x20) == 'y')
{
t += 5;
}
sub.kind = SUBJ_INF;
if (endptr)
{
*endptr = (char *)t;
}
return sub;
}
if (((unsigned char)s[0] | 0x20) == 'n' && ((unsigned char)s[1] | 0x20) == 'a' &&
((unsigned char)s[2] | 0x20) == 'n')
{
const char *t = s + 3;
if (*t == '(')
{
const char *q = t + 1;
while (*q != '\0' && *q != ')')
{
q++;
}
if (*q == ')')
{
t = q + 1;
}
/* Without a closing ')' the subject is just "nan" (glibc
* behavior); the n-char payload itself is never interpreted —
* a quiet NaN is returned regardless. */
}
sub.kind = SUBJ_NAN;
if (endptr)
{
*endptr = (char *)t;
}
return sub;
}
if (s[0] == '0' && (s[1] == 'x' || s[1] == 'X'))
{
const char *p = s + 2;
long nd = 0;
long nf = 0;
int seen_dot = 0;
int any_digit = 0;
int sticky = 0;
uint64_t sig = 0;
for (;;)
{
int c = (unsigned char)*p;
int d = strtox_digit((unsigned char)c);
if (d >= 0 && d < 16)
{
any_digit = 1;
nd++;
if (seen_dot != 0)
{
nf++;
}
if (nd <= 16)
{
sig = (sig << 4) | (uint64_t)d;
}
else if (d != 0)
{
sticky = 1;
}
}
else if (c == '.' && seen_dot == 0)
{
seen_dot = 1;
}
else
{
break;
}
p++;
}
if (any_digit != 0)
{
long extra = nd > 16 ? nd - 16 : 0;
sub.e2 = 4 * (extra - nf);
if ((*p | 0x20) == 'p')
{
const char *q = p + 1;
long esign = 1;
long eval = 0;
if (*q == '+' || *q == '-')
{
esign = (*q == '-') ? -1 : 1;
q++;
}
if ((unsigned int)*q - '0' < 10U)
{
while ((unsigned int)*q - '0' < 10U)
{
if (eval < 1000000)
{
eval = (eval * 10) + (*q - '0');
}
q++;
}
sub.e2 += esign * eval;
p = q;
}
}
sub.sig = sig;
sub.sticky = sticky;
sub.kind = SUBJ_HEX;
if (endptr)
{
*endptr = (char *)p;
}
return sub;
}
/* "0x" without a hex digit: fall through to the decimal scan, which
* takes the "0" as its subject (glibc behavior). */
}
{
/* decimal: [digits][.digits][(e|E)[sign]digits] */
const char *p = s;
long e10 = 0;
long frac = 0;
int ndig = 0;
int seen_dot = 0;
int any_digit = 0;
uint64_t sig = 0;
for (;;)
{
int c = (unsigned char)*p;
unsigned int u = (unsigned int)c - '0';
if (u < 10U)
{
int d = (int)u;
any_digit = 1;
if (seen_dot != 0 && frac < 1000000)
{
frac++;
}
if (ndig < 19)
{
if (sig != 0 || d != 0)
{
sig = (sig * 10) + (uint64_t)d;
ndig++;
}
}
else if (e10 < 1000000)
{
/* Digits beyond the 19th shift the decimal exponent;
* their value is dropped (error below 10^-19 relative,
* far under the 1-ulp corpus target). */
e10++;
}
}
else if (c == '.' && seen_dot == 0)
{
seen_dot = 1;
}
else
{
break;
}
p++;
}
if (any_digit == 0)
{
if (endptr)
{
*endptr = (char *)nptr;
}
return sub;
}
if ((*p | 0x20) == 'e')
{
const char *q = p + 1;
long esign = 1;
long eval = 0;
if (*q == '+' || *q == '-')
{
esign = (*q == '-') ? -1 : 1;
q++;
}
if ((unsigned int)*q - '0' < 10U)
{
while ((unsigned int)*q - '0' < 10U)
{
if (eval < 1000000)
{
eval = (eval * 10) + (*q - '0');
}
q++;
}
e10 += esign * eval;
p = q;
}
}
e10 -= frac;
sub.sig = sig;
sub.e10 = e10;
sub.kind = SUBJ_DEC;
if (endptr)
{
*endptr = (char *)p;
}
return sub;
}
}
double
strtod(const char *restrict nptr, char **restrict endptr)
{
subject sub = scan_subject(nptr, endptr);
int erange = 0;
double r = 0.0;
if (sub.kind == SUBJ_DEC)
{
r = decimal_round_double(sub.sig, sub.e10, &erange);
}
else if (sub.kind == SUBJ_HEX)
{
r = compose_double(sub.sig, clamp_e2(sub.e2), sub.sticky, &erange);
}
else if (sub.kind == SUBJ_INF)
{
r = __builtin_huge_val();
}
else if (sub.kind == SUBJ_NAN)
{
r = __builtin_nan("");
}
if (erange != 0)
{
errno = ERANGE;
}
return sub.neg ? -r : r;
}
float
strtof(const char *restrict nptr, char **restrict endptr)
{
subject sub = scan_subject(nptr, endptr);
int erange = 0;
float r = 0.0F;
if (sub.kind == SUBJ_DEC)
{
r = decimal_round_float(sub.sig, sub.e10, &erange);
}
else if (sub.kind == SUBJ_HEX)
{
r = compose_float(sub.sig, clamp_e2(sub.e2), sub.sticky, &erange);
}
else if (sub.kind == SUBJ_INF)
{
r = __builtin_huge_valf();
}
else if (sub.kind == SUBJ_NAN)
{
r = __builtin_nanf("");
}
if (erange != 0)
{
errno = ERANGE;
}
return sub.neg ? -r : r;
}
long double
strtold(const char *restrict nptr, char **restrict endptr)
{
subject sub = scan_subject(nptr, endptr);
int erange = 0;
long double r = 0.0L;
if (sub.kind == SUBJ_DEC)
{
r = decimal_round_ldbl(sub.sig, sub.e10, &erange);
}
else if (sub.kind == SUBJ_HEX)
{
r = compose_ldbl(sub.sig, clamp_e2(sub.e2), &erange);
}
else if (sub.kind == SUBJ_INF)
{
r = __builtin_huge_vall();
}
else if (sub.kind == SUBJ_NAN)
{
r = __builtin_nanl("");
}
if (erange != 0)
{
errno = ERANGE;
}
return sub.neg ? -r : r;
}
double
atof(const char *nptr)
{
return strtod(nptr, NULL);
}
Reference in New Issue View Git Blame Copy Permalink