feat(math): exp/log/pow/sqrt/cbrt/hypot/log1p/expm1
This commit is contained in:
+236
@@ -0,0 +1,236 @@
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <config.h>
|
||||
#endif
|
||||
#include <math.h>
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <errno.h>
|
||||
#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;
|
||||
}
|
||||
+661
@@ -0,0 +1,661 @@
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <config.h>
|
||||
#endif
|
||||
#include <math.h>
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <errno.h>
|
||||
#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;
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,265 @@
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <config.h>
|
||||
#endif
|
||||
#include <math.h>
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <errno.h>
|
||||
#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;
|
||||
}
|
||||
+457
@@ -0,0 +1,457 @@
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <config.h>
|
||||
#endif
|
||||
#include <math.h>
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <errno.h>
|
||||
#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);
|
||||
}
|
||||
@@ -0,0 +1,427 @@
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <config.h>
|
||||
#endif
|
||||
#include <math.h>
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <errno.h>
|
||||
#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);
|
||||
}
|
||||
+1347
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,56 @@
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <config.h>
|
||||
#endif
|
||||
#include <math.h>
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include <errno.h>
|
||||
#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;
|
||||
}
|
||||
Reference in New Issue
Block a user