diff --git a/Makefile.am b/Makefile.am index e5ce8ef..45b237c 100644 --- a/Makefile.am +++ b/Makefile.am @@ -187,8 +187,10 @@ VLIBC_CORE_SRCS = \ src/mman/mmap.c \ src/passwd/passwd.c \ src/group/group.c \ + src/math/cbrt.c \ src/math/ceil.c \ src/math/copysign.c \ + src/math/exp.c \ src/math/fabs.c \ src/math/fdim.c \ src/math/floor.c \ @@ -196,20 +198,25 @@ VLIBC_CORE_SRCS = \ src/math/fmin.c \ src/math/fmod.c \ src/math/frexp.c \ + src/math/hypot.c \ src/math/ilogb.c \ src/math/ldexp.c \ src/math/llrint.c \ src/math/llround.c \ + src/math/log.c \ + src/math/log1p.c \ src/math/logb.c \ src/math/lrint.c \ src/math/lround.c \ src/math/modf.c \ src/math/nearbyint.c \ + src/math/pow.c \ src/math/remainder.c \ src/math/remquo.c \ src/math/rint.c \ src/math/round.c \ src/math/scalbn.c \ + src/math/sqrt.c \ src/math/trunc.c # Profile-gated sources, EXISTING gates preserved (strlcpy/strlcat stay L2 diff --git a/include/math.h b/include/math.h index 395a776..0403597 100644 --- a/include/math.h +++ b/include/math.h @@ -463,6 +463,170 @@ llround(double x); long long llroundl(long double x); +/* + * e raised to the power x (C23 7.12.6.1). exp(+0) == exp(-0) == 1; + * exp(-Inf) == +0; exp(+Inf) == +Inf; exp(NaN) == NaN. A result too + * large or too small is a range error and sets errno to ERANGE, so no + * const attribute. + */ +double +exp(double x); + +float +expf(float x); + +long double +expl(long double x); + +/* + * 2 raised to the power x (C23 7.12.6.2). exp2(+0) == exp2(-0) == 1; + * exp2(-Inf) == +0; exp2(+Inf) == +Inf; exp2(NaN) == NaN. A result too + * large or too small is a range error and sets errno to ERANGE, so no + * const attribute. + */ +double +exp2(double x); + +float +exp2f(float x); + +long double +exp2l(long double x); + +/* + * e raised to the power x minus 1 (C23 7.12.6.3). expm1(+0) == +0; + * expm1(-0) == -0; expm1(-Inf) == -1; expm1(+Inf) == +Inf; + * expm1(NaN) == NaN. A result too large is a range error and sets + * errno to ERANGE, so no const attribute. + */ +double +expm1(double x); + +float +expm1f(float x); + +long double +expm1l(long double x); + +/* + * Natural logarithm of x (C23 7.12.6.7). log(+0) == log(-0) == -Inf; + * log(1) == 0; log(-x) == NaN (domain error); log(+Inf) == +Inf; + * log(NaN) == NaN. A negative argument is a domain error setting errno + * to EDOM and a zero argument is a range error setting errno to ERANGE, + * so no const attribute. + */ +double +log(double x); + +float +logf(float x); + +long double +logl(long double x); + +/* + * Base-2 logarithm of x (C23 7.12.6.8). log2(1) == 0; log2(2) == 1. + * Special values, domain errors, and range errors as for log, so no + * const attribute. + */ +double +log2(double x); + +float +log2f(float x); + +long double +log2l(long double x); + +/* + * Base-10 logarithm of x (C23 7.12.6.9). log10(1) == 0; log10(10) == 1. + * Special values, domain errors, and range errors as for log, so no + * const attribute. + */ +double +log10(double x); + +float +log10f(float x); + +long double +log10l(long double x); + +/* + * Natural logarithm of 1 + x (C23 7.12.6.13). log1p(+-0) == +-0; + * log1p(-1) == -Inf (a pole error setting errno to ERANGE); log1p(x) for + * x < -1 is a domain error setting errno to EDOM; log1p(+Inf) == +Inf; + * log1p(NaN) == NaN. Accurate for x close to 0, so no const attribute. + */ +double +log1p(double x); + +float +log1pf(float x); + +long double +log1pl(long double x); + +/* + * Principal square root of x (C23 7.12.5.4). sqrt(+-0) == +-0; sqrt(x) + * for x < 0 is a domain error returning NaN and setting errno to EDOM; + * sqrt(+Inf) == +Inf; sqrt(NaN) == NaN. Correctly rounded; no const + * attribute because of the errno path. + */ +double +sqrt(double x); + +float +sqrtf(float x); + +long double +sqrtl(long double x); + +/* + * Cube root of x (C23 7.12.5.5). cbrt(+-0) == +-0; cbrt(+-Inf) == +-Inf; + * cbrt(NaN) == NaN; cbrt(-x) == -cbrt(x). Defined for all reals and never + * sets errno. + */ +double +cbrt(double x); + +float +cbrtf(float x); + +long double +cbrtl(long double x); + +/* + * Square root of x*x + y*y without undue overflow or underflow + * (C23 7.12.5.6). hypot(+-0, +-0) == +0; hypot(+-Inf, y) == +Inf even for + * NaN y; hypot(x, NaN) == NaN for finite x. The internal scaling avoids + * spurious overflow, so no errno path exists. + */ +double +hypot(double x, double y); + +float +hypotf(float x, float y); + +long double +hypotl(long double x, long double y); + +/* + * x raised to the power y (C23 7.12.6.11). pow(x, 0) == 1; pow(1, y) == 1 + * even for NaN y; pow(NaN, 0) == 1; pow(+-0, negative y) == +-Inf (a pole + * error setting errno to ERANGE); pow(negative finite x, non-integer y) is + * a domain error setting errno to EDOM; overflow and underflow set errno to + * ERANGE. No const attribute. + */ +double +pow(double x, double y); + +float +powf(float x, float y); + +long double +powl(long double x, long double y); + #ifdef __cplusplus } #endif diff --git a/src/math/cbrt.c b/src/math/cbrt.c new file mode 100644 index 0000000..1afc91c --- /dev/null +++ b/src/math/cbrt.c @@ -0,0 +1,236 @@ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifdef HAVE_CONFIG_H +#define SET_ERANGE() (errno = ERANGE) +#define SET_EDOM() (errno = EDOM) +#else +#define SET_ERANGE() \ + do \ + { \ + } while (0) +#define SET_EDOM() \ + do \ + { \ + } while (0) +#endif + +/* cbrt(x) = cube root of x, after musl src/math/cbrt.c (FreeBSD msun + * s_cbrt.c, Bruce D. Evans): a 5-bit bit-hack estimate, a degree-4 + * polynomial refinement to 23 bits, and one Newton step to 53 bits. */ +static const unsigned int B1 = 715094163; /* B1 = (1023-1023/3-0.03306235651)*2**20 */ +static const unsigned int B2 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */ + +/* |1/cbrt(x) - p(x)| < 2**-23.5 (~[-7.93e-8, 7.929e-8]). */ +static const double P0 = 1.87595182427177009643; /* 0x3ffe03e60f61e692 */ +static const double P1 = -1.88497979543377169875; /* 0xbffe28e092f02420 */ +static const double P2 = 1.621429720105354466140; /* 0x3ff9f1604a49d6c2 */ +static const double P3 = -0.758397934778766047437; /* 0xbfe844cbbee751d9 */ +static const double P4 = 0.145996192886612446982; /* 0x3fc2b000d4e4edd7 */ + +double +cbrt(double x) +{ + union + { + double f; + unsigned long long i; + } u; + double r; + double s; + double t; + double w; + unsigned int hx; + + u.f = x; + hx = (unsigned int)(u.i >> 32) & 0x7fffffffU; + if (hx >= 0x7ff00000U) + return x + x; /* cbrt(NaN,Inf) is itself */ + + /* rough cbrt to 5 bits */ + if (hx < 0x00100000U) + { /* zero or subnormal? */ + u.f = x * 0x1p54; + hx = (unsigned int)(u.i >> 32) & 0x7fffffffU; + if (hx == 0) + return x; /* cbrt(0) is itself */ + hx = hx / 3 + B2; + } + else + hx = hx / 3 + B1; + u.i &= 1ULL << 63; + u.i |= (unsigned long long)hx << 32; + t = u.f; + + /* new cbrt to 23 bits */ + r = (t * t) * (t / x); + t = t * ((P0 + r * (P1 + r * P2)) + ((r * r) * r) * (P3 + r * P4)); + + /* round t away from zero to 23 bits, larger in magnitude than cbrt(x) */ + u.f = t; + u.i = (u.i + 0x80000000ULL) & 0xffffffffc0000000ULL; + t = u.f; + + /* one Newton iteration to 53 bits with error < 0.667 ulps */ + s = t * t; /* t*t is exact */ + r = x / s; /* error <= 0.5 ulps; |r| < |t| */ + w = t + t; /* t+t is exact */ + r = (r - t) / (w + r); + t = t + t * r; + return t; +} + +/* cbrtf(x) = cube root of x, after musl src/math/cbrtf.c: the 5-bit + * bit-hack plus two Newton steps evaluated in double precision. */ +static const unsigned int B1f = 709958130; /* B1f = (127-127.0/3-0.03306235651)*2**23 */ +static const unsigned int B2f = 642849266; /* B2f = (127-127.0/3-24/3-0.03306235651)*2**23 */ + +float +cbrtf(float x) +{ + union + { + float f; + unsigned int i; + } u; + double r; + double T; + unsigned int hx; + + u.f = x; + hx = u.i & 0x7fffffffU; + if (hx >= 0x7f800000U) + return x + x; /* cbrt(NaN,Inf) is itself */ + + /* rough cbrt to 5 bits */ + if (hx < 0x00800000U) + { /* zero or subnormal? */ + if (hx == 0) + return x; /* cbrt(+-0) is itself */ + u.f = x * 0x1p24f; + hx = u.i & 0x7fffffffU; + hx = hx / 3 + B2f; + } + else + hx = hx / 3 + B1f; + u.i &= 0x80000000U; + u.i |= hx; + + /* first Newton step (solving t*t - x/t == 0) to 16 bits, in double */ + T = u.f; + r = T * T * T; + T = T * ((double)x + x + r) / (x + r + r); + + /* second Newton step to 47 bits; rounding to 24 bits is exact in + * round-to-nearest mode */ + r = T * T * T; + T = T * ((double)x + x + r) / (x + r + r); + return (float)T; +} + +/* cbrtl(x) = cube root of x, after musl src/math/cbrtl.c (ld80 branch, + * S. Kargl / B. Evans / D. Schultz): reduce x to [1,8) via e%3 and run the + * cbrtf-style estimates in double, then one long-double Newton step. */ +static const unsigned int B1l = 709958130; /* B1l = (127-127.0/3-0.03306235651)*2**23 */ + +long double +cbrtl(long double x) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } u; + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } v; + union + { + float f; + unsigned int i; + } uft; + long double r; + long double s; + long double t; + long double w; + double dr; + double dt; + double dx; + float ft; + int e; + int sign; + + u.f = x; + e = u.p.se & 0x7fff; + sign = u.p.se & 0x8000; + if (e == 0x7fff) + return x + x; /* cbrt(+-Inf) = +-Inf, cbrt(NaN) = NaN */ + if (e == 0) + { /* adjust subnormal numbers */ + u.f *= 0x1p120L; + e = u.p.se & 0x7fff; + if (e == 0) + return x; /* cbrt(+-0) = +-0 */ + e -= 120; + } + e -= 0x3fff; + u.p.se = 0x3fff; + x = u.f; + switch (e % 3) + { + case 1: + case -2: + x *= 2; + e--; + break; + case 2: + case -1: + x *= 4; + e -= 2; + break; + } + v.f = 1.0L; + v.p.se = (unsigned short)(sign | (0x3fff + e / 3)); + + /* ~5-bit estimate */ + uft.f = x; + uft.i = (uft.i & 0x7fffffffU) / 3 + B1l; + ft = uft.f; + + /* ~16-bit estimate */ + dx = x; + dt = ft; + dr = dt * dt * dt; + dt = dt * (dx + dx + dr) / (dx + dr + dr); + + /* ~47-bit estimate */ + dr = dt * dt * dt; + dt = dt * (dx + dx + dr) / (dx + dr + dr); + + /* round dt away from zero to 32 bits so t*t is exact */ + t = dt + (0x1.0p32L + 0x1.0p-31L) - 0x1.0p32; + + /* final Newton iteration to 64 bits with error < 0.667 ulps */ + s = t * t; /* t*t is exact */ + r = x / s; /* error <= 0.5 ulps; |r| < |t| */ + w = t + t; /* t+t is exact */ + r = (r - t) / (w + r); + t = t + t * r; + t *= v.f; + return t; +} diff --git a/src/math/exp.c b/src/math/exp.c new file mode 100644 index 0000000..d1ed2d3 --- /dev/null +++ b/src/math/exp.c @@ -0,0 +1,661 @@ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifdef HAVE_CONFIG_H +#define SET_ERANGE() (errno = ERANGE) +#else +#define SET_ERANGE() \ + do \ + { \ + } while (0) +#endif + +static const double ln2hi = 6.93147180369123816490e-01; /* 0x3fe62e42fee00000 */ +static const double ln2lo = 1.90821492927058770002e-10; /* 0x3dea39ef35793c76 */ +static const double invln2 = 1.44269504088896338700e+00; /* 0x3ff71547652b82fe */ +static const double P1 = 1.66666666666666019037e-01; /* 0x3FC555555555553E */ +static const double P2 = -2.77777777770155933842e-03; /* 0xBF66C16C16BEBD93 */ +static const double P3 = 6.61375632143793436117e-05; /* 0x3F11566AAF25DE2C */ +static const double P4 = -1.65339022054652515390e-06; /* 0xBEBBBD41C5D26BF1 */ +static const double P5 = 4.13813679705723846039e-08; /* 0x3E66376972BEA4D0 */ +static const double o_threshold = 709.782712893384; +static const double u_threshold = -745.133219101941; +static const double twom1000 = 9.33263618503218878990e-302; +static const double tiny = 1.0e-300; + +static double +vlibc_exp_core(double x) +{ + union + { + double f; + unsigned long long i; + } u; + + u.f = x; + unsigned long long sign = u.i >> 63; + unsigned int hx = (unsigned int)((u.i >> 32) & 0x7fffffffULL); + + /* NaN or +/-Inf: exp(-inf)=+0, exp(+inf)=+inf, exp(NaN)=NaN */ + if ((u.i & 0x7fffffffffffffffULL) >= 0x7ff0000000000000ULL) + { + if ((u.i & 0x000fffffffffffffULL) != 0ULL) + return x + x; /* NaN */ + return sign != 0ULL ? 0.0 : __builtin_huge_val(); + } + + if (x > o_threshold) + return __builtin_huge_val(); /* caller marks ERANGE */ + if (x < u_threshold) + return 0.0; /* caller marks ERANGE */ + + double hi = 0.0; + double lo = 0.0; + int k = 0; + + if (hx > 0x3fd62e42U) + { /* |x| > 0.5 ln2 */ + if (hx < 0x3ff0a2b2U) + { /* |x| < 1.5 ln2 */ + if (sign == 0ULL) + { + hi = x - ln2hi; + lo = ln2lo; + k = 1; + } + else + { + hi = x + ln2hi; + lo = -ln2lo; + k = -1; + } + } + else + { + k = (int)(invln2 * x + (sign != 0ULL ? -0.5 : 0.5)); + double t = (double)k; + hi = x - t * ln2hi; + lo = t * ln2lo; + } + x = hi - lo; + } + else if (hx < 0x3e300000U) + { /* |x| < 2^-28: exp(x) ~ 1+x */ + (void)tiny; + return 1.0 + x; + } + else + { + k = 0; + } + + double t = x * x; + double c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); + + if (k == 0) + return 1.0 - ((x * c) / (c - 2.0) - x); + + double y = 1.0 - ((lo - (x * c) / (2.0 - c)) - hi); + + /* scale y by 2^k by adding k to the exponent field */ + if (k >= -1021 && k <= 1023) + { + union + { + double f; + unsigned long long i; + } b; + b.f = y; + b.i += (unsigned long long)k << 52; + return b.f; + } + if (k == 1024) + return y * 2.0 * 0x1p1023; + /* k <= -1022: scale up into the normal range, then multiply down */ + union + { + double f; + unsigned long long i; + } b; + b.f = y; + b.i += (unsigned long long)(k + 1000) << 52; + return b.f * twom1000; +} + +double +exp(double x) +{ + if (x > o_threshold) + SET_ERANGE(); + else if (x < u_threshold) + SET_ERANGE(); + return vlibc_exp_core(x); +} + +#ifndef exp +double +exp(double); +#endif + +double +exp2(double x) +{ + union + { + double f; + unsigned long long i; + } u; + + u.f = x; + + /* NaN or +/-Inf: 2^(-inf)=+0, 2^(+inf)=+inf, 2^(NaN)=NaN */ + if ((u.i & 0x7fffffffffffffffULL) >= 0x7ff0000000000000ULL) + { + if ((u.i & 0x000fffffffffffffULL) != 0ULL) + return x; /* NaN */ + return (u.i >> 63) != 0ULL ? 0.0 : __builtin_huge_val(); + } + + if (x >= 1024.0) + { + SET_ERANGE(); + return __builtin_huge_val(); + } + if (x <= -1075.0) + { + SET_ERANGE(); + return 0.0; + } + + /* round to nearest integer; |x| < 1024 so the conversion is safe */ + long long k = (long long)(x + (x > 0.0 ? 0.5 : -0.5)); + /* exact by Sterbenz: k is within 2^52 of x */ + double r = x - (double)k; + /* e^(r ln2) = 2^r */ + double y = vlibc_exp_core(r * 0.6931471805599453); + + /* scale y by 2^k by adding k to the exponent field */ + if (k >= -1021 && k <= 1023) + { + union + { + double f; + unsigned long long i; + } b; + b.f = y; + b.i += (unsigned long long)k << 52; + return b.f; + } + if (k == 1024) + return y * 2.0 * 0x1p1023; + /* k <= -1022: scale up into the normal range, then multiply down */ + union + { + double f; + unsigned long long i; + } b; + b.f = y; + b.i += (unsigned long long)(k + 1000) << 52; + return b.f * twom1000; +} + +float +exp2f(float x) +{ + if (x >= 128.0f) + return __builtin_huge_valf(); + if (x <= -150.0f) + return 0.0f; + return (float)exp2((double)x); +} + +long double +exp2l(long double x) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } u; + + u.f = x; + if ((u.p.se & 0x7fffU) == 0x7fffU) + { + if (u.p.m != 0x8000000000000000ULL) + return x; /* NaN */ + return (u.p.se & 0x8000U) != 0 ? 0.0L : x; /* +/-Inf */ + } + + if (x >= 16384.0L) + { + SET_ERANGE(); + return __builtin_huge_vall(); + } + if (x <= -16445.0L) + { + SET_ERANGE(); + return 0.0L; + } + + /* round to nearest integer; |x| < 16384 so the conversion is safe */ + long long k = (long long)(x + (x > 0.0L ? 0.5L : -0.5L)); + /* exact by Sterbenz: k is within 2^63 of x */ + long double r = x - (long double)k; + long double y = r * 0.693147180559945309417232121458176568L; + + /* e^y, |y| <= ln2/2: Taylor series with compensated summation */ + long double s; + if (y > -1e-8L && y < 1e-8L) + s = 1.0L + y; + else + { + long double term = y; + long double comp = 0.0L; + s = 1.0L + y; + for (int i = 2; i <= 60; i++) + { + term *= y / (long double)i; + long double yk = term - comp; + long double t = s + yk; + comp = (t - s) - yk; + s = t; + if (term > -1e-30L && term < 1e-30L) + break; + } + } + + if (k < -16445) + return x > 0.0L ? 0.0L : -0.0L; /* unreachable; defensive */ + + /* scale s by 2^k: construct 2^k in long double and multiply */ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } two; + + if (k >= -16382 && k <= 16383) + { + two.p.m = 0x8000000000000000ULL; + two.p.se = (unsigned short)(16383 + (int)k); + return s * two.f; + } + if (k == 16384) + { + /* x in [16383.5, 16384): s = 2^r with r in [-0.5, 0) < 1, so + * s * 2^16383 * 2 stays finite */ + two.p.m = 0x8000000000000000ULL; + two.p.se = (unsigned short)(16383 + 16383); + return s * two.f * 2.0L; + } + /* k in [-16445, -16383]: 2^k is a pseudo-denormal, int bit clear */ + two.p.m = 1ULL << (unsigned int)(k + 16445); + two.p.se = 0; + return s * two.f; +} + +static const double Q1 = -3.33333333333331316428e-02; /* 0xbfa1111111110f54 */ +static const double Q2 = 1.58730158725481460165e-03; /* 0x3f5a01a019fe5585 */ +static const double Q3 = -7.93650757867487942473e-05; /* 0xbf14ce199eaadbb7 */ +static const double Q4 = 4.00821782732936260352e-06; /* 0x3ed0cfca86e65239 */ +static const double Q5 = -2.01099218183624371326e-07; /* 0xbe8afdb76e09c32d */ + +/* expm1(x) = e^x - 1, accurate for small x (fdlibm s_expm1). */ +double +expm1(double x) +{ + union + { + double f; + unsigned long long i; + } u; + double hi = 0.0; + double lo = 0.0; + double c = 0.0; + unsigned int sign; + unsigned int hx; + int k = 0; + + u.f = x; + sign = (unsigned int)(u.i >> 63); + hx = (unsigned int)((u.i >> 32) & 0x7fffffffULL); + + /* filter out huge and non-finite argument */ + if (hx >= 0x7ff00000U) + { + if ((u.i & 0xfffffffffffffULL) != 0ULL) + return x; /* NaN */ + return sign != 0U ? -1.0 : x; /* expm1(+-inf) = {-1, +inf} */ + } + if (x == 0.0) + return x; /* preserve -0.0 */ + + if (x > o_threshold) + { + SET_ERANGE(); + return __builtin_huge_val(); /* overflow */ + } + if (sign != 0U && hx >= 0x4043687AU) + return -1.0; /* x < -56*ln2, exp(x) < 2^-56 -> -1.0 */ + + /* argument reduction */ + if (hx > 0x3fd62e42U) + { + if (hx < 0x3ff0a2b2U) + { + if (sign == 0U) + { + hi = x - ln2hi; + lo = ln2lo; + k = 1; + } + else + { + hi = x + ln2hi; + lo = -ln2lo; + k = -1; + } + } + else + { + k = (int)(invln2 * x + (sign != 0U ? -0.5 : 0.5)); + { + double t = (double)k; + hi = x - t * ln2hi; /* t*ln2hi is exact */ + lo = t * ln2lo; + } + } + x = hi - lo; + c = (hi - x) - lo; + } + else if (hx < 0x3c900000U) + { + return x; /* |x| < 2^-54: expm1(x) == x */ + } + else + { + k = 0; + } + + /* x is now in primary range */ + { + double hfx = 0.5 * x; + double hxs = x * hfx; + double r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5)))); + double t = 3.0 - r1 * hfx; + double e = hxs * ((r1 - t) / (6.0 - x * t)); + + if (k == 0) + return x - (x * e - hxs); /* c is 0 */ + + e = (x * (e - c) - c) - hxs; + if (k == -1) + return 0.5 * (x - e) - 0.5; + if (k == 1) + { + if (x < -0.25) + return -2.0 * (e - (x + 0.5)); + return 1.0 + 2.0 * (x - e); + } + /* general: expm1(x) = 2^k * (1 - (e - x)) - 1 */ + { + double y = 1.0 - (e - x); + union + { + double f; + unsigned long long i; + } b; + if (k == 1024) + return y * 2.0 * 0x1p1023 - 1.0; + b.f = y; + if (k >= -1021) + b.i += (unsigned long long)k << 52; + else + { + b.i += (unsigned long long)(k + 1000) << 52; + return b.f * twom1000 - 1.0; + } + return b.f - 1.0; + } + } +} + +float +expm1f(float x) +{ + if (x >= 88.72283905206835f && !isinf(x)) + { + SET_ERANGE(); + return __builtin_huge_valf(); /* overflow */ + } + return (float)expm1((double)x); +} + +/* 80-bit x87 expm1l */ +long double +expm1l(long double x) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } u; + long double r; + long double sum; + long long k; + + u.f = x; + if ((u.p.se & 0x7fffU) == 0x7fffU) + { + if (u.p.m != 0x8000000000000000ULL) + return x; /* NaN */ + return (u.p.se & 0x8000U) != 0U ? -1.0L : x; /* expm1(+-inf) = {-1, +inf} */ + } + if (x == 0.0L) + return x; + + if (x > 11356.523406294143L) + { + SET_ERANGE(); + return __builtin_huge_vall(); + } + if (x < -16447.0L) + return -1.0L; + + /* 2^x = 2^k * e^r with r in [-ln2/2, ln2/2]; ln2 split hi/lo so that + * k*ln2hi is exact and r keeps full 64-bit precision */ + { + const long double ln2hi_l = 6.93147180369123816490e-01L; + const long double ln2lo_l = 1.90821492927058770002e-10L; + const long double invln2l = 1.442695040888963407359924681001892137L; + k = (long long)(x * invln2l + (x > 0.0L ? 0.5L : -0.5L)); + r = x - (long double)k * ln2hi_l - (long double)k * ln2lo_l; + } + + /* expm1(r) = e^r - 1 via Kahan-compensated Taylor sum */ + { + long double term = r; + long double comp = 0.0L; + int i; + sum = r; + for (i = 2; i <= 60; i++) + { + long double yt; + long double t; + term *= r / (long double)i; + if (term == 0.0L) + break; + yt = term - comp; + t = sum + yt; + comp = (t - sum) - yt; + sum = t; + } + } + + if (k == 0) + return sum; + + /* scale (sum + 1) by 2^k and subtract 1 */ + { + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } b; + long double z = sum + 1.0L; + + if (k >= 16384) + return z * 2.0L * 0x1p16383L - 1.0L; + if (k >= -16382) + { + b.p.m = 0x8000000000000000ULL; + b.p.se = (unsigned short)(16383 + k); + return z * b.f - 1.0L; + } + /* k in [-16445, -16383]: e^x < 2^-16382, so e^x - 1 rounds to -1 */ + return -1.0L; + } +} + +float +expf(float x) +{ + if (isinf(x)) + return x > 0.0f ? x : 0.0f; /* exp(+inf) = +inf, exp(-inf) = +0 */ + if (x >= 88.72283905206835f) + { + SET_ERANGE(); + return __builtin_huge_valf(); /* overflow */ + } + if (x < -103.97207708219999f) + { + SET_ERANGE(); + return 0.0f; /* underflow to zero */ + } + { + double y = exp((double)x); + if (y < 0x1p-149) /* glibc flags ERANGE below FLT_TRUE_MIN */ + SET_ERANGE(); + return (float)y; + } +} + +/* 80-bit x87 expl: exp(x) = 2^k * e^r, r in [-ln2/2, ln2/2] */ +long double +expl(long double x) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } u; + long double r; + long double sum; + long long k; + + u.f = x; + if ((u.p.se & 0x7fffU) == 0x7fffU) + { + if (u.p.m != 0x8000000000000000ULL) + return x; /* NaN */ + return (u.p.se & 0x8000U) != 0U ? 0.0L : x; /* exp(+-inf) = {+0, inf} */ + } + if (x == 0.0L) + return 1.0L; + + if (x >= 0xb.17217f7d1cf79acp+10L) /* first 80-bit x where e^x rounds to +inf */ + { + SET_ERANGE(); + return __builtin_huge_vall(); + } + if (x <= -11399.8484L) /* exp(x) < 2^-16446: rounds to +0 */ + { + SET_ERANGE(); + return 0.0L; + } + + /* ln2 split so k*ln2hi is exact (15-bit hi) and r keeps full + * 64-bit precision; k fits long long since |x| < 11356.6 */ + { + const long double ln2hi_l = 6.9314575195312500000000e-1L; + const long double ln2lo_l = 1.4286068203094172321214581765680755e-6L; + const long double invln2l = 1.442695040888963407359924681001892137L; + k = (long long)(x * invln2l + (x > 0.0L ? 0.5L : -0.5L)); + r = x - (long double)k * ln2hi_l; + r -= (long double)k * ln2lo_l; + } + + /* e^r via Kahan-compensated Taylor */ + { + long double term = r; + long double comp = 0.0L; + int i; + sum = 1.0L + r; + for (i = 2; i <= 60; i++) + { + long double yt; + long double t; + term *= r / (long double)i; + if (term == 0.0L) + break; + yt = term - comp; + t = sum + yt; + comp = (t - sum) - yt; + sum = t; + } + } + + /* scale sum by 2^k */ + { + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } b; + + if (k == 0) + return sum; + if (k >= 16384) + return sum * 2.0L * 0x1p16383L; /* overflow path */ + if (k <= -16382) + { + /* result is subnormal: scale up by 2^(k+16382), then by 2^-16382 */ + b.p.m = 0x8000000000000000ULL; + b.p.se = (unsigned short)(16383 + k + 16382); + sum *= b.f; + return sum * 0x1p-16382L; + } + b.p.m = 0x8000000000000000ULL; + b.p.se = (unsigned short)(16383 + k); + return sum * b.f; + } +} diff --git a/src/math/hypot.c b/src/math/hypot.c new file mode 100644 index 0000000..140dbc5 --- /dev/null +++ b/src/math/hypot.c @@ -0,0 +1,265 @@ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifdef HAVE_CONFIG_H +#define SET_ERANGE() (errno = ERANGE) +#define SET_EDOM() (errno = EDOM) +#else +#define SET_ERANGE() \ + do \ + { \ + } while (0) +#define SET_EDOM() \ + do \ + { \ + } while (0) +#endif + +/* hypot(x,y) = sqrt(x*x + y*y) without spurious overflow/underflow, after + * musl src/math/hypot.c: |x| >= |y| via union sort, inf/NaN specials, a + * power-of-two scaling window, and an exact double-double square so the + * final sqrt argument is precise. No errno: the scaling avoids overflow. */ +static const double split = 0x1p27 + 1.0; /* FLT_EVAL_METHOD <= 1 on x86-64 */ + +static double +vlibc_hypot_core(double x, double y) +{ + union + { + double f; + unsigned long long i; + } ux, uy, ut; + int ex; + int ey; + double hx; + double lx; + double hy; + double ly; + double z; + double xc; + double xh; + double xl; + + /* arrange |x| >= |y| */ + ux.f = x; + uy.f = y; + ux.i &= 0x7fffffffffffffffULL; + uy.i &= 0x7fffffffffffffffULL; + if (ux.i < uy.i) + { + ut = ux; + ux = uy; + uy = ut; + } + ex = (int)(ux.i >> 52); + ey = (int)(uy.i >> 52); + x = ux.f; + y = uy.f; + /* note: hypot(inf,nan) == inf */ + if (ey == 0x7ff) + return y; + if (ex == 0x7ff || uy.i == 0) + return x; + /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */ + if (ex - ey > 64) + return x + y; + + /* precise sqrt argument in nearest rounding mode without overflow: + * xh*xh must not overflow and xl*xl must not underflow */ + z = 1.0; + if (ex > 0x3ff + 510) + { + z = 0x1p700; + x *= 0x1p-700; + y *= 0x1p-700; + } + else if (ey < 0x3ff - 450) + { + z = 0x1p-700; + x *= 0x1p700; + y *= 0x1p700; + } + /* hx+lx == x*x exactly (double-double square) */ + xc = x * split; + xh = x - xc + xc; + xl = x - xh; + hx = x * x; + lx = xh * xh - hx + 2.0 * xh * xl + xl * xl; + xc = y * split; + xh = y - xc + xc; + xl = y - xh; + hy = y * y; + ly = xh * xh - hy + 2.0 * xh * xl + xl * xl; + x = ly + lx + hy + hx; + __asm__("sqrtsd %1, %0" : "=x"(y) : "x"(x)); + return z * y; +} + +/* hypotf after musl src/math/hypotf.c: same shape, scaled into a safe + * window and evaluated in double for the final square and sqrt. */ +static float +vlibc_hypotf_core(float x, float y) +{ + union + { + float f; + unsigned int i; + } ux, uy, ut; + float z; + + ux.f = x; + uy.f = y; + ux.i &= 0x7fffffffU; + uy.i &= 0x7fffffffU; + if (ux.i < uy.i) + { + ut = ux; + ux = uy; + uy = ut; + } + x = ux.f; + y = uy.f; + if (uy.i == 0xffU << 23) + return y; /* hypot(inf,nan) == inf */ + if (ux.i >= 0xffU << 23 || uy.i == 0 || ux.i - uy.i >= 25U << 23) + return x + y; + z = 1.0f; + if (ux.i >= (0x7fU + 60) << 23) + { + z = 0x1p90f; + x *= 0x1p-90f; + y *= 0x1p-90f; + } + else if (uy.i < (0x7fU - 60) << 23) + { + z = 0x1p-90f; + x *= 0x1p90f; + y *= 0x1p90f; + } + x = (float)((double)x * x + (double)y * y); + __asm__("sqrtss %1, %0" : "=x"(y) : "x"(x)); + return z * y; +} + +/* hypotl after musl src/math/hypotl.c (ld80 branch): the full 64-bit + * mantissa path with a 2^10000 scaling window and an exact double-double + * square. */ +static long double +vlibc_hypotl_core(long double x, long double y) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } ux; + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } uy; + int ex; + int ey; + long double hx; + long double lx; + long double hy; + long double ly; + long double z; + long double xc; + long double xh; + long double xl; + + ux.f = x; + uy.f = y; + ux.p.se &= 0x7fff; + uy.p.se &= 0x7fff; + if (ux.p.se < uy.p.se) + { + ex = uy.p.se; + ey = ux.p.se; + x = uy.f; + y = ux.f; + } + else + { + ex = ux.p.se; + ey = uy.p.se; + x = ux.f; + y = uy.f; + } + if (ex == 0x7fff && __builtin_isinf(y)) + return y; /* hypot(inf,nan) == inf */ + if (ex == 0x7fff || y == 0.0L) + return x; + if (ex - ey > 64) + return x + y; + z = 1.0L; + if (ex > 0x3fff + 8000) + { + z = 0x1p10000L; + x *= 0x1p-10000L; + y *= 0x1p-10000L; + } + else if (ey < 0x3fff - 8000) + { + z = 0x1p-10000L; + x *= 0x1p10000L; + y *= 0x1p10000L; + } + /* hx+lx == x*x exactly (double-double square) */ + xc = x * (0x1p32L + 1.0L); + xh = x - xc + xc; + xl = x - xh; + hx = x * x; + lx = xh * xh - hx + 2.0L * xh * xl + xl * xl; + xc = y * (0x1p32L + 1.0L); + xh = y - xc + xc; + xl = y - xh; + hy = y * y; + ly = xh * xh - hy + 2.0L * xh * xl + xl * xl; + x = ly + lx + hy + hx; + __asm__("fsqrt" : "+t"(x)); + return z * x; +} + +double +hypot(double x, double y) +{ + double z = vlibc_hypot_core(x, y); + + if (z == __builtin_huge_val() && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_ERANGE(); /* genuine overflow: glibc parity */ + return z; +} + +float +hypotf(float x, float y) +{ + float z = vlibc_hypotf_core(x, y); + + if (z == __builtin_huge_valf() && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_ERANGE(); + return z; +} + +long double +hypotl(long double x, long double y) +{ + long double z = vlibc_hypotl_core(x, y); + + if (z == __builtin_huge_vall() && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_ERANGE(); + return z; +} diff --git a/src/math/log.c b/src/math/log.c new file mode 100644 index 0000000..3b6f8a3 --- /dev/null +++ b/src/math/log.c @@ -0,0 +1,457 @@ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifdef HAVE_CONFIG_H +#define SET_ERANGE() (errno = ERANGE) +#define SET_EDOM() (errno = EDOM) +#else +#define SET_ERANGE() \ + do \ + { \ + } while (0) +#define SET_EDOM() \ + do \ + { \ + } while (0) +#endif + +static const double two54 = 1.80143985094819840000e+16; /* 0x4350000000000000 */ +static const double ln2_hi = 6.93147180369123816490e-01; /* 0x3fe62e42fee00000 */ +static const double ln2_lo = 1.90821492927058770002e-10; /* 0x3dea39ef35793c76 */ +static const double Lg1 = 6.666666666666735130e-01; /* 0x3FE5555555555593 */ +static const double Lg2 = 3.999999999940941908e-01; /* 0x3FD999999997FA04 */ +static const double Lg3 = 2.857142874366239149e-01; /* 0x3FD2492494229359 */ +static const double Lg4 = 2.222219843214978396e-01; /* 0x3FCC71C51D8E78AF */ +static const double Lg5 = 1.818357216161805012e-01; /* 0x3FC7466496CB03DE */ +static const double Lg6 = 1.531383769920937332e-01; /* 0x3FC39A09D078C69F */ +static const double Lg7 = 1.479819860511658591e-01; /* 0x3FC2F112DF3E5244 */ + +/* 80-bit ld80 ln2 split: ln2hi_l is a 15-bit hi word (0xb.172p-4L) so that + * k * ln2hi_l is exact for |k| < 2^15, and ln2hi_l + ln2lo_l == ln2. */ +static const long double ln2hi_l = 6.9314575195312500000000e-1L; /* 0xb.172p-4L */ +static const long double ln2lo_l = 1.4286068203094172321214581765680755e-6L; +static const long double ivln2_l = 0x1.71547652b82fe178p+0L; /* 1/ln2, 64-bit */ +static const long double ivln10_l = 0x1.bcb7b1526e50e32ap-2L; /* 1/ln10, 64-bit */ + +/* log(x) = k * ln2 + log(1+f), fdlibm e_log with the double split into + * ln2_hi/ln2_lo so k*ln2_hi is exact for |k| <= 1024. */ +static double +vlibc_log_core(double x) +{ + union + { + double f; + unsigned long long i; + } u; + unsigned int hx; + unsigned int lx; + int k; + int i; + int j; + double f; + double s; + double z; + double R; + double w; + double t1; + double t2; + double hfsq; + double dk; + + u.f = x; + hx = (unsigned int)(u.i >> 32); + lx = (unsigned int)u.i; + k = 0; + + if (hx >= 0x80000000U) + { /* sign bit set: -0 or negative */ + if (((hx & 0x7fffffffU) | lx) == 0U) + return -__builtin_huge_val(); /* log(-0) = -inf */ + return __builtin_nan(""); /* log(negative) = NaN */ + } + if ((hx | lx) == 0U) + return -__builtin_huge_val(); /* log(+0) = -inf */ + if (hx >= 0x7ff00000U) + return x + x; /* +inf -> +inf, NaN -> NaN */ + + if (hx < 0x00100000U) + { /* subnormal: scale into the normal range */ + k -= 54; + x *= two54; + u.f = x; + hx = (unsigned int)(u.i >> 32); + lx = (unsigned int)u.i; + } + k += (int)(hx >> 20) - 1023; + hx &= 0x000fffffU; + i = (int)((hx + 0x95f64U) & 0x100000U); + u.i = ((unsigned long long)(hx | (i ^ 0x3ff00000U)) << 32) | (unsigned long long)lx; + x = u.f; + k += i >> 20; + f = x - 1.0; + + if ((0x000fffffU & (2U + hx)) < 3U) + { /* |f| < 2^-20 */ + if (f == 0.0) + { + if (k == 0) + return 0.0; + dk = (double)k; + return dk * ln2_hi + dk * ln2_lo; + } + R = f * f * (0.5 - 0.33333333333333333333 * f); + if (k == 0) + return f - R; + dk = (double)k; + return dk * ln2_hi - ((R - dk * ln2_lo) - f); + } + s = f / (2.0 + f); + dk = (double)k; + z = s * s; + i = (int)hx - 0x6147a; + w = z * z; + j = 0x6b851 - (int)hx; + t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); + t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); + i |= j; + R = t2 + t1; + if (i > 0) + { + hfsq = 0.5 * f * f; + if (k == 0) + return f - (hfsq - s * (hfsq + R)); + return dk * ln2_hi - ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) - f); + } + if (k == 0) + return f - s * (f - R); + return dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f); +} + +double +log(double x) +{ + if (x < 0.0) + SET_EDOM(); + else if (x == 0.0) + SET_ERANGE(); + return vlibc_log_core(x); +} + +float +logf(float x) +{ + return (float)log((double)x); +} + +/* z = A/(2^64 + A) for 0 <= A <= 2^62, correctly rounded to long double, + * via integer long division of A^2/(2^64 + A) */ +static long double +vlibc_logl_ratio(unsigned long long A) +{ + unsigned long long a0; + unsigned long long a1; + unsigned long long p0; + unsigned long long p1; + unsigned long long p2; + unsigned long long lo; + unsigned long long hi; + unsigned long long q; + unsigned long long rhi; + unsigned long long rlo; + unsigned long long zq; + int msb; + int e; + int j; + unsigned long long dj; + unsigned long long shi; + unsigned long long slo; + unsigned long long M; + unsigned long long c2; + unsigned long long l2; + int i; + union + { + long double f; + struct + { + unsigned long long mm; + unsigned short se; + } p; + } two; + + a0 = A & 0xffffffffULL; + a1 = A >> 32; + p0 = a0 * a0; + p1 = (a1 * a0) << 1; + p2 = a1 * a1; + lo = p0 + (p1 << 32); + hi = p2 + (p1 >> 32) + (lo < p0); + + /* q = floor(A^2 / (2^64 + A)); the top dividend bits shift in with no + * quotient bits (q < 2^60), so start with R = hi and shift in lo */ + q = 0; + rhi = 0; + rlo = hi; + for (i = 63; i >= 0; i--) + { + unsigned long long nhi = (rhi << 1) | (rlo >> 63); + unsigned long long nlo = (rlo << 1) | ((lo >> i) & 1ULL); + if (nhi > 1 || (nhi == 1 && nlo >= A)) + { + unsigned long long nnlo = nlo - A; + unsigned long long borrow = nnlo > nlo; + nlo = nnlo; + nhi -= 1 + borrow; + if (i < 60) + q |= 1ULL << i; + } + rhi = nhi; + rlo = nlo; + } + + /* z*2^64 = zq - R/D with R = (rhi,rlo) in [0, D) */ + zq = A - q; + if (zq == 0) + return 0.0L; + if (A == 1) + { /* z = (2^64 - 1) * 2^-128: below the general path's reach */ + two.p.mm = 0xffffffffffffffffULL; + two.p.se = (unsigned short)(16383 - 65); + return two.f; + } + /* e = floor(log2 z): z*2^64 is in (zq-1, zq]; significand is + * M = round(z * 2^(63-e)) = round(z*2^64 * 2^j) with j = -1-e */ + msb = 63 - __builtin_clzll((rhi != 0 || rlo != 0) ? zq - 1 : zq); + e = msb - 64; + j = -1 - e; + /* floor(R * 2^j / D) via j more long-division steps */ + dj = 0; + shi = rhi; + slo = rlo; + for (i = j - 1; i >= 0; i--) + { + unsigned long long nhi = (shi << 1) | (slo >> 63); + unsigned long long nlo = slo << 1; + if (nhi > 1 || (nhi == 1 && nlo >= A)) + { + unsigned long long nnlo = nlo - A; + unsigned long long borrow = nnlo > nlo; + nlo = nnlo; + nhi -= 1 + borrow; + dj |= 1ULL << i; + } + shi = nhi; + slo = nlo; + } + /* M = zq*2^j - dj - r_j/D; round to nearest even */ + M = (zq << j) - dj; + c2 = (shi << 1) + (slo >> 63); + l2 = slo << 1; + if (c2 > 1 || (c2 == 1 && l2 > A) || (c2 == 1 && l2 == A && (M & 1ULL))) + M--; + two.p.mm = 0x8000000000000000ULL; + two.p.se = (unsigned short)(16383 + e - 63); + return (long double)M * two.f; +} + +/* 80-bit x87 log: normalize m to [2^63, 2^64), then + * log(x) = k * ln2 + 2*atanh(f/(2+f)). The k*ln2 part is exact and the + * atanh sum is accumulated as a double-double pair so the result keeps + * full 64-bit precision even when k*ln2 and ln(1+f) cancel. The caller + * scales by c to get log, log2, or log10. */ +static long double +vlibc_logl_scaled(long double x, long double c) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } u; + unsigned long long m; + unsigned short se; + int k; + int s; + int n; + long double f; + long double z; + long double w; + long double term; + long double sh; + long double sl; + long double t; + long double e; + long double s2; + long double zz; + long double a1; + long double a2; + long double b1; + long double b2; + long double hi; + long double lo; + + u.f = x; + se = u.p.se; + m = u.p.m; + if (m == 0) + return -__builtin_huge_vall(); /* +-0 -> -inf */ + if (se != 0U && m < 0x8000000000000000ULL) + return __builtin_nanl(""); /* unnormal: invalid operand, glibc NaN */ + if ((se & 0x8000U) != 0U) + return __builtin_nanl(""); /* negative incl -inf -> NaN */ + if ((se & 0x7fffU) == 0x7fffU) + return x; /* +inf, NaN */ + + if (m < 0x8000000000000000ULL) + { /* subnormal or pseudo-denormal: normalize */ + s = __builtin_clzll(m); + m <<= s; + } + else + s = 0; + k = (se == 0U ? -16382 : (int)(se & 0x7fffU) - 16383) - s; + + f = (long double)m * 0x1p-63L - 1.0L; /* exact, f in [0,1) */ + if (f >= 0.5L) + { /* fold [0.5,1) into [-0.25,0) so k*ln2 + ln(1+f) does not cancel: + * z = -v/(1-v) = -(v + v^2 + v^3 + ...) with v = (2^64 - m)*2^-65 in + * (0,1/8], accumulated with TwoSum so the fast convergence keeps z + * accurate (a plain division would round its 65-bit denominator) */ + unsigned long long C = 0ULL - m; + long double v = (long double)C * 0x1p-65L; + long double term; + long double sh; + long double sl; + long double t; + long double e; + long double s2; + long double zz; + int n; + + k += 1; + sh = -v; + sl = 0.0L; + term = v; + for (n = 2; n <= 43; n++) + { + term *= v; + if (term == 0.0L) + break; + t = -term; + s2 = sh + t; + zz = s2 - sh; + e = (sh - (s2 - zz)) + (t - zz); + sh = s2; + sl += e; + } + z = sh + sl; + } + else + z = vlibc_logl_ratio(m - 0x8000000000000000ULL); + w = z * z; + term = z; + sh = z; + sl = 0.0L; + for (n = 3; n <= 43; n += 2) + { /* TwoSum accumulation: (sh,sl) == sum of terms to ~2^-126 */ + term *= w; + if (term == 0.0L) + break; + t = term / (long double)n; + s2 = sh + t; + zz = s2 - sh; + e = (sh - (s2 - zz)) + (t - zz); + sh = s2; + sl += e; + } + a1 = (long double)k * ln2hi_l; /* exact: |k| < 2^15, ln2hi_l is 15 bits */ + a2 = (long double)k * ln2lo_l; + b1 = 2.0L * sh; + b2 = 2.0L * sl; + /* double-double add (a1,a2) + (b1,b2), then renormalize */ + s2 = a1 + b1; + zz = s2 - a1; + e = (a1 - (s2 - zz)) + (b1 - zz); + hi = s2; + lo = e + a2 + b2; + s2 = hi + lo; + lo -= s2 - hi; + hi = s2; + return hi * c + lo * c; +} + +static long double +vlibc_logl_core(long double x) +{ + return vlibc_logl_scaled(x, 1.0L); +} + +long double +logl(long double x) +{ + if (x < 0.0L) + SET_EDOM(); + else if (x == 0.0L) + SET_ERANGE(); + return vlibc_logl_core(x); +} + +double +log2(double x) +{ + if (x < 0.0) + SET_EDOM(); + else if (x == 0.0) + SET_ERANGE(); + return (double)vlibc_logl_scaled((long double)x, ivln2_l); +} + +float +log2f(float x) +{ + return (float)log2((double)x); +} + +long double +log2l(long double x) +{ + if (x < 0.0L) + SET_EDOM(); + else if (x == 0.0L) + SET_ERANGE(); + return vlibc_logl_scaled(x, ivln2_l); +} + +double +log10(double x) +{ + if (x < 0.0) + SET_EDOM(); + else if (x == 0.0) + SET_ERANGE(); + return (double)vlibc_logl_scaled((long double)x, ivln10_l); +} + +float +log10f(float x) +{ + return (float)log10((double)x); +} + +long double +log10l(long double x) +{ + if (x < 0.0L) + SET_EDOM(); + else if (x == 0.0L) + SET_ERANGE(); + return vlibc_logl_scaled(x, ivln10_l); +} diff --git a/src/math/log1p.c b/src/math/log1p.c new file mode 100644 index 0000000..c20e203 --- /dev/null +++ b/src/math/log1p.c @@ -0,0 +1,427 @@ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifdef HAVE_CONFIG_H +#define SET_ERANGE() (errno = ERANGE) +#define SET_EDOM() (errno = EDOM) +#else +#define SET_ERANGE() \ + do \ + { \ + } while (0) +#define SET_EDOM() \ + do \ + { \ + } while (0) +#endif + +static const double ln2_hi = 6.93147180369123816490e-01; /* 0x3fe62e42fee00000 */ +static const double ln2_lo = 1.90821492927058770002e-10; /* 0x3dea39ef35793c76 */ +static const double Lg1 = 6.666666666666735130e-01; /* 0x3FE5555555555593 */ +static const double Lg2 = 3.999999999940941908e-01; /* 0x3FD999999997FA04 */ +static const double Lg3 = 2.857142874366239149e-01; /* 0x3FD2492494229359 */ +static const double Lg4 = 2.222219843214978396e-01; /* 0x3FCC71C51D8E78AF */ +static const double Lg5 = 1.818357216161805012e-01; /* 0x3FC7466496CB03DE */ +static const double Lg6 = 1.531383769920937332e-01; /* 0x3FC39A09D078C69F */ +static const double Lg7 = 1.479819860511658591e-01; /* 0x3FC2F112DF3E5244 */ + +/* 80-bit ld80 ln2 split: ln2hi_l is a 15-bit hi word (0xb.172p-4L) so that + * k * ln2hi_l is exact for |k| < 2^15, and ln2hi_l + ln2lo_l == ln2. */ +static const long double ln2hi_l = 6.9314575195312500000000e-1L; /* 0xb.172p-4L */ +static const long double ln2lo_l = 1.4286068203094172321214581765680755e-6L; + +/* log1p(x) = log(1+x) for the double case: fdlibm e_log1p with the argument + * reduced into [sqrt(2)/2-1, sqrt(2)-1] so no precision is lost for small x. */ +static double +vlibc_log1p_core(double x) +{ + union + { + double f; + unsigned long long i; + } u; + unsigned int hx; + unsigned int hu; + int k; + double hfsq; + double f; + double c; + double s; + double z; + double R; + double w; + double t1; + double t2; + double dk; + + u.f = x; + hx = (unsigned int)(u.i >> 32); + k = 1; + if (hx < 0x3fda827aU || (hx >> 31) != 0U) + { /* 1+x < sqrt(2)+, or x < 0 */ + if (hx >= 0xbff00000U) + { /* x <= -1.0 */ + if (x == -1.0) + return -__builtin_huge_val(); /* log1p(-1) = -inf */ + return __builtin_nan(""); /* log1p(x<-1) = NaN */ + } + if (hx << 1 < 0x3ca00000U << 1) + return x; /* |x| < 2^-53: log1p(x) == x */ + if (hx <= 0xbfd2bec4U) + { /* x in [-sqrt(2)/2+1, sqrt(2)-1]: no reduction needed */ + k = 0; + c = 0.0; + f = x; + } + } + else if (hx >= 0x7ff00000U) + return x; /* +inf -> +inf, NaN -> NaN */ + + if (k != 0) + { /* reduce 1+x into [sqrt(2)/2, sqrt(2)] */ + u.f = 1.0 + x; + hu = (unsigned int)(u.i >> 32); + hu += 0x3ff00000U - 0x3fe6a09eU; + k = (int)(hu >> 20) - 0x3ff; + if (k < 54) + { + c = k >= 2 ? 1.0 - (u.f - x) : x - (u.f - 1.0); + c /= u.f; + } + else + c = 0.0; + hu = (hu & 0x000fffffU) + 0x3fe6a09eU; + u.i = ((unsigned long long)hu << 32) | (u.i & 0xffffffffULL); + f = u.f - 1.0; + } + hfsq = 0.5 * f * f; + s = f / (2.0 + f); + z = s * s; + w = z * z; + t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); + t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); + R = t2 + t1; + dk = (double)k; + return s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi; +} + +double +log1p(double x) +{ + if (x < -1.0) + SET_EDOM(); + else if (x == -1.0) + SET_ERANGE(); + return vlibc_log1p_core(x); +} + +float +log1pf(float x) +{ + return (float)log1p((double)x); +} + +/* z = A/(2^64 + A) for 0 <= A <= 2^62, correctly rounded to long double, + * via integer long division of A^2/(2^64 + A) */ +static long double +vlibc_logl_ratio(unsigned long long A) +{ + unsigned long long a0; + unsigned long long a1; + unsigned long long p0; + unsigned long long p1; + unsigned long long p2; + unsigned long long lo; + unsigned long long hi; + unsigned long long q; + unsigned long long rhi; + unsigned long long rlo; + unsigned long long zq; + int msb; + int e; + int j; + unsigned long long dj; + unsigned long long shi; + unsigned long long slo; + unsigned long long M; + unsigned long long c2; + unsigned long long l2; + int i; + union + { + long double f; + struct + { + unsigned long long mm; + unsigned short se; + } p; + } two; + + a0 = A & 0xffffffffULL; + a1 = A >> 32; + p0 = a0 * a0; + p1 = (a1 * a0) << 1; + p2 = a1 * a1; + lo = p0 + (p1 << 32); + hi = p2 + (p1 >> 32) + (lo < p0); + + /* q = floor(A^2 / (2^64 + A)); the top dividend bits shift in with no + * quotient bits (q < 2^60), so start with R = hi and shift in lo */ + q = 0; + rhi = 0; + rlo = hi; + for (i = 63; i >= 0; i--) + { + unsigned long long nhi = (rhi << 1) | (rlo >> 63); + unsigned long long nlo = (rlo << 1) | ((lo >> i) & 1ULL); + if (nhi > 1 || (nhi == 1 && nlo >= A)) + { + unsigned long long nnlo = nlo - A; + unsigned long long borrow = nnlo > nlo; + nlo = nnlo; + nhi -= 1 + borrow; + if (i < 60) + q |= 1ULL << i; + } + rhi = nhi; + rlo = nlo; + } + + /* z*2^64 = zq - R/D with R = (rhi,rlo) in [0, D) */ + zq = A - q; + if (zq == 0) + return 0.0L; + if (A == 1) + { /* z = (2^64 - 1) * 2^-128: below the general path's reach */ + two.p.mm = 0xffffffffffffffffULL; + two.p.se = (unsigned short)(16383 - 65); + return two.f; + } + /* e = floor(log2 z): z*2^64 is in (zq-1, zq]; significand is + * M = round(z * 2^(63-e)) = round(z*2^64 * 2^j) with j = -1-e */ + msb = 63 - __builtin_clzll((rhi != 0 || rlo != 0) ? zq - 1 : zq); + e = msb - 64; + j = -1 - e; + /* floor(R * 2^j / D) via j more long-division steps */ + dj = 0; + shi = rhi; + slo = rlo; + for (i = j - 1; i >= 0; i--) + { + unsigned long long nhi = (shi << 1) | (slo >> 63); + unsigned long long nlo = slo << 1; + if (nhi > 1 || (nhi == 1 && nlo >= A)) + { + unsigned long long nnlo = nlo - A; + unsigned long long borrow = nnlo > nlo; + nlo = nnlo; + nhi -= 1 + borrow; + dj |= 1ULL << i; + } + shi = nhi; + slo = nlo; + } + /* M = zq*2^j - dj - r_j/D; round to nearest even */ + M = (zq << j) - dj; + c2 = (shi << 1) + (slo >> 63); + l2 = slo << 1; + if (c2 > 1 || (c2 == 1 && l2 > A) || (c2 == 1 && l2 == A && (M & 1ULL))) + M--; + two.p.mm = 0x8000000000000000ULL; + two.p.se = (unsigned short)(16383 + e - 63); + return (long double)M * two.f; +} + +/* 80-bit x87 log: normalize m to [2^63, 2^64), then + * log(x) = k * ln2 + 2*atanh(f/(2+f)). The k*ln2 part is exact and the + * atanh sum is accumulated as a double-double pair so the result keeps + * full 64-bit precision even when k*ln2 and ln(1+f) cancel. */ +static long double +vlibc_logl_core(long double x) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } u; + unsigned long long m; + unsigned short se; + int k; + int s; + int n; + long double f; + long double z; + long double w; + long double term; + long double sh; + long double sl; + long double t; + long double e; + long double s2; + long double zz; + long double a1; + long double a2; + long double b1; + long double b2; + long double hi; + long double lo; + + u.f = x; + se = u.p.se; + m = u.p.m; + if (m == 0) + return -__builtin_huge_vall(); /* +-0 -> -inf */ + if (se != 0U && m < 0x8000000000000000ULL) + return __builtin_nanl(""); /* unnormal: invalid operand, glibc NaN */ + if ((se & 0x8000U) != 0U) + return __builtin_nanl(""); /* negative incl -inf -> NaN */ + if ((se & 0x7fffU) == 0x7fffU) + return x; /* +inf, NaN */ + + if (m < 0x8000000000000000ULL) + { /* subnormal or pseudo-denormal: normalize */ + s = __builtin_clzll(m); + m <<= s; + } + else + s = 0; + k = (se == 0U ? -16382 : (int)(se & 0x7fffU) - 16383) - s; + + f = (long double)m * 0x1p-63L - 1.0L; /* exact, f in [0,1) */ + if (f >= 0.5L) + { /* fold [0.5,1) into [-0.25,0) so k*ln2 + ln(1+f) does not cancel: + * z = -v/(1-v) = -(v + v^2 + v^3 + ...) with v = (2^64 - m)*2^-65 in + * (0,1/8], accumulated with TwoSum so the fast convergence keeps z + * accurate (a plain division would round its 65-bit denominator) */ + unsigned long long C = 0ULL - m; + long double v = (long double)C * 0x1p-65L; + long double term; + long double sh; + long double sl; + long double t; + long double e; + long double s2; + long double zz; + int n; + + k += 1; + sh = -v; + sl = 0.0L; + term = v; + for (n = 2; n <= 43; n++) + { + term *= v; + if (term == 0.0L) + break; + t = -term; + s2 = sh + t; + zz = s2 - sh; + e = (sh - (s2 - zz)) + (t - zz); + sh = s2; + sl += e; + } + z = sh + sl; + } + else + z = vlibc_logl_ratio(m - 0x8000000000000000ULL); + w = z * z; + term = z; + sh = z; + sl = 0.0L; + for (n = 3; n <= 43; n += 2) + { /* TwoSum accumulation: (sh,sl) == sum of terms to ~2^-126 */ + term *= w; + if (term == 0.0L) + break; + t = term / (long double)n; + s2 = sh + t; + zz = s2 - sh; + e = (sh - (s2 - zz)) + (t - zz); + sh = s2; + sl += e; + } + a1 = (long double)k * ln2hi_l; /* exact: |k| < 2^15, ln2hi_l is 15 bits */ + a2 = (long double)k * ln2lo_l; + b1 = 2.0L * sh; + b2 = 2.0L * sl; + /* double-double add (a1,a2) + (b1,b2), then renormalize */ + s2 = a1 + b1; + zz = s2 - a1; + e = (a1 - (s2 - zz)) + (b1 - zz); + hi = s2; + lo = e + a2 + b2; + s2 = hi + lo; + lo -= s2 - hi; + hi = s2; + return hi + lo; +} + +/* 80-bit x87 log1p: y = 1+x as an exact double-double pair (yh,yl), then + * log1p(x) = log(yh) + log(1 + yl/yh) = log(yh) + t with |t| <= 2^-64, so + * 1+x never rounds away the small argument. */ +static long double +vlibc_log1pl_core(long double x) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } u; + unsigned long long m; + unsigned short se; + long double yh; + long double yl; + long double t; + + u.f = x; + se = u.p.se; + m = u.p.m; + if (m == 0) + return x; /* +-0 -> +-0 */ + if (se != 0U && m < 0x8000000000000000ULL) + return __builtin_nanl(""); /* unnormal: invalid operand, glibc NaN */ + if ((se & 0x7fffU) == 0x7fffU) + { + if (m == 0x8000000000000000ULL) + return (se & 0x8000U) != 0U ? __builtin_nanl("") : x; /* -inf, +inf */ + return x; /* NaN */ + } + if (x <= -1.0L) + return x == -1.0L ? -__builtin_huge_vall() : __builtin_nanl(""); + if (x > -0x1p-65L && x < 0x1p-65L) + { + if (se == 0U && m >= 0x8000000000000000ULL) + { /* pseudo-denormal: glibc returns the (m, 1) encoding */ + u.p.se = 1; + return u.f; + } + return x; /* |x| < 2^-65: log1p(x) == x */ + } + yh = 1.0L + x; + yl = x - (yh - 1.0L); /* exact residual of the 1+x rounding for |x| <= 1 */ + t = yl / yh; + return vlibc_logl_core(yh) + t; +} + +long double +log1pl(long double x) +{ + if (x < -1.0L) + SET_EDOM(); + else if (x == -1.0L) + SET_ERANGE(); + return vlibc_log1pl_core(x); +} diff --git a/src/math/pow.c b/src/math/pow.c new file mode 100644 index 0000000..e815d8f --- /dev/null +++ b/src/math/pow.c @@ -0,0 +1,1347 @@ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifdef HAVE_CONFIG_H +#define SET_ERANGE() (errno = ERANGE) +#define SET_EDOM() (errno = EDOM) +#else +#define SET_ERANGE() \ + do \ + { \ + } while (0) +#define SET_EDOM() \ + do \ + { \ + } while (0) +#endif + +/* pow(x,y) = x**y, after FreeBSD msun e_pow.c (Sun fdlibm, 2004): + * log2(x) in two pieces, y*log2(x) = n+y' by simulated multi-precision + * arithmetic, then x**y = 2**n * exp(y'*log2). */ +static const double bp[] = {1.0, 1.5}; +static const double dp_h[] = {0.0, 5.84962487220764160156e-01}; /* 0x3FE2B80340000000 */ +static const double dp_l[] = {0.0, 1.35003920212974897128e-08}; /* 0x3E4CFDEB43CFD006 */ +static const double zero = 0.0; +static const double half = 0.5; +static const double qrtr = 0.25; +static const double thrd = 3.3333333333333331e-01; /* 0x3fd5555555555555 */ +static const double one = 1.0; +static const double two = 2.0; +static const double two53 = 9007199254740992.0; /* 0x4340000000000000 */ +static const double huge = 1.0e300; +static const double tiny = 1.0e-300; +static const double L1 = 5.99999999999994648725e-01; /* 0x3FE3333333333303 */ +static const double L2 = 4.28571428578550184252e-01; /* 0x3FDB6DB6DB6FABFF */ +static const double L3 = 3.33333329818377432918e-01; /* 0x3FD55555518F264D */ +static const double L4 = 2.72728123808534006489e-01; /* 0x3FD17460A91D4101 */ +static const double L5 = 2.30660745775561754067e-01; /* 0x3FCD864A93C9DB65 */ +static const double L6 = 2.06975017800338417784e-01; /* 0x3FCA7E284A454EEF */ +static const double P1 = 1.66666666666666019037e-01; /* 0x3FC555555555553E */ +static const double P2 = -2.77777777770155933842e-03; /* 0xBF66C16C16BEBD93 */ +static const double P3 = 6.61375632143793436117e-05; /* 0x3F11566AAF25DE2C */ +static const double P4 = -1.65339022054652515390e-06; /* 0xBEBBBD41C5D26BF1 */ +static const double P5 = 4.13813679705723846039e-08; /* 0x3E66376972BEA4D0 */ +static const double lg2 = 6.93147180559945286227e-01; /* 0x3FE62E42FEFA39EF */ +static const double lg2_h = 6.93147182464599609375e-01; /* 0x3FE62E4300000000 */ +static const double lg2_l = -1.90465429995776804525e-09; /* 0xBE205C610CA86C39 */ +static const double ovt = 8.0085662595372944372e-17; /* -(1024-log2(ovfl+.5ulp)) */ +static const double cp = 9.61796693925975554329e-01; /* 0x3FEEC709DC3A03FD =2/(3ln2) */ +static const double cp_h = 9.61796700954437255859e-01; /* 0x3FEEC709E0000000 =(float)cp */ +static const double cp_l = -7.02846165095275826516e-09; /* 0xBE3E2FE0145B01F5 =tail of cp_h */ +static const double ivln2 = 1.44269504088896338700e+00; /* 0x3FF71547652B82FE =1/ln2 */ +static const double ivln2_h = 1.44269502162933349609e+00; /* 0x3FF7154760000000 =24b 1/ln2 */ +static const double ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0BF85DDF44 =1/ln2 tail */ + +static double +vlibc_pow_core(double x, double y) +{ + union + { + double f; + unsigned long long i; + } u; + double z; + double ax; + double z_h; + double z_l; + double p_h; + double p_l; + double y1; + double t1; + double t2; + double r; + double s; + double t; + double uu; + double v; + double w; + int i; + int j; + int k; + int yisint; + int n; + int hx; + int hy; + int ix; + int iy; + unsigned int lx; + unsigned int ly; + + u.f = x; + hx = (int)(u.i >> 32); + lx = (unsigned int)u.i; + u.f = y; + hy = (int)(u.i >> 32); + ly = (unsigned int)u.i; + ix = hx & 0x7fffffff; + iy = hy & 0x7fffffff; + + /* y==zero: x**0 = 1 */ + if ((iy | (int)ly) == 0) + return one; + + /* x==1: 1**y = 1, even if y is NaN */ + if (hx == 0x3ff00000 && lx == 0) + return one; + + /* y!=zero: result is NaN if either arg is NaN */ + if (ix > 0x7ff00000 || ((ix == 0x7ff00000) && (lx != 0)) || iy > 0x7ff00000 || + ((iy == 0x7ff00000) && (ly != 0))) + return __builtin_nan(""); + + /* determine if y is an odd int when x < 0 */ + yisint = 0; + if (hx < 0) + { + if (iy >= 0x43400000) + yisint = 2; /* even integer y */ + else if (iy >= 0x3ff00000) + { + k = (iy >> 20) - 0x3ff; /* exponent */ + if (k > 20) + { + j = (int)(ly >> (52 - k)); + if (((unsigned int)j << (52 - k)) == ly) + yisint = 2 - (j & 1); + } + else if (ly == 0) + { + j = iy >> (20 - k); + if ((j << (20 - k)) == iy) + yisint = 2 - (j & 1); + } + } + } + + /* special value of y */ + if (ly == 0) + { + if (iy == 0x7ff00000) + { /* y is +-inf */ + if (((ix - 0x3ff00000) | (int)lx) == 0) + return one; /* (-1)**+-inf is 1 */ + else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */ + return (hy >= 0) ? y : zero; + else /* (|x|<1)**-,+inf = inf,0 */ + return (hy < 0) ? -y : zero; + } + if (iy == 0x3ff00000) + { /* y is +-1 */ + if (hy < 0) + return one / x; + else + return x; + } + if (hy == 0x40000000) + return x * x; /* y is 2 */ + if (hy == 0x3fe00000) + { /* y is 0.5 */ + if (hx >= 0) /* x >= +0 */ + { + double sq; + __asm__("sqrtsd %1, %0" : "=x"(sq) : "x"(x)); + return sq; + } + } + } + + ax = __builtin_fabs(x); + /* special value of x */ + if (lx == 0) + { + if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) + { + z = ax; /* x is +-0,+-inf,+-1 */ + if (hy < 0) + z = one / z; /* z = (1/|x|) */ + if (hx < 0) + { + if (((ix - 0x3ff00000) | yisint) == 0) + { + z = (z - z) / (z - z); /* (-1)**non-int is NaN */ + } + else if (yisint == 1) + z = -z; /* (x<0)**odd = -(|x|**odd) */ + } + return z; + } + } + + n = (int)((unsigned int)hx >> 31) - 1; + + /* (x<0)**(non-int) is NaN */ + if ((n | yisint) == 0) + return (x - x) / (x - x); + + s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ + if ((n | (yisint - 1)) == 0) + s = -one; /* (-ve)**(odd int) */ + + /* |y| is huge */ + if (iy > 0x41e00000) + { /* if |y| > 2**31 */ + if (iy > 0x43f00000) + { /* if |y| > 2**64, must o/uflow */ + if (ix <= 0x3fefffff) + return (hy < 0) ? huge * huge : tiny * tiny; + if (ix >= 0x3ff00000) + return (hy > 0) ? huge * huge : tiny * tiny; + } + /* over/underflow if x is not close to one */ + if (ix < 0x3fefffff) + return (hy < 0) ? s * huge * huge : s * tiny * tiny; + if (ix > 0x3ff00000) + return (hy > 0) ? s * huge * huge : s * tiny * tiny; + /* now |1-x| is tiny <= 2**-20, suffice to compute + log(x) by x-x^2/2+x^3/3-x^4/4 */ + t = ax - one; /* t has 20 trailing zeros */ + w = (t * t) * (half - t * (thrd - t * qrtr)); + uu = ivln2_h * t; /* ivln2_h has 21 sig. bits */ + v = t * ivln2_l - w * ivln2; + t1 = uu + v; + u.f = t1; + u.i &= 0xffffffff00000000ULL; + t1 = u.f; + t2 = v - (t1 - uu); + } + else + { + double ss; + double s2; + double s_h; + double s_l; + double t_h; + double t_l; + + n = 0; + /* take care subnormal number */ + if (ix < 0x00100000) + { + ax *= two53; + n -= 53; + u.f = ax; + ix = (int)(u.i >> 32); + } + n += (ix >> 20) - 0x3ff; + j = ix & 0x000fffff; + /* determine interval */ + ix = j | 0x3ff00000; /* normalize ix */ + if (j <= 0x3988E) + k = 0; /* |x|> 1) | 0x20000000) + 0x00080000 + (k << 18)) + << 32); + t_h = u.f; + t_l = ax - (t_h - bp[k]); + s_l = v * ((uu - s_h * t_h) - s_h * t_l); + /* compute log(ax) */ + s2 = ss * ss; + r = s2 * s2 * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6))))); + r += s_l * (s_h + ss); + s2 = s_h * s_h; + t_h = 3 + s2 + r; + u.f = t_h; + u.i &= 0xffffffff00000000ULL; + t_h = u.f; + t_l = r - ((t_h - 3) - s2); + /* uu+v = ss*(1+...) */ + uu = s_h * t_h; + v = s_l * t_h + t_l * ss; + /* 2/(3log2)*(ss+...) */ + p_h = uu + v; + u.f = p_h; + u.i &= 0xffffffff00000000ULL; + p_h = u.f; + p_l = v - (p_h - uu); + z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */ + z_l = cp_l * p_h + p_l * cp + dp_l[k]; + /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */ + t = n; + t1 = (((z_h + z_l) + dp_h[k]) + t); + u.f = t1; + u.i &= 0xffffffff00000000ULL; + t1 = u.f; + t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); + } + + /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ + y1 = y; + u.f = y1; + u.i &= 0xffffffff00000000ULL; + y1 = u.f; + p_l = (y - y1) * t1 + y * t2; + p_h = y1 * t1; + z = p_l + p_h; + u.f = z; + j = (int)(u.i >> 32); + i = (int)u.i; + if (j >= 0x40900000) + { /* z >= 1024 */ + if (((j - 0x40900000) | i) != 0) /* if z > 1024 */ + return s * huge * huge; /* overflow */ + else + { + if (p_l + ovt > z - p_h) + return s * huge * huge; /* overflow */ + } + } + else if ((j & 0x7fffffff) >= 0x4090cc00) + { /* z <= -1075 */ + if (((j - (int)0xc090cc00) | i) != 0) /* z < -1075 */ + return s * tiny * tiny; /* underflow */ + else + { + if (p_l <= z - p_h) + return s * tiny * tiny; /* underflow */ + } + } + /* compute 2**(p_h+p_l) */ + i = j & 0x7fffffff; + k = (i >> 20) - 0x3ff; + n = 0; + if (i > 0x3fe00000) + { /* if |z| > 0.5, set n = [z+0.5] */ + n = j + (0x00100000 >> (k + 1)); + k = ((n & 0x7fffffff) >> 20) - 0x3ff; /* new k for n */ + t = zero; + u.f = t; + u.i = (u.i & 0xffffffffULL) | + ((unsigned long long)(unsigned int)(n & ~(0x000fffff >> k)) << 32); + t = u.f; + n = ((n & 0x000fffff) | 0x00100000) >> (20 - k); + if (j < 0) + n = -n; + p_h -= t; + } + t = p_l + p_h; + u.f = t; + u.i &= 0xffffffff00000000ULL; + t = u.f; + uu = t * lg2_h; + v = (p_l - (t - p_h)) * lg2 + t * lg2_l; + z = uu + v; + w = v - (z - uu); + t = z * z; + t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); + r = (z * t1) / (t1 - two) - (w + z * w); + z = one - (r - z); + u.f = z; + j = (int)(u.i >> 32); + j += n << 20; + if ((j >> 20) <= 0) + { /* subnormal output: scale z by 2**n in two exact steps */ + union + { + double f; + unsigned long long i; + } b; + b.f = z; + if (n >= -1021) + { + b.i += (unsigned long long)n << 52; + z = b.f; + } + else + { + b.i += (unsigned long long)(n + 1022) << 52; + return s * b.f * 0x1p-1022; + } + } + else + { + u.f = z; + u.i = (u.i & 0xffffffffULL) | ((unsigned long long)(unsigned int)j << 32); + z = u.f; + } + return s * z; +} + +double +pow(double x, double y) +{ + double z = vlibc_pow_core(x, y); + + /* errno parity with glibc's wrapper: domain, pole/overflow, underflow */ + if (__builtin_isnan(z) && x < 0.0 && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_EDOM(); + else if (__builtin_isinf(z) && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_ERANGE(); + else if (z == 0.0 && x != 0.0 && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_ERANGE(); + return z; +} + +/* powf(x,y) = x**y, after FreeBSD msun e_powf.c (Sun fdlibm, Ian Lance + * Taylor's float conversion). */ +static const float bpf[] = {1.0, 1.5}; +static const float dp_hf[] = {0.0, 5.84960938e-01}; /* 0x3f15c000 */ +static const float dp_lf[] = {0.0, 1.56322085e-06}; /* 0x35d1cfdc */ +static const float zerof = 0.0; +static const float halff = 0.5; +static const float qrtrf = 0.25; +static const float thrdf = 3.33333343e-01; /* 0x3eaaaaab */ +static const float onef = 1.0; +static const float twof = 2.0; +static const float two24 = 16777216.0; /* 0x4b800000 */ +static const float hugef = 1.0e30; +static const float tinyf = 1.0e-30; +static const float L1f = 6.0000002384e-01; /* 0x3f19999a */ +static const float L2f = 4.2857143283e-01; /* 0x3edb6db7 */ +static const float L3f = 3.3333334327e-01; /* 0x3eaaaaab */ +static const float L4f = 2.7272811532e-01; /* 0x3e8ba305 */ +static const float L5f = 2.3066075146e-01; /* 0x3e6c3255 */ +static const float L6f = 2.0697501302e-01; /* 0x3e53f142 */ +static const float P1f = 1.6666667163e-01; /* 0x3e2aaaab */ +static const float P2f = -2.7777778450e-03; /* 0xbb360b61 */ +static const float P3f = 6.6137559770e-05; /* 0x388ab355 */ +static const float P4f = -1.6533901999e-06; /* 0xb5ddea0e */ +static const float P5f = 4.1381369442e-08; /* 0x3331bb4c */ +static const float lg2f = 6.9314718246e-01; /* 0x3f317218 */ +static const float lg2_hf = 6.93145752e-01; /* 0x3f317200 */ +static const float lg2_lf = 1.42860654e-06; /* 0x35bfbe8c */ +static const float ovtf = 4.2995665694e-08; /* -(128-log2(ovfl+.5ulp)) */ +static const float cpf = 9.6179670095e-01; /* 0x3f76384f =2/(3ln2) */ +static const float cp_hf = 9.6191406250e-01; /* 0x3f764000 =12b cp */ +static const float cp_lf = -1.1736857402e-04; /* 0xb8f623c6 =tail of cp_h */ +static const float ivln2f = 1.4426950216e+00; /* 0x3fb8aa3b =1/ln2 */ +static const float ivln2_hf = 1.4426879883e+00; /* 0x3fb8aa00 =16b 1/ln2 */ +static const float ivln2_lf = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail */ + +static float +vlibc_powf_core(float x, float y) +{ + union + { + float f; + unsigned int i; + } u; + float z; + float ax; + float z_h; + float z_l; + float p_h; + float p_l; + float y1; + float t1; + float t2; + float r; + float s; + float sn; + float t; + float uu; + float v; + float w; + int i; + int j; + int k; + int yisint; + int n; + int hx; + int hy; + int ix; + int iy; + int is; + + u.f = x; + hx = (int)u.i; + u.f = y; + hy = (int)u.i; + ix = hx & 0x7fffffff; + iy = hy & 0x7fffffff; + + /* y==zero: x**0 = 1 */ + if (iy == 0) + return onef; + + /* x==1: 1**y = 1, even if y is NaN */ + if (hx == 0x3f800000) + return onef; + + /* y!=zero: result is NaN if either arg is NaN */ + if (ix > 0x7f800000 || iy > 0x7f800000) + return __builtin_nanf(""); + + /* determine if y is an odd int when x < 0 */ + yisint = 0; + if (hx < 0) + { + if (iy >= 0x4b800000) + yisint = 2; /* even integer y */ + else if (iy >= 0x3f800000) + { + k = (iy >> 23) - 0x7f; /* exponent */ + j = iy >> (23 - k); + if ((j << (23 - k)) == iy) + yisint = 2 - (j & 1); + } + } + + /* special value of y */ + if (iy == 0x7f800000) + { /* y is +-inf */ + if (ix == 0x3f800000) + return onef; /* (-1)**+-inf is 1 */ + else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */ + return (hy >= 0) ? y : zerof; + else /* (|x|<1)**-,+inf = inf,0 */ + return (hy < 0) ? -y : zerof; + } + if (iy == 0x3f800000) + { /* y is +-1 */ + if (hy < 0) + return onef / x; + else + return x; + } + if (hy == 0x40000000) + return x * x; /* y is 2 */ + if (hy == 0x3f000000) + { /* y is 0.5 */ + if (hx >= 0) /* x >= +0 */ + { + float sq; + __asm__("sqrtss %1, %0" : "=x"(sq) : "x"(x)); + return sq; + } + } + + ax = __builtin_fabsf(x); + /* special value of x */ + if (ix == 0x7f800000 || ix == 0 || ix == 0x3f800000) + { + z = ax; /* x is +-0,+-inf,+-1 */ + if (hy < 0) + z = onef / z; /* z = (1/|x|) */ + if (hx < 0) + { + if (((ix - 0x3f800000) | yisint) == 0) + { + z = (z - z) / (z - z); /* (-1)**non-int is NaN */ + } + else if (yisint == 1) + z = -z; /* (x<0)**odd = -(|x|**odd) */ + } + return z; + } + + n = (int)((unsigned int)hx >> 31) - 1; + + /* (x<0)**(non-int) is NaN */ + if ((n | yisint) == 0) + return (x - x) / (x - x); + + sn = onef; /* s (sign of result -ve**odd) = -1 else = 1 */ + if ((n | (yisint - 1)) == 0) + sn = -onef; /* (-ve)**(odd int) */ + + /* |y| is huge */ + if (iy > 0x4d000000) + { /* if |y| > 2**27 */ + /* over/underflow if x is not close to one */ + if (ix < 0x3f7ffff6) + return (hy < 0) ? sn * hugef * hugef : sn * tinyf * tinyf; + if (ix > 0x3f800007) + return (hy > 0) ? sn * hugef * hugef : sn * tinyf * tinyf; + /* now |1-x| is tiny <= 2**-20, suffice to compute + log(x) by x-x^2/2+x^3/3-x^4/4 */ + t = ax - 1; /* t has 20 trailing zeros */ + w = (t * t) * (halff - t * (thrdf - t * qrtrf)); + uu = ivln2_hf * t; /* ivln2_h has 16 sig. bits */ + v = t * ivln2_lf - w * ivln2f; + t1 = uu + v; + u.f = t1; + u.i &= 0xfffff000U; + t1 = u.f; + t2 = v - (t1 - uu); + } + else + { + float s2; + float s_h; + float s_l; + float t_h; + float t_l; + + n = 0; + /* take care subnormal number (may need more than one 2^24 step) */ + if (ix < 0x00800000) + { + do + { + ax *= two24; + n -= 24; + u.f = ax; + ix = (int)u.i; + } while (ix < 0x00800000 && ix != 0); + } + n += (ix >> 23) - 0x7f; + j = ix & 0x007fffff; + /* determine interval */ + ix = j | 0x3f800000; /* normalize ix */ + if (j <= 0x1cc471) + k = 0; /* |x|> 1) & 0xfffff000) | 0x20000000; + u.f = t_h; + u.i = (unsigned int)(is + 0x00400000 + (k << 21)); + t_h = u.f; + t_l = ax - (t_h - bpf[k]); + s_l = v * ((uu - s_h * t_h) - s_h * t_l); + /* compute log(ax) */ + s2 = s * s; + r = s2 * s2 * (L1f + s2 * (L2f + s2 * (L3f + s2 * (L4f + s2 * (L5f + s2 * L6f))))); + r += s_l * (s_h + s); + s2 = s_h * s_h; + t_h = 3 + s2 + r; + u.f = t_h; + u.i &= 0xfffff000U; + t_h = u.f; + t_l = r - ((t_h - 3) - s2); + /* uu+v = s*(1+...) */ + uu = s_h * t_h; + v = s_l * t_h + t_l * s; + /* 2/(3log2)*(s+...) */ + p_h = uu + v; + u.f = p_h; + u.i &= 0xfffff000U; + p_h = u.f; + p_l = v - (p_h - uu); + z_h = cp_hf * p_h; /* cp_h+cp_l = 2/(3*log2) */ + z_l = cp_lf * p_h + p_l * cpf + dp_lf[k]; + /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ + t = n; + t1 = (((z_h + z_l) + dp_hf[k]) + t); + u.f = t1; + u.i &= 0xfffff000U; + t1 = u.f; + t2 = z_l - (((t1 - t) - dp_hf[k]) - z_h); + } + + /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ + u.f = y; + is = (int)u.i; + u.i &= 0xfffff000U; + y1 = u.f; + p_l = (y - y1) * t1 + y * t2; + p_h = y1 * t1; + z = p_l + p_h; + u.f = z; + j = (int)u.i; + if (j > 0x43000000) /* if z > 128 */ + return sn * hugef * hugef; /* overflow */ + else if (j == 0x43000000) + { /* if z == 128 */ + if (p_l + ovtf > z - p_h) + return sn * hugef * hugef; /* overflow */ + } + else if ((j & 0x7fffffff) > 0x43160000) /* z <= -150 */ + return sn * tinyf * tinyf; /* underflow */ + else if (j == (int)0xc3160000) + { /* z == -150 */ + if (p_l <= z - p_h) + return sn * tinyf * tinyf; /* underflow */ + } + /* compute 2**(p_h+p_l) */ + i = j & 0x7fffffff; + k = (i >> 23) - 0x7f; + n = 0; + if (i > 0x3f000000) + { /* if |z| > 0.5, set n = [z+0.5] */ + n = j + (0x00800000 >> (k + 1)); + k = ((n & 0x7fffffff) >> 23) - 0x7f; /* new k for n */ + u.f = t; + u.i = (unsigned int)(n & ~(0x007fffff >> k)); + t = u.f; + n = ((n & 0x007fffff) | 0x00800000) >> (23 - k); + if (j < 0) + n = -n; + p_h -= t; + } + t = p_l + p_h; + u.f = t; + u.i &= 0xffff8000U; + t = u.f; + uu = t * lg2_hf; + v = (p_l - (t - p_h)) * lg2f + t * lg2_lf; + z = uu + v; + w = v - (z - uu); + t = z * z; + t1 = z - t * (P1f + t * (P2f + t * (P3f + t * (P4f + t * P5f)))); + r = (z * t1) / (t1 - twof) - (w + z * w); + z = onef - (r - z); + u.f = z; + j = (int)u.i; + j += n << 23; + if ((j >> 23) <= 0) + { /* subnormal output: scale z by 2**n in two exact steps */ + union + { + float f; + unsigned int i; + } b; + b.f = z; + if (n >= -125) + { + b.i += (unsigned int)n << 23; + z = b.f; + } + else + { + b.i += (unsigned int)(n + 126) << 23; + return sn * b.f * 0x1p-126f; + } + } + else + { + u.f = z; + u.i = (unsigned int)j; + z = u.f; + } + return sn * z; +} + +float +powf(float x, float y) +{ + float z = vlibc_powf_core(x, y); + + if (__builtin_isnan(z) && x < 0.0f && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_EDOM(); + else if (__builtin_isinf(z) && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_ERANGE(); + else if (z == 0.0f && x != 0.0f && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_ERANGE(); + return z; +} + +/* powl(x,y) = x**y, after musl src/math/powl.c (OpenBSD ld80 e_powl.c, + * Stephen L. Moshier, 2008): a 2^-i/32 antilog table with pseudo extended + * precision for the logarithm and the exponential. */ +#define NXT 32 + +/* log(1+x) = x - .5x^2 + x^3 P(z)/Q(z) on 2^(-1/32)-1 <= x <= 2^(1/32)-1 */ +static const long double Pl[] = { + 8.3319510773868690346226E-4L, + 4.9000050881978028599627E-1L, + 1.7500123722550302671919E0L, + 1.4000100839971580279335E0L, +}; +static const long double Ql[] = { + 5.2500282295834889175431E0L, + 8.4000598057587009834666E0L, + 4.2000302519914740834728E0L, +}; +/* A[i] = 2^(-i/32) rounded to ld80; A[i] + B[i/2] adds accuracy for even i. */ +static const long double Al[33] = { + 1.0000000000000000000000E0L, 9.7857206208770013448287E-1L, 9.5760328069857364691013E-1L, + 9.3708381705514995065011E-1L, 9.1700404320467123175367E-1L, 8.9735453750155359320742E-1L, + 8.7812608018664974155474E-1L, 8.5930964906123895780165E-1L, 8.4089641525371454301892E-1L, + 8.2287773907698242225554E-1L, 8.0524516597462715409607E-1L, 7.8799042255394324325455E-1L, + 7.7110541270397041179298E-1L, 7.5458221379671136985669E-1L, 7.3841307296974965571198E-1L, + 7.2259040348852331001267E-1L, 7.0710678118654752438189E-1L, 6.9195494098191597746178E-1L, + 6.7712777346844636413344E-1L, 6.6261832157987064729696E-1L, 6.4841977732550483296079E-1L, + 6.3452547859586661129850E-1L, 6.2092890603674202431705E-1L, 6.0762367999023443907803E-1L, + 5.9460355750136053334378E-1L, 5.8186242938878875689693E-1L, 5.6939431737834582684856E-1L, + 5.5719337129794626814472E-1L, 5.4525386633262882960438E-1L, 5.3357020033841180906486E-1L, + 5.2213689121370692017331E-1L, 5.1094857432705833910408E-1L, 5.0000000000000000000000E-1L, +}; +static const long double Bl[17] = { + 0.0000000000000000000000E0L, 2.6176170809902549338711E-20L, -1.0126791927256478897086E-20L, + 1.3438228172316276937655E-21L, 1.2207982955417546912101E-20L, -6.3084814358060867200133E-21L, + 1.3164426894366316434230E-20L, -1.8527916071632873716786E-20L, 1.8950325588932570796551E-20L, + 1.5564775779538780478155E-20L, 6.0859793637556860974380E-21L, -2.0208749253662532228949E-20L, + 1.4966292219224761844552E-20L, 3.3540909728056476875639E-21L, -8.6987564101742849540743E-22L, + -1.2327176863327626135542E-20L, 0.0000000000000000000000E0L, +}; +/* log2(e) - 1, split hi/lo so the log's base-2 conversion keeps full + * precision (a single ld80 constant rounds at 2^-65) */ +static const long double LOG2EA_hi = + 0.44269504088896340736054672848620583636147785000503063201904296875L; +static const long double LOG2EA_lo = + -6.2204748431369893483189585204469788359356181889078081881492e-22L; + +static const long double MAXLOGL = 1.1356523406294143949492E4L; +static const long double MINLOGL = -1.13994985314888605586758E4L; +static const long double LOGE2L = 6.9314718055994530941723E-1L; +static const long double hugel = 0x1p10000L; +static const long double twom10000l = 0x1p-10000L; + +/* floor(x) for long double, via union bit surgery */ +static long double +vlibc_floorl(long double x) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } u; + int e; + + u.f = x; + e = (int)(u.p.se & 0x7fffU) - 16383; + if (e < 0) + { + if ((u.p.se & 0x8000U) != 0U && u.p.m != 0) + return -1.0L; /* negative fraction */ + return (u.p.m == 0 && (u.p.se & 0x8000U) != 0U) ? x : 0.0L; + } + if (e >= 63) + return x; + u.p.m &= ~((1ULL << (63 - e)) - 1); + if ((u.p.se & 0x8000U) != 0U && u.f != x) + u.f -= 1.0L; + return u.f; +} + +/* z * 2^i for any i the exponent arithmetic can produce, exact (power of 2) */ +static long double +vlibc_scalbnl(long double z, int i) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } two; + + if (i >= -16382 && i <= 16383) + { + two.p.m = 0x8000000000000000ULL; + two.p.se = (unsigned short)(16383 + i); + return z * two.f; + } + if (i > 16383) + { + two.p.m = 0x8000000000000000ULL; + two.p.se = 0x7ffe; /* 2^16383 */ + return z * two.f * (i - 16383 == 1 ? 2.0L : 4.0L); + } + /* i < -16382: 2^i = 2^(i+16445) * 2^-16445 (min subnormal) */ + two.p.m = 1; + two.p.se = 0; + return z * two.f * (i + 16445 == 1 ? 2.0L : i + 16445 == 2 ? 4.0L : 0.5L); +} + +/* x > 0 raised to an integer power (binary exponentiation) */ +static long double +vlibc_powil(long double x, int nn) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } u; + long double ww; + long double y; + long double s; + int n; + int e; + int sign; + int lx; + + if (nn == 0) + return 1.0L; + + if (nn < 0) + { + sign = -1; + n = -nn; + } + else + { + sign = 1; + n = nn; + } + + /* overflow detection: approximate logarithm of the answer */ + s = x; + u.f = s; + if (u.p.se == 0) + { /* subnormal: normalize mantissa first (frexpl semantics) */ + int sh = __builtin_clzll(u.p.m); + u.p.m <<= sh; + lx = -16381 - sh; + } + else + { + lx = (int)(u.p.se & 0x7fffU) - 16383 + 1; + u.p.se = 0x3ffe; + } + s = u.f; + e = (lx - 1) * n; + if ((e == 0) || (e > 64) || (e < -64)) + { + s = (s - 7.0710678118654752e-1L) / (s + 7.0710678118654752e-1L); + s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L; + } + else + { + s = LOGE2L * e; + } + + if (s > MAXLOGL) + return hugel * hugel; /* overflow */ + if (s < MINLOGL) + return twom10000l * twom10000l; /* underflow */ + if (s < -MAXLOGL + 2.0L) + { + x = 1.0L / x; + sign = -sign; + } + + /* first bit of the power */ + if (n & 1) + y = x; + else + y = 1.0L; + + ww = x; + n >>= 1; + while (n) + { + ww = ww * ww; /* arg to the 2-to-the-kth power */ + if (n & 1) /* if that bit is set, then include in product */ + y *= ww; + n >>= 1; + } + + if (sign < 0) + y = 1.0L / y; + return y; +} + +static long double +vlibc_powl_core(long double x, long double y) +{ + union + { + long double f; + struct + { + unsigned long long m; + unsigned short se; + } p; + } u; + int i; + int nflg; + int iyflg; + int yoddint; + long e; + long double z; + long double zlo; + long double xlo; + long double w; + long double ya; + long double yb; + + /* make sure no invalid exception is raised by nan comparison */ + if (__builtin_isnan(x)) + { + if (!__builtin_isnan(y) && y == 0.0L) + return 1.0L; + return x; + } + if (__builtin_isnan(y)) + { + if (x == 1.0L) + return 1.0L; + return y; + } + if (x == 1.0L) + return 1.0L; /* 1**y = 1, even if y is nan */ + if (y == 0.0L) + return 1.0L; /* x**0 = 1, even if x is nan */ + if (y == 1.0L) + return x; + if (__builtin_fabsl(y) > 2 * (-LDBL_MIN_EXP + LDBL_MANT_DIG + 1) / LDBL_EPSILON) + { + /* y is not an odd int */ + if (x == -1.0L) + return 1.0L; + if (y == __builtin_huge_vall()) + { + if (x > 1.0L || x < -1.0L) + return __builtin_huge_vall(); + return 0.0L; + } + if (y == -__builtin_huge_vall()) + { + if (x > 1.0L || x < -1.0L) + return 0.0L; + return __builtin_huge_vall(); + } + if ((x > 1.0L || x < -1.0L) == (y > 0.0L)) + return hugel * hugel; + return twom10000l * twom10000l; + } + if (x == __builtin_huge_vall()) + { + if (y > 0.0L) + return __builtin_huge_vall(); + return 0.0L; + } + + w = vlibc_floorl(y); + + /* set iyflg to 1 if y is an integer */ + iyflg = 0; + if (w == y) + iyflg = 1; + + /* test for odd integer y */ + yoddint = 0; + if (iyflg) + { + ya = __builtin_fabsl(y); + ya = vlibc_floorl(0.5L * ya); + yb = 0.5L * __builtin_fabsl(w); + if (ya != yb) + yoddint = 1; + } + + if (x == -__builtin_huge_vall()) + { + if (y > 0.0L) + { + if (yoddint) + return -__builtin_huge_vall(); + return __builtin_huge_vall(); + } + if (y < 0.0L) + { + if (yoddint) + return -0.0L; + return 0.0L; + } + } + nflg = 0; /* (x<0)**(odd int) */ + if (x <= 0.0L) + { + if (x == 0.0L) + { + u.f = x; + if (y < 0.0L) + { + if ((u.p.se & 0x8000U) != 0U && yoddint) + return -__builtin_huge_vall(); /* (-0.0)**(-odd int) */ + return __builtin_huge_vall(); /* (+-0.0)**(negative) */ + } + if ((u.p.se & 0x8000U) != 0U && yoddint) + return -0.0L; + return 0.0L; + } + if (iyflg == 0) + return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */ + /* (x<0)**(integer) */ + if (yoddint) + nflg = 1; /* negate result */ + x = -x; + } + /* (+integer)**(integer) */ + if (iyflg && vlibc_floorl(x) == x && __builtin_fabsl(y) < 32768.0L) + { + w = vlibc_powil(x, (int)y); + return nflg ? -w : w; + } + + /* separate significand from exponent: x = frac * 2^i, frac in [0.5,1) */ + u.f = x; + if (u.p.se == 0) + { /* subnormal: normalize mantissa first (frexpl semantics) */ + int sh = __builtin_clzll(u.p.m); + u.p.m <<= sh; + e = (long)(-16381 - sh); + } + else + { + e = (long)(u.p.se & 0x7fffU) - 16383 + 1; + u.p.se = 0x3ffe; + } + x = u.f; + i = (int)e; + + /* find significand in antilog table A[] */ + i = 1; + if (x <= Al[17]) + i = 17; + if (x <= Al[i + 8]) + i += 8; + if (x <= Al[i + 4]) + i += 4; + if (x <= Al[i + 2]) + i += 2; + if (x >= Al[1]) + i = -1; + i += 1; + + /* find v = (x - A[i] - B[i/2]) / A[i] in double-double to compute + * log(x/A[i]) = log(1+v): the naive division would round at 2^-64 of + * v and cost ~|y| ulps of the result */ + x -= Al[i]; + x -= Bl[i / 2]; + { + const long double split = 0x1p32L + 1.0L; + long double q = x / Al[i]; + long double qh = q * split; + long double ql; + long double ah = Al[i] * split; + long double al; + long double ph; + long double pl; + long double r; + + qh = qh - (qh - q); + ql = q - qh; + ah = ah - (ah - Al[i]); + al = Al[i] - ah; + ph = qh * ah; /* 32x32: exact */ + pl = qh * al + ql * ah + ql * al; + r = (x - ph) - pl; + x = q; + xlo = r / Al[i]; + } + + /* rational approximation for log(1+v) = v - v**2/2 + v**3 P(v)/Q(v), + * plus the first-order correction for the v_lo tail */ + z = x * x; + { + long double pv = Pl[0]; + long double qv = x + Ql[0]; + int m; + + for (m = 1; m <= 3; m++) + pv = pv * x + Pl[m]; + for (m = 1; m < 3; m++) + qv = qv * x + Ql[m]; + w = x * (z * pv / qv); + } + w = w - 0.5L * z; + w += xlo * (1.0L - x + z); /* v_lo * d/dv log(1+v) */ + + /* convert to base 2 logarithm: multiply by log2(e) = 1 + LOG2EA, + * accumulating the fraction as a double-double so the log keeps full + * precision for any |x| (a plain sum would lose ~2^-66 absolute) */ + { + long double t1 = x; + long double zl = 0.0L; + long double t; + long double s; + long double zz; + long double e; + + t = LOG2EA_hi * x; + s = t1 + t; + zz = s - t1; + e = (t1 - (s - zz)) + (t - zz); + t1 = s; + zl += e; + s = t1 + w; + zz = s - t1; + e = (t1 - (s - zz)) + (w - zz); + t1 = s; + zl += e; + t = LOG2EA_hi * w; + s = t1 + t; + zz = s - t1; + e = (t1 - (s - zz)) + (t - zz); + t1 = s; + zl += e; + zl += (x + w) * LOG2EA_lo; /* hi/lo correction of the constant */ + z = t1; + zlo = zl; + } + + /* compute exponent term of the base 2 logarithm */ + w = -i; + w /= NXT; + w += e; + /* now base 2 log of x is w + z + zlo */ + + /* t = y * (w + z + zlo) as a double-double via Veltkamp 32-bit + * splitting: the products of 32-bit halves are exact, so t keeps full + * precision even for huge |log2(x)| where the naive product loses ulps */ + { + const long double split = 0x1p32L + 1.0L; + long double yh = y * split; + long double yl; + long double zh; + long double zl; + long double p; + long double p2; + long double th; + long double tl; + long double sh; + long double sl; + long double s; + long double zz; + long double e2; + long double rr; + long double term; + long double exp_h; + long double exp_l; + int m; + + /* sh + sl = w + z + zlo: the log2(x) as a double-double */ + s = w + z; + zz = s - w; + e2 = (w - (s - zz)) + (z - zz); + sh = s; + sl = e2 + zlo; + + yh = yh - (yh - y); + yl = y - yh; + + /* th = fl(y*sh); tl keeps the exact residual of the product */ + th = yh * sh; + zh = sh * split; + zh = zh - (zh - sh); + zl = sh - zh; + p = yh * zh; /* 32x32: exact */ + p2 = yh * zl; /* 32x32: exact */ + s = p + p2; + zz = s - p; + tl = (p - (s - zz)) + (p2 - zz); + tl += yl * sh; + tl += y * sl; + + /* overflow/underflow before the integer conversion */ + if (th > 16384.0L) + return nflg ? -hugel * hugel : hugel * hugel; + if (th < -16445.0L) + return nflg ? -0.0L : 0.0L; + + /* N = round(t), r = t - N in (-0.5, 0.5]; when the fractional part + * lives entirely in tl (integer y, |y*z| > 0.5) a single +-1 patch + * is not enough, so loop until the remainder is in range */ + m = (int)(th >= 0.0L ? th + 0.5L : th - 0.5L); + rr = (th - (long double)m) + tl; + while (rr > 0.5L) + { + rr -= 1.0L; + m += 1; + } + while (rr <= -0.5L) + { + rr += 1.0L; + m -= 1; + } + + /* 2^t = 2^m * e^(rr*ln2) via a TwoSum-accumulated Taylor sum */ + rr *= 0x1.62e42fefa39ef35793c7673007e5edp-1L; /* ln2 */ + term = rr; + s = 1.0L + rr; + zz = s - 1.0L; + exp_h = s; + exp_l = (1.0L - (s - zz)) + (rr - zz); + for (i = 2; i <= 32; i++) + { + term *= rr / (long double)i; + if (term > -1e-30L && term < 1e-30L) + break; + s = exp_h + term; + zz = s - exp_h; + e2 = (exp_h - (s - zz)) + (term - zz); + exp_h = s; + exp_l += e2; + } + z = vlibc_scalbnl(exp_h + exp_l, m); + } + + if (nflg) + z = -z; + return z; +} + +long double +powl(long double x, long double y) +{ + long double z = vlibc_powl_core(x, y); + + if (__builtin_isnan(z) && x < 0.0L && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_EDOM(); + else if (__builtin_isinf(z) && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_ERANGE(); + else if (z == 0.0L && x != 0.0L && __builtin_isfinite(x) && __builtin_isfinite(y)) + SET_ERANGE(); + return z; +} diff --git a/src/math/sqrt.c b/src/math/sqrt.c new file mode 100644 index 0000000..a608c2e --- /dev/null +++ b/src/math/sqrt.c @@ -0,0 +1,56 @@ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifdef HAVE_CONFIG_H +#define SET_ERANGE() (errno = ERANGE) +#define SET_EDOM() (errno = EDOM) +#else +#define SET_ERANGE() \ + do \ + { \ + } while (0) +#define SET_EDOM() \ + do \ + { \ + } while (0) +#endif + +/* sqrt is the one function IEEE 754 requires to be correctly rounded. The + * sqrtsd/sqrtss/fsqrt hardware instructions are; inline asm keeps the call + * self-contained (GCC's __builtin_sqrt would emit a libm call for NaN + * inputs and a false recursion warning inside a function named sqrt). */ +double +sqrt(double x) +{ + double y; + + if (x < 0.0) + SET_EDOM(); + __asm__("sqrtsd %1, %0" : "=x"(y) : "x"(x)); + return y; +} + +float +sqrtf(float x) +{ + float y; + + if (x < 0.0f) + SET_EDOM(); + __asm__("sqrtss %1, %0" : "=x"(y) : "x"(x)); + return y; +} + +long double +sqrtl(long double x) +{ + if (x < 0.0L) + SET_EDOM(); + __asm__("fsqrt" : "+t"(x)); + return x; +}