text:
#include <fenv.h>
#include "libm.h"

#if LDBL_MANT_DIG==64 && LDBL_MAX_EXP==16384
/* exact add, assumes exponent_x >= exponent_y */
static void add(long double *hi, long double *lo, long double x, long double y)
{
	long double r;

	r = x + y;
	*hi = r;
	r -= x;
	*lo = y - r;
}

/* exact mul, assumes no over/underflow */
static void mul(long double *hi, long double *lo, long double x, long double y)
{
	static const long double c = 1.0 + 0x1p32L;
	long double cx, xh, xl, cy, yh, yl;

	cx = c*x;
	xh = (x - cx) + cx;
	xl = x - xh;
	cy = c*y;
	yh = (y - cy) + cy;
	yl = y - yh;
	*hi = x*y;
	*lo = (xh*yh - *hi) + xh*yl + xl*yh + xl*yl;
}

/*
assume (long double)(hi+lo) == hi
return an adjusted hi so that rounding it to double (or less) precision is correct
*/
static long double adjust(long double hi, long double lo)
{
	union ldshape uhi, ulo;

	if (lo == 0)
		return hi;
	uhi.f = hi;
	if (uhi.i.m & 0x3ff)
		return hi;
	ulo.f = lo;
	if ((uhi.i.se & 0x8000) == (ulo.i.se & 0x8000))
		uhi.i.m++;
	else {
		/* handle underflow and take care of ld80 implicit msb */
		if (uhi.i.m << 1 == 0) {
			uhi.i.m = 0;
			uhi.i.se--;
		}
		uhi.i.m--;
	}
	return uhi.f;
}

/* adjusted add so the result is correct when rounded to double (or less) precision */
static long double dadd(long double x, long double y)
{
	add(&x, &y, x, y);
	return adjust(x, y);
}

/* adjusted mul so the result is correct when rounded to double (or less) precision */
static long double dmul(long double x, long double y)
{
	mul(&x, &y, x, y);
	return adjust(x, y);
}

static int getexp(long double x)
{
	union ldshape u;
	u.f = x;
	return u.i.se & 0x7fff;
}

double fma(double x, double y, double z)
{
	#pragma STDC FENV_ACCESS ON
	long double hi, lo1, lo2, xy;
	int round, ez, exy;

	/* handle +-inf,nan */
	if (!isfinite(x) || !isfinite(y))
		return x*y + z;
	if (!isfinite(z))
		return z;
	/* handle +-0 */
	if (x == 0.0 || y == 0.0)
		return x*y + z;
	round = fegetround();
	if (z == 0.0) {
		if (round == FE_TONEAREST)
			return dmul(x, y);
		return x*y;
	}

	/* exact mul and add require nearest rounding */
	/* spurious inexact exceptions may be raised */
	fesetround(FE_TONEAREST);
	mul(&xy, &lo1, x, y);
	exy = getexp(xy);
	ez = getexp(z);
	if (ez > exy) {
		add(&hi, &lo2, z, xy);
	} else if (ez > exy - 12) {
		add(&hi, &lo2, xy, z);
		if (hi == 0) {
			/*
			xy + z is 0, but it should be calculated with the
			original rounding mode so the sign is correct, if the
			compiler does not support FENV_ACCESS ON it does not
			know about the changed rounding mode and eliminates
			the xy + z below without the volatile memory access
			*/
			volatile double z_;
			fesetround(round);
			z_ = z;
			return (xy + z_) + lo1;
		}
	} else {
		/*
		ez <= exy - 12
		the 12 extra bits (1guard, 11round+sticky) are needed so with
			lo = dadd(lo1, lo2)
		elo <= ehi - 11, and we use the last 10 bits in adjust so
			dadd(hi, lo)
		gives correct result when rounded to double
		*/
		hi = xy;
		lo2 = z;
	}
	/*
	the result is stored before return for correct precision and exceptions

	one corner case is when the underflow flag should be raised because
	the precise result is an inexact subnormal double, but the calculated
	long double result is an exact subnormal double
	(so rounding to double does not raise exceptions)

	in nearest rounding mode dadd takes care of this: the last bit of the
	result is adjusted so rounding sees an inexact value when it should

	in non-nearest rounding mode fenv is used for the workaround
	*/
	fesetround(round);
	if (round == FE_TONEAREST)
		z = dadd(hi, dadd(lo1, lo2));
	else {
#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
		int e = fetestexcept(FE_INEXACT);
		feclearexcept(FE_INEXACT);
#endif
		z = hi + (lo1 + lo2);
#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
		if (getexp(z) < 0x3fff-1022 && fetestexcept(FE_INEXACT))
			feraiseexcept(FE_UNDERFLOW);
		else if (e)
			feraiseexcept(FE_INEXACT);
#endif
	}
	return z;
}
#else
/* origin: FreeBSD /usr/src/lib/msun/src/s_fma.c */
/*-
 * Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

/*
 * A struct dd represents a floating-point number with twice the precision
 * of a double.  We maintain the invariant that "hi" stores the 53 high-order
 * bits of the result.
 */
struct dd {
	double hi;
	double lo;
};

/*
 * Compute a+b exactly, returning the exact result in a struct dd.  We assume
 * that both a and b are finite, but make no assumptions about their relative
 * magnitudes.
 */
static inline struct dd dd_add(double a, double b)
{
	struct dd ret;
	double s;

	ret.hi = a + b;
	s = ret.hi - a;
	ret.lo = (a - (ret.hi - s)) + (b - s);
	return (ret);
}

/*
 * Compute a+b, with a small tweak:  The least significant bit of the
 * result is adjusted into a sticky bit summarizing all the bits that
 * were lost to rounding.  This adjustment negates the effects of double
 * rounding when the result is added to another number with a higher
 * exponent.  For an explanation of round and sticky bits, see any reference
 * on FPU design, e.g.,
 *
 *     J. Coonen.  An Implementation Guide to a Proposed Standard for
 *     Floating-Point Arithmetic.  Computer, vol. 13, no. 1, Jan 1980.
 */
static inline double add_adjusted(double a, double b)
{
	struct dd sum;
	union {double f; uint64_t i;} uhi, ulo;

	sum = dd_add(a, b);
	if (sum.lo != 0) {
		uhi.f = sum.hi;
		if ((uhi.i & 1) == 0) {
			/* hibits += (int)copysign(1.0, sum.hi * sum.lo) */
			ulo.f = sum.lo;
			uhi.i += 1 - ((uhi.i ^ ulo.i) >> 62);
			sum.hi = uhi.f;
		}
	}
	return (sum.hi);
}

/*
 * Compute ldexp(a+b, scale) with a single rounding error. It is assumed
 * that the result will be subnormal, and care is taken to ensure that
 * double rounding does not occur.
 */
static inline double add_and_denormalize(double a, double b, int scale)
{
	struct dd sum;
	union {double f; uint64_t i;} uhi, ulo;
	int bits_lost;

	sum = dd_add(a, b);

	/*
	 * If we are losing at least two bits of accuracy to denormalization,
	 * then the first lost bit becomes a round bit, and we adjust the
	 * lowest bit of sum.hi to make it a sticky bit summarizing all the
	 * bits in sum.lo. With the sticky bit adjusted, the hardware will
	 * break any ties in the correct direction.
	 *
	 * If we are losing only one bit to denormalization, however, we must
	 * break the ties manually.
	 */
	if (sum.lo != 0) {
		uhi.f = sum.hi;
		bits_lost = -((int)(uhi.i >> 52) & 0x7ff) - scale + 1;
		if (bits_lost != 1 ^ (int)(uhi.i & 1)) {
			/* hibits += (int)copysign(1.0, sum.hi * sum.lo) */
			ulo.f = sum.lo;
			uhi.i += 1 - (((uhi.i ^ ulo.i) >> 62) & 2);
			sum.hi = uhi.f;
		}
	}
	return scalbn(sum.hi, scale);
}

/*
 * Compute a*b exactly, returning the exact result in a struct dd.  We assume
 * that both a and b are normalized, so no underflow or overflow will occur.
 * The current rounding mode must be round-to-nearest.
 */
static inline struct dd dd_mul(double a, double b)
{
	static const double split = 0x1p27 + 1.0;
	struct dd ret;
	double ha, hb, la, lb, p, q;

	p = a * split;
	ha = a - p;
	ha += p;
	la = a - ha;

	p = b * split;
	hb = b - p;
	hb += p;
	lb = b - hb;

	p = ha * hb;
	q = ha * lb + la * hb;

	ret.hi = p + q;
	ret.lo = p - ret.hi + q + la * lb;
	return (ret);
}

/*
 * Fused multiply-add: Compute x * y + z with a single rounding error.
 *
 * We use scaling to avoid overflow/underflow, along with the
 * canonical precision-doubling technique adapted from:
 *
 *      Dekker, T.  A Floating-Point Technique for Extending the
 *      Available Precision.  Numer. Math. 18, 224-242 (1971).
 *
 * This algorithm is sensitive to the rounding precision.  FPUs such
 * as the i387 must be set in double-precision mode if variables are
 * to be stored in FP registers in order to avoid incorrect results.
 * This is the default on FreeBSD, but not on many other systems.
 *
 * Hardware instructions should be used on architectures that support it,
 * since this implementation will likely be several times slower.
 */
double fma(double x, double y, double z)
{
	#pragma STDC FENV_ACCESS ON
	double xs, ys, zs, adj;
	struct dd xy, r;
	int oround;
	int ex, ey, ez;
	int spread;

	/*
	 * Handle special cases. The order of operations and the particular
	 * return values here are crucial in handling special cases involving
	 * infinities, NaNs, overflows, and signed zeroes correctly.
	 */
	if (!isfinite(x) || !isfinite(y))
		return (x * y + z);
	if (!isfinite(z))
		return (z);
	if (x == 0.0 || y == 0.0)
		return (x * y + z);
	if (z == 0.0)
		return (x * y);

	xs = frexp(x, &ex);
	ys = frexp(y, &ey);
	zs = frexp(z, &ez);
	oround = fegetround();
	spread = ex + ey - ez;

	/*
	 * If x * y and z are many orders of magnitude apart, the scaling
	 * will overflow, so we handle these cases specially.  Rounding
	 * modes other than FE_TONEAREST are painful.
	 */
	if (spread < -DBL_MANT_DIG) {
#ifdef FE_INEXACT
		feraiseexcept(FE_INEXACT);
#endif
#ifdef FE_UNDERFLOW
		if (!isnormal(z))
			feraiseexcept(FE_UNDERFLOW);
#endif
		switch (oround) {
		default: /* FE_TONEAREST */
			return (z);
#ifdef FE_TOWARDZERO
		case FE_TOWARDZERO:
			if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
				return (z);
			else
				return (nextafter(z, 0));
#endif
#ifdef FE_DOWNWARD
		case FE_DOWNWARD:
			if (x > 0.0 ^ y < 0.0)
				return (z);
			else
				return (nextafter(z, -INFINITY));
#endif
#ifdef FE_UPWARD
		case FE_UPWARD:
			if (x > 0.0 ^ y < 0.0)
				return (nextafter(z, INFINITY));
			else
				return (z);
#endif
		}
	}
	if (spread <= DBL_MANT_DIG * 2)
		zs = scalbn(zs, -spread);
	else
		zs = copysign(DBL_MIN, zs);

	fesetround(FE_TONEAREST);

	/*
	 * Basic approach for round-to-nearest:
	 *
	 *     (xy.hi, xy.lo) = x * y           (exact)
	 *     (r.hi, r.lo)   = xy.hi + z       (exact)
	 *     adj = xy.lo + r.lo               (inexact; low bit is sticky)
	 *     result = r.hi + adj              (correctly rounded)
	 */
	xy = dd_mul(xs, ys);
	r = dd_add(xy.hi, zs);

	spread = ex + ey;

	if (r.hi == 0.0) {
		/*
		 * When the addends cancel to 0, ensure that the result has
		 * the correct sign.
		 */
		fesetround(oround);
		volatile double vzs = zs; /* XXX gcc CSE bug workaround */
		return xy.hi + vzs + scalbn(xy.lo, spread);
	}

	if (oround != FE_TONEAREST) {
		/*
		 * There is no need to worry about double rounding in directed
		 * rounding modes.
		 * But underflow may not be raised properly, example in downward rounding:
		 * fma(0x1.000000001p-1000, 0x1.000000001p-30, -0x1p-1066)
		 */
		double ret;
#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
		int e = fetestexcept(FE_INEXACT);
		feclearexcept(FE_INEXACT);
#endif
		fesetround(oround);
		adj = r.lo + xy.lo;
		ret = scalbn(r.hi + adj, spread);
#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
		if (ilogb(ret) < -1022 && fetestexcept(FE_INEXACT))
			feraiseexcept(FE_UNDERFLOW);
		else if (e)
			feraiseexcept(FE_INEXACT);
#endif
		return ret;
	}

	adj = add_adjusted(r.lo, xy.lo);
	if (spread + ilogb(r.hi) > -1023)
		return scalbn(r.hi + adj, spread);
	else
		return add_and_denormalize(r.hi, adj, spread);
}
#endif


text:
#include "libm.h"

float tanhf(float x)
{
	union {float f; uint32_t i;} u = {.f = x};
	uint32_t w;
	int sign;
	float t;

	/* x = |x| */
	sign = u.i >> 31;
	u.i &= 0x7fffffff;
	x = u.f;
	w = u.i;

	if (w > 0x3f0c9f54) {
		/* |x| > log(3)/2 ~= 0.5493 or nan */
		if (w > 0x41200000) {
			/* |x| > 10 */
			t = 1 + 0/x;
		} else {
			t = expm1f(2*x);
			t = 1 - 2/(t+2);
		}
	} else if (w > 0x3e82c578) {
		/* |x| > log(5/3)/2 ~= 0.2554 */
		t = expm1f(2*x);
		t = t/(t+2);
	} else if (w >= 0x00800000) {
		/* |x| >= 0x1p-126 */
		t = expm1f(-2*x);
		t = -t/(t+2);
	} else {
		/* |x| is subnormal */
		FORCE_EVAL(x*x);
		t = x;
	}
	return sign ? -t : t;
}

repo_language: C

mime_type: text/x-c


text:
#include <math.h>
#include <float.h>

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double fminl(long double x, long double y)
{
	return fmin(x, y);
}
#else
long double fminl(long double x, long double y)
{
	if (isnan(x))
		return y;
	if (isnan(y))
		return x;
	/* handle signed zeros, see C99 Annex F.9.9.2 */
	if (signbit(x) != signbit(y))
		return signbit(x) ? x : y;
	return x < y ? x : y;
}
#endif


text:
#include "libm.h"

/* cosh(x) = (exp(x) + 1/exp(x))/2
 *         = 1 + 0.5*(exp(x)-1)*(exp(x)-1)/exp(x)
 *         = 1 + x*x/2 + o(x^4)
 */
double cosh(double x)
{
	union {double f; uint64_t i;} u = {.f = x};
	uint32_t w;
	double t;

	/* |x| */
	u.i &= (uint64_t)-1/2;
	x = u.f;
	w = u.i >> 32;

	/* |x| < log(2) */
	if (w < 0x3fe62e42) {
		if (w < 0x3ff00000 - (26<<20)) {
			/* raise inexact if x!=0 */
			FORCE_EVAL(x + 0x1p120f);
			return 1;
		}
		t = expm1(x);
		return 1 + t*t/(2*(1+t));
	}

	/* |x| < log(DBL_MAX) */
	if (w < 0x40862e42) {
		t = exp(x);
		/* note: if x>log(0x1p26) then the 1/t is not needed */
		return 0.5*(t + 1/t);
	}

	/* |x| > log(DBL_MAX) or nan */
	/* note: the result is stored to handle overflow */
	t = __expo2(x);
	return t;
}


text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_j1.c */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */
/* j1(x), y1(x)
 * Bessel function of the first and second kinds of order zero.
 * Method -- j1(x):
 *      1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
 *      2. Reduce x to |x| since j1(x)=-j1(-x),  and
 *         for x in (0,2)
 *              j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
 *         (precision:  |j1/x - 1/2 - R0/S0 |<2**-61.51 )
 *         for x in (2,inf)
 *              j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
 *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
 *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
 *         as follow:
 *              cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
 *                      =  1/sqrt(2) * (sin(x) - cos(x))
 *              sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 *                      = -1/sqrt(2) * (sin(x) + cos(x))
 *         (To avoid cancellation, use
 *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 *          to compute the worse one.)
 *
 *      3 Special cases
 *              j1(nan)= nan
 *              j1(0) = 0
 *              j1(inf) = 0
 *
 * Method -- y1(x):
 *      1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
 *      2. For x<2.
 *         Since
 *              y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
 *         therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
 *         We use the following function to approximate y1,
 *              y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
 *         where for x in [0,2] (abs err less than 2**-65.89)
 *              U(z) = U0[0] + U0[1]*z + ... + U0[4]*z^4
 *              V(z) = 1  + v0[0]*z + ... + v0[4]*z^5
 *         Note: For tiny x, 1/x dominate y1 and hence
 *              y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
 *      3. For x>=2.
 *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
 *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
 *         by method mentioned above.
 */

#include "libm.h"

static double pone(double), qone(double);

static const double
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi       = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */

static double common(uint32_t ix, double x, int y1, int sign)
{
	double z,s,c,ss,cc;

	/*
	 * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x-3pi/4)-q1(x)*sin(x-3pi/4))
	 * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x-3pi/4)+q1(x)*cos(x-3pi/4))
	 *
	 * sin(x-3pi/4) = -(sin(x) + cos(x))/sqrt(2)
	 * cos(x-3pi/4) = (sin(x) - cos(x))/sqrt(2)
	 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
	 */
	s = sin(x);
	if (y1)
		s = -s;
	c = cos(x);
	cc = s-c;
	if (ix < 0x7fe00000) {
		/* avoid overflow in 2*x */
		ss = -s-c;
		z = cos(2*x);
		if (s*c > 0)
			cc = z/ss;
		else
			ss = z/cc;
		if (ix < 0x48000000) {
			if (y1)
				ss = -ss;
			cc = pone(x)*cc-qone(x)*ss;
		}
	}
	if (sign)
		cc = -cc;
	return invsqrtpi*cc/sqrt(x);
}

/* R0/S0 on [0,2] */
static const double
r00 = -6.25000000000000000000e-02, /* 0xBFB00000, 0x00000000 */
r01 =  1.40705666955189706048e-03, /* 0x3F570D9F, 0x98472C61 */
r02 = -1.59955631084035597520e-05, /* 0xBEF0C5C6, 0xBA169668 */
r03 =  4.96727999609584448412e-08, /* 0x3E6AAAFA, 0x46CA0BD9 */
s01 =  1.91537599538363460805e-02, /* 0x3F939D0B, 0x12637E53 */
s02 =  1.85946785588630915560e-04, /* 0x3F285F56, 0xB9CDF664 */
s03 =  1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */
s04 =  5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
s05 =  1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */

double j1(double x)
{
	double z,r,s;
	uint32_t ix;
	int sign;

	GET_HIGH_WORD(ix, x);
	sign = ix>>31;
	ix &= 0x7fffffff;
	if (ix >= 0x7ff00000)
		return 1/(x*x);
	if (ix >= 0x40000000)  /* |x| >= 2 */
		return common(ix, fabs(x), 0, sign);
	if (ix >= 0x38000000) {  /* |x| >= 2**-127 */
		z = x*x;
		r = z*(r00+z*(r01+z*(r02+z*r03)));
		s = 1+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
		z = r/s;
	} else
		/* avoid underflow, raise inexact if x!=0 */
		z = x;
	return (0.5 + z)*x;
}

static const double U0[5] = {
 -1.96057090646238940668e-01, /* 0xBFC91866, 0x143CBC8A */
  5.04438716639811282616e-02, /* 0x3FA9D3C7, 0x76292CD1 */
 -1.91256895875763547298e-03, /* 0xBF5F55E5, 0x4844F50F */
  2.35252600561610495928e-05, /* 0x3EF8AB03, 0x8FA6B88E */
 -9.19099158039878874504e-08, /* 0xBE78AC00, 0x569105B8 */
};
static const double V0[5] = {
  1.99167318236649903973e-02, /* 0x3F94650D, 0x3F4DA9F0 */
  2.02552581025135171496e-04, /* 0x3F2A8C89, 0x6C257764 */
  1.35608801097516229404e-06, /* 0x3EB6C05A, 0x894E8CA6 */
  6.22741452364621501295e-09, /* 0x3E3ABF1D, 0x5BA69A86 */
  1.66559246207992079114e-11, /* 0x3DB25039, 0xDACA772A */
};

double y1(double x)
{
	double z,u,v;
	uint32_t ix,lx;

	EXTRACT_WORDS(ix, lx, x);
	/* y1(nan)=nan, y1(<0)=nan, y1(0)=-inf, y1(inf)=0 */
	if ((ix<<1 | lx) == 0)
		return -1/0.0;
	if (ix>>31)
		return 0/0.0;
	if (ix >= 0x7ff00000)
		return 1/x;

	if (ix >= 0x40000000)  /* x >= 2 */
		return common(ix, x, 1, 0);
	if (ix < 0x3c900000)  /* x < 2**-54 */
		return -tpi/x;
	z = x*x;
	u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
	v = 1+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
	return x*(u/v) + tpi*(j1(x)*log(x)-1/x);
}

/* For x >= 8, the asymptotic expansions of pone is
 *      1 + 15/128 s^2 - 4725/2^15 s^4 - ...,   where s = 1/x.
 * We approximate pone by
 *      pone(x) = 1 + (R/S)
 * where  R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
 *        S = 1 + ps0*s^2 + ... + ps4*s^10
 * and
 *      | pone(x)-1-R/S | <= 2  ** ( -60.06)
 */

static const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
  1.17187499999988647970e-01, /* 0x3FBDFFFF, 0xFFFFFCCE */
  1.32394806593073575129e+01, /* 0x402A7A9D, 0x357F7FCE */
  4.12051854307378562225e+02, /* 0x4079C0D4, 0x652EA590 */
  3.87474538913960532227e+03, /* 0x40AE457D, 0xA3A532CC */
  7.91447954031891731574e+03, /* 0x40BEEA7A, 0xC32782DD */
};
static const double ps8[5] = {
  1.14207370375678408436e+02, /* 0x405C8D45, 0x8E656CAC */
  3.65093083420853463394e+03, /* 0x40AC85DC, 0x964D274F */
  3.69562060269033463555e+04, /* 0x40E20B86, 0x97C5BB7F */
  9.76027935934950801311e+04, /* 0x40F7D42C, 0xB28F17BB */
  3.08042720627888811578e+04, /* 0x40DE1511, 0x697A0B2D */
};

static const double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
  1.31990519556243522749e-11, /* 0x3DAD0667, 0xDAE1CA7D */
  1.17187493190614097638e-01, /* 0x3FBDFFFF, 0xE2C10043 */
  6.80275127868432871736e+00, /* 0x401B3604, 0x6E6315E3 */
  1.08308182990189109773e+02, /* 0x405B13B9, 0x452602ED */
  5.17636139533199752805e+02, /* 0x40802D16, 0xD052D649 */
  5.28715201363337541807e+02, /* 0x408085B8, 0xBB7E0CB7 */
};
static const double ps5[5] = {
  5.92805987221131331921e+01, /* 0x404DA3EA, 0xA8AF633D */
  9.91401418733614377743e+02, /* 0x408EFB36, 0x1B066701 */
  5.35326695291487976647e+03, /* 0x40B4E944, 0x5706B6FB */
  7.84469031749551231769e+03, /* 0x40BEA4B0, 0xB8A5BB15 */
  1.50404688810361062679e+03, /* 0x40978030, 0x036F5E51 */
};

static const double pr3[6] = {
  3.02503916137373618024e-09, /* 0x3E29FC21, 0xA7AD9EDD */
  1.17186865567253592491e-01, /* 0x3FBDFFF5, 0x5B21D17B */
  3.93297750033315640650e+00, /* 0x400F76BC, 0xE85EAD8A */
  3.51194035591636932736e+01, /* 0x40418F48, 0x9DA6D129 */
  9.10550110750781271918e+01, /* 0x4056C385, 0x4D2C1837 */
  4.85590685197364919645e+01, /* 0x4048478F, 0x8EA83EE5 */
};
static const double ps3[5] = {
  3.47913095001251519989e+01, /* 0x40416549, 0xA134069C */
  3.36762458747825746741e+02, /* 0x40750C33, 0x07F1A75F */
  1.04687139975775130551e+03, /* 0x40905B7C, 0x5037D523 */
  8.90811346398256432622e+02, /* 0x408BD67D, 0xA32E31E9 */
  1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */
};

static const double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
  1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */
  1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */
  2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */
  1.22426109148261232917e+01, /* 0x40287C37, 0x7F71A964 */
  1.76939711271687727390e+01, /* 0x4031B1A8, 0x177F8EE2 */
  5.07352312588818499250e+00, /* 0x40144B49, 0xA574C1FE */
};
static const double ps2[5] = {
  2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */
  1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */
  2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */
  1.17679373287147100768e+02, /* 0x405D6B7A, 0xDA1884A9 */
  8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
};

static double pone(double x)
{
	const double *p,*q;
	double_t z,r,s;
	uint32_t ix;

	GET_HIGH_WORD(ix, x);
	ix &= 0x7fffffff;
	if      (ix >= 0x40200000){p = pr8; q = ps8;}
	else if (ix >= 0x40122E8B){p = pr5; q = ps5;}
	else if (ix >= 0x4006DB6D){p = pr3; q = ps3;}
	else /*ix >= 0x40000000*/ {p = pr2; q = ps2;}
	z = 1.0/(x*x);
	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
	s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
	return 1.0+ r/s;
}

/* For x >= 8, the asymptotic expansions of qone is
 *      3/8 s - 105/1024 s^3 - ..., where s = 1/x.
 * We approximate pone by
 *      qone(x) = s*(0.375 + (R/S))
 * where  R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
 *        S = 1 + qs1*s^2 + ... + qs6*s^12
 * and
 *      | qone(x)/s -0.375-R/S | <= 2  ** ( -61.13)
 */

static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
 -1.02539062499992714161e-01, /* 0xBFBA3FFF, 0xFFFFFDF3 */
 -1.62717534544589987888e+01, /* 0xC0304591, 0xA26779F7 */
 -7.59601722513950107896e+02, /* 0xC087BCD0, 0x53E4B576 */
 -1.18498066702429587167e+04, /* 0xC0C724E7, 0x40F87415 */
 -4.84385124285750353010e+04, /* 0xC0E7A6D0, 0x65D09C6A */
};
static const double qs8[6] = {
  1.61395369700722909556e+02, /* 0x40642CA6, 0xDE5BCDE5 */
  7.82538599923348465381e+03, /* 0x40BE9162, 0xD0D88419 */
  1.33875336287249578163e+05, /* 0x4100579A, 0xB0B75E98 */
  7.19657723683240939863e+05, /* 0x4125F653, 0x72869C19 */
  6.66601232617776375264e+05, /* 0x412457D2, 0x7719AD5C */
 -2.94490264303834643215e+05, /* 0xC111F969, 0x0EA5AA18 */
};

static const double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
 -2.08979931141764104297e-11, /* 0xBDB6FA43, 0x1AA1A098 */
 -1.02539050241375426231e-01, /* 0xBFBA3FFF, 0xCB597FEF */
 -8.05644828123936029840e+00, /* 0xC0201CE6, 0xCA03AD4B */
 -1.83669607474888380239e+02, /* 0xC066F56D, 0x6CA7B9B0 */
 -1.37319376065508163265e+03, /* 0xC09574C6, 0x6931734F */
 -2.61244440453215656817e+03, /* 0xC0A468E3, 0x88FDA79D */
};
static const double qs5[6] = {
  8.12765501384335777857e+01, /* 0x405451B2, 0xFF5A11B2 */
  1.99179873460485964642e+03, /* 0x409F1F31, 0xE77BF839 */
  1.74684851924908907677e+04, /* 0x40D10F1F, 0x0D64CE29 */
  4.98514270910352279316e+04, /* 0x40E8576D, 0xAABAD197 */
  2.79480751638918118260e+04, /* 0x40DB4B04, 0xCF7C364B */
 -4.71918354795128470869e+03, /* 0xC0B26F2E, 0xFCFFA004 */
};

static const double qr3[6] = {
 -5.07831226461766561369e-09, /* 0xBE35CFA9, 0xD38FC84F */
 -1.02537829820837089745e-01, /* 0xBFBA3FEB, 0x51AEED54 */
 -4.61011581139473403113e+00, /* 0xC01270C2, 0x3302D9FF */
 -5.78472216562783643212e+01, /* 0xC04CEC71, 0xC25D16DA */
 -2.28244540737631695038e+02, /* 0xC06C87D3, 0x4718D55F */
 -2.19210128478909325622e+02, /* 0xC06B66B9, 0x5F5C1BF6 */
};
static const double qs3[6] = {
  4.76651550323729509273e+01, /* 0x4047D523, 0xCCD367E4 */
  6.73865112676699709482e+02, /* 0x40850EEB, 0xC031EE3E */
  3.38015286679526343505e+03, /* 0x40AA684E, 0x448E7C9A */
  5.54772909720722782367e+03, /* 0x40B5ABBA, 0xA61D54A6 */
  1.90311919338810798763e+03, /* 0x409DBC7A, 0x0DD4DF4B */
 -1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */
};

static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
 -1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */
 -1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */
 -2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */
 -1.96636162643703720221e+01, /* 0xC033A9E2, 0xC168907F */
 -4.23253133372830490089e+01, /* 0xC04529A3, 0xDE104AAA */
 -2.13719211703704061733e+01, /* 0xC0355F36, 0x39CF6E52 */
};
static const double qs2[6] = {
  2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */
  2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */
  7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */
  7.39393205320467245656e+02, /* 0x40871B25, 0x48D4C029 */
  1.55949003336666123687e+02, /* 0x40637E5E, 0x3C3ED8D4 */
 -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
};

static double qone(double x)
{
	const double *p,*q;
	double_t s,r,z;
	uint32_t ix;

	GET_HIGH_WORD(ix, x);
	ix &= 0x7fffffff;
	if      (ix >= 0x40200000){p = qr8; q = qs8;}
	else if (ix >= 0x40122E8B){p = qr5; q = qs5;}
	else if (ix >= 0x4006DB6D){p = qr3; q = qs3;}
	else /*ix >= 0x40000000*/ {p = qr2; q = qs2;}
	z = 1.0/(x*x);
	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
	s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
	return (.375 + r/s)/x;
}


text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_powf.c */
/*
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
 */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */

#include "libm.h"

static const float
bp[]   = {1.0, 1.5,},
dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
two24  =  16777216.0,  /* 0x4b800000 */
huge   =  1.0e30,
tiny   =  1.0e-30,
/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
L1 =  6.0000002384e-01, /* 0x3f19999a */
L2 =  4.2857143283e-01, /* 0x3edb6db7 */
L3 =  3.3333334327e-01, /* 0x3eaaaaab */
L4 =  2.7272811532e-01, /* 0x3e8ba305 */
L5 =  2.3066075146e-01, /* 0x3e6c3255 */
L6 =  2.0697501302e-01, /* 0x3e53f142 */
P1 =  1.6666667163e-01, /* 0x3e2aaaab */
P2 = -2.7777778450e-03, /* 0xbb360b61 */
P3 =  6.6137559770e-05, /* 0x388ab355 */
P4 = -1.6533901999e-06, /* 0xb5ddea0e */
P5 =  4.1381369442e-08, /* 0x3331bb4c */
lg2     =  6.9314718246e-01, /* 0x3f317218 */
lg2_h   =  6.93145752e-01,   /* 0x3f317200 */
lg2_l   =  1.42860654e-06,   /* 0x35bfbe8c */
ovt     =  4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
cp      =  9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
cp_h    =  9.6191406250e-01, /* 0x3f764000 =12b cp */
cp_l    = -1.1736857402e-04, /* 0xb8f623c6 =tail of cp_h */
ivln2   =  1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */
ivln2_h =  1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
ivln2_l =  7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/

float powf(float x, float y)
{
	float z,ax,z_h,z_l,p_h,p_l;
	float y1,t1,t2,r,s,sn,t,u,v,w;
	int32_t i,j,k,yisint,n;
	int32_t hx,hy,ix,iy,is;

	GET_FLOAT_WORD(hx, x);
	GET_FLOAT_WORD(hy, y);
	ix = hx & 0x7fffffff;
	iy = hy & 0x7fffffff;

	/* x**0 = 1, even if x is NaN */
	if (iy == 0)
		return 1.0f;
	/* 1**y = 1, even if y is NaN */
	if (hx == 0x3f800000)
		return 1.0f;
	/* NaN if either arg is NaN */
	if (ix > 0x7f800000 || iy > 0x7f800000)
		return x + y;

	/* determine if y is an odd int when x < 0
	 * yisint = 0       ... y is not an integer
	 * yisint = 1       ... y is an odd int
	 * yisint = 2       ... y is an even int
	 */
	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)      /* (-1)**+-inf is 1 */
			return 1.0f;
		else if (ix > 0x3f800000)  /* (|x|>1)**+-inf = inf,0 */
			return hy >= 0 ? y : 0.0f;
		else if (ix != 0)          /* (|x|<1)**+-inf = 0,inf if x!=0 */
			return hy >= 0 ? 0.0f: -y;
	}
	if (iy == 0x3f800000)    /* y is +-1 */
		return hy >= 0 ? x : 1.0f/x;
	if (hy == 0x40000000)    /* y is 2 */
		return x*x;
	if (hy == 0x3f000000) {  /* y is  0.5 */
		if (hx >= 0)     /* x >= +0 */
			return sqrtf(x);
	}

	ax = fabsf(x);
	/* special value of x */
	if (ix == 0x7f800000 || ix == 0 || ix == 0x3f800000) { /* x is +-0,+-inf,+-1 */
		z = ax;
		if (hy < 0)  /* z = (1/|x|) */
			z = 1.0f/z;
		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;
	}

	sn = 1.0f; /* sign of result */
	if (hx < 0) {
		if (yisint == 0) /* (x<0)**(non-int) is NaN */
			return (x-x)/(x-x);
		if (yisint == 1) /* (x<0)**(odd int) */
			sn = -1.0f;
	}

	/* |y| is huge */
	if (iy > 0x4d000000) { /* if |y| > 2**27 */
		/* over/underflow if x is not close to one */
		if (ix < 0x3f7ffff8)
			return hy < 0 ? sn*huge*huge : sn*tiny*tiny;
		if (ix > 0x3f800007)
			return hy > 0 ? sn*huge*huge : sn*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 - 1;     /* t has 20 trailing zeros */
		w = (t*t)*(0.5f - t*(0.333333333333f - t*0.25f));
		u = ivln2_h*t;  /* ivln2_h has 16 sig. bits */
		v = t*ivln2_l - w*ivln2;
		t1 = u + v;
		GET_FLOAT_WORD(is, t1);
		SET_FLOAT_WORD(t1, is & 0xfffff000);
		t2 = v - (t1-u);
	} else {
		float s2,s_h,s_l,t_h,t_l;
		n = 0;
		/* take care subnormal number */
		if (ix < 0x00800000) {
			ax *= two24;
			n -= 24;
			GET_FLOAT_WORD(ix, ax);
		}
		n += ((ix)>>23) - 0x7f;
		j = ix & 0x007fffff;
		/* determine interval */
		ix = j | 0x3f800000;     /* normalize ix */
		if (j <= 0x1cc471)       /* |x|<sqrt(3/2) */
			k = 0;
		else if (j < 0x5db3d7)   /* |x|<sqrt(3)   */
			k = 1;
		else {
			k = 0;
			n += 1;
			ix -= 0x00800000;
		}
		SET_FLOAT_WORD(ax, ix);

		/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
		u = ax - bp[k];   /* bp[0]=1.0, bp[1]=1.5 */
		v = 1.0f/(ax+bp[k]);
		s = u*v;
		s_h = s;
		GET_FLOAT_WORD(is, s_h);
		SET_FLOAT_WORD(s_h, is & 0xfffff000);
		/* t_h=ax+bp[k] High */
		is = ((ix>>1) & 0xfffff000) | 0x20000000;
		SET_FLOAT_WORD(t_h, is + 0x00400000 + (k<<21));
		t_l = ax - (t_h - bp[k]);
		s_l = v*((u - s_h*t_h) - s_h*t_l);
		/* compute log(ax) */
		s2 = s*s;
		r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
		r += s_l*(s_h+s);
		s2 = s_h*s_h;
		t_h = 3.0f + s2 + r;
		GET_FLOAT_WORD(is, t_h);
		SET_FLOAT_WORD(t_h, is & 0xfffff000);
		t_l = r - ((t_h - 3.0f) - s2);
		/* u+v = s*(1+...) */
		u = s_h*t_h;
		v = s_l*t_h + t_l*s;
		/* 2/(3log2)*(s+...) */
		p_h = u + v;
		GET_FLOAT_WORD(is, p_h);
		SET_FLOAT_WORD(p_h, is & 0xfffff000);
		p_l = v - (p_h - u);
		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) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
		t = (float)n;
		t1 = (((z_h + z_l) + dp_h[k]) + t);
		GET_FLOAT_WORD(is, t1);
		SET_FLOAT_WORD(t1, is & 0xfffff000);
		t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
	}

	/* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
	GET_FLOAT_WORD(is, y);
	SET_FLOAT_WORD(y1, is & 0xfffff000);
	p_l = (y-y1)*t1 + y*t2;
	p_h = y1*t1;
	z = p_l + p_h;
	GET_FLOAT_WORD(j, z);
	if (j > 0x43000000)          /* if z > 128 */
		return sn*huge*huge;  /* overflow */
	else if (j == 0x43000000) {  /* if z == 128 */
		if (p_l + ovt > z - p_h)
			return sn*huge*huge;  /* overflow */
	} else if ((j&0x7fffffff) > 0x43160000)  /* z < -150 */ // FIXME: check should be  (uint32_t)j > 0xc3160000
		return sn*tiny*tiny;  /* underflow */
	else if (j == 0xc3160000) {  /* z == -150 */
		if (p_l <= z-p_h)
			return sn*tiny*tiny;  /* 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 */
		SET_FLOAT_WORD(t, n & ~(0x007fffff>>k));
		n = ((n&0x007fffff)|0x00800000)>>(23-k);
		if (j < 0)
			n = -n;
		p_h -= t;
	}
	t = p_l + p_h;
	GET_FLOAT_WORD(is, t);
	SET_FLOAT_WORD(t, is & 0xffff8000);
	u = t*lg2_h;
	v = (p_l-(t-p_h))*lg2 + t*lg2_l;
	z = u + v;
	w = v - (z - u);
	t = z*z;
	t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
	r = (z*t1)/(t1-2.0f) - (w+z*w);
	z = 1.0f - (r - z);
	GET_FLOAT_WORD(j, z);
	j += n<<23;
	if ((j>>23) <= 0)  /* subnormal output */
		z = scalbnf(z, n);
	else
		SET_FLOAT_WORD(z, j);
	return sn*z;
}


text:
#include "libm.h"

/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
double atanh(double x)
{
	union {double f; uint64_t i;} u = {.f = x};
	unsigned e = u.i >> 52 & 0x7ff;
	unsigned s = u.i >> 63;
	double_t y;

	/* |x| */
	u.i &= (uint64_t)-1/2;
	y = u.f;

	if (e < 0x3ff - 1) {
		if (e < 0x3ff - 32) {
			/* handle underflow */
			if (e == 0)
				FORCE_EVAL((float)y);
		} else {
			/* |x| < 0.5, up to 1.7ulp error */
			y = 0.5*log1p(2*y + 2*y*y/(1-y));
		}
	} else {
		/* avoid overflow */
		y = 0.5*log1p(2*(y/(1-y)));
	}
	return s ? -y : y;
}


text:
/* origin: FreeBSD /usr/src/lib/msun/src/s_exp2.c */
/*-
 * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

#include "libm.h"

#define TBLSIZE 256

static const double
redux = 0x1.8p52 / TBLSIZE,
P1    = 0x1.62e42fefa39efp-1,
P2    = 0x1.ebfbdff82c575p-3,
P3    = 0x1.c6b08d704a0a6p-5,
P4    = 0x1.3b2ab88f70400p-7,
P5    = 0x1.5d88003875c74p-10;

static const double tbl[TBLSIZE * 2] = {
/*  exp2(z + eps)          eps     */
  0x1.6a09e667f3d5dp-1,  0x1.9880p-44,
  0x1.6b052fa751744p-1,  0x1.8000p-50,
  0x1.6c012750bd9fep-1, -0x1.8780p-45,
  0x1.6cfdcddd476bfp-1,  0x1.ec00p-46,
  0x1.6dfb23c651a29p-1, -0x1.8000p-50,
  0x1.6ef9298593ae3p-1, -0x1.c000p-52,
  0x1.6ff7df9519386p-1, -0x1.fd80p-45,
  0x1.70f7466f42da3p-1, -0x1.c880p-45,
  0x1.71f75e8ec5fc3p-1,  0x1.3c00p-46,
  0x1.72f8286eacf05p-1, -0x1.8300p-44,
  0x1.73f9a48a58152p-1, -0x1.0c00p-47,
  0x1.74fbd35d7ccfcp-1,  0x1.f880p-45,
  0x1.75feb564267f1p-1,  0x1.3e00p-47,
  0x1.77024b1ab6d48p-1, -0x1.7d00p-45,
  0x1.780694fde5d38p-1, -0x1.d000p-50,
  0x1.790b938ac1d00p-1,  0x1.3000p-49,
  0x1.7a11473eb0178p-1, -0x1.d000p-49,
  0x1.7b17b0976d060p-1,  0x1.0400p-45,
  0x1.7c1ed0130c133p-1,  0x1.0000p-53,
  0x1.7d26a62ff8636p-1, -0x1.6900p-45,
  0x1.7e2f336cf4e3bp-1, -0x1.2e00p-47,
  0x1.7f3878491c3e8p-1, -0x1.4580p-45,
  0x1.80427543e1b4ep-1,  0x1.3000p-44,
  0x1.814d2add1071ap-1,  0x1.f000p-47,
  0x1.82589994ccd7ep-1, -0x1.1c00p-45,
  0x1.8364c1eb942d0p-1,  0x1.9d00p-45,
  0x1.8471a4623cab5p-1,  0x1.7100p-43,
  0x1.857f4179f5bbcp-1,  0x1.2600p-45,
  0x1.868d99b4491afp-1, -0x1.2c40p-44,
  0x1.879cad931a395p-1, -0x1.3000p-45,
  0x1.88ac7d98a65b8p-1, -0x1.a800p-45,
  0x1.89bd0a4785800p-1, -0x1.d000p-49,
  0x1.8ace5422aa223p-1,  0x1.3280p-44,
  0x1.8be05bad619fap-1,  0x1.2b40p-43,
  0x1.8cf3216b54383p-1, -0x1.ed00p-45,
  0x1.8e06a5e08664cp-1, -0x1.0500p-45,
  0x1.8f1ae99157807p-1,  0x1.8280p-45,
  0x1.902fed0282c0ep-1, -0x1.cb00p-46,
  0x1.9145b0b91ff96p-1, -0x1.5e00p-47,
  0x1.925c353aa2ff9p-1,  0x1.5400p-48,
  0x1.93737b0cdc64ap-1,  0x1.7200p-46,
  0x1.948b82b5f98aep-1, -0x1.9000p-47,
  0x1.95a44cbc852cbp-1,  0x1.5680p-45,
  0x1.96bdd9a766f21p-1, -0x1.6d00p-44,
  0x1.97d829fde4e2ap-1, -0x1.1000p-47,
  0x1.98f33e47a23a3p-1,  0x1.d000p-45,
  0x1.9a0f170ca0604p-1, -0x1.8a40p-44,
  0x1.9b2bb4d53ff89p-1,  0x1.55c0p-44,
  0x1.9c49182a3f15bp-1,  0x1.6b80p-45,
  0x1.9d674194bb8c5p-1, -0x1.c000p-49,
  0x1.9e86319e3238ep-1,  0x1.7d00p-46,
  0x1.9fa5e8d07f302p-1,  0x1.6400p-46,
  0x1.a0c667b5de54dp-1, -0x1.5000p-48,
  0x1.a1e7aed8eb8f6p-1,  0x1.9e00p-47,
  0x1.a309bec4a2e27p-1,  0x1.ad80p-45,
  0x1.a42c980460a5dp-1, -0x1.af00p-46,
  0x1.a5503b23e259bp-1,  0x1.b600p-47,
  0x1.a674a8af46213p-1,  0x1.8880p-44,
  0x1.a799e1330b3a7p-1,  0x1.1200p-46,
  0x1.a8bfe53c12e8dp-1,  0x1.6c00p-47,
  0x1.a9e6b5579fcd2p-1, -0x1.9b80p-45,
  0x1.ab0e521356fb8p-1,  0x1.b700p-45,
  0x1.ac36bbfd3f381p-1,  0x1.9000p-50,
  0x1.ad5ff3a3c2780p-1,  0x1.4000p-49,
  0x1.ae89f995ad2a3p-1, -0x1.c900p-45,
  0x1.afb4ce622f367p-1,  0x1.6500p-46,
  0x1.b0e07298db790p-1,  0x1.fd40p-45,
  0x1.b20ce6c9a89a9p-1,  0x1.2700p-46,
  0x1.b33a2b84f1a4bp-1,  0x1.d470p-43,
  0x1.b468415b747e7p-1, -0x1.8380p-44,
  0x1.b59728de5593ap-1,  0x1.8000p-54,
  0x1.b6c6e29f1c56ap-1,  0x1.ad00p-47,
  0x1.b7f76f2fb5e50p-1,  0x1.e800p-50,
  0x1.b928cf22749b2p-1, -0x1.4c00p-47,
  0x1.ba5b030a10603p-1, -0x1.d700p-47,
  0x1.bb8e0b79a6f66p-1,  0x1.d900p-47,
  0x1.bcc1e904bc1ffp-1,  0x1.2a00p-47,
  0x1.bdf69c3f3a16fp-1, -0x1.f780p-46,
  0x1.bf2c25bd71db8p-1, -0x1.0a00p-46,
  0x1.c06286141b2e9p-1, -0x1.1400p-46,
  0x1.c199bdd8552e0p-1,  0x1.be00p-47,
  0x1.c2d1cd9fa64eep-1, -0x1.9400p-47,
  0x1.c40ab5fffd02fp-1, -0x1.ed00p-47,
  0x1.c544778fafd15p-1,  0x1.9660p-44,
  0x1.c67f12e57d0cbp-1, -0x1.a100p-46,
  0x1.c7ba88988c1b6p-1, -0x1.8458p-42,
  0x1.c8f6d9406e733p-1, -0x1.a480p-46,
  0x1.ca3405751c4dfp-1,  0x1.b000p-51,
  0x1.cb720dcef9094p-1,  0x1.1400p-47,
  0x1.ccb0f2e6d1689p-1,  0x1.0200p-48,
  0x1.cdf0b555dc412p-1,  0x1.3600p-48,
  0x1.cf3155b5bab3bp-1, -0x1.6900p-47,
  0x1.d072d4a0789bcp-1,  0x1.9a00p-47,
  0x1.d1b532b08c8fap-1, -0x1.5e00p-46,
  0x1.d2f87080d8a85p-1,  0x1.d280p-46,
  0x1.d43c8eacaa203p-1,  0x1.1a00p-47,
  0x1.d5818dcfba491p-1,  0x1.f000p-50,
  0x1.d6c76e862e6a1p-1, -0x1.3a00p-47,
  0x1.d80e316c9834ep-1, -0x1.cd80p-47,
  0x1.d955d71ff6090p-1,  0x1.4c00p-48,
  0x1.da9e603db32aep-1,  0x1.f900p-48,
  0x1.dbe7cd63a8325p-1,  0x1.9800p-49,
  0x1.dd321f301b445p-1, -0x1.5200p-48,
  0x1.de7d5641c05bfp-1, -0x1.d700p-46,
  0x1.dfc97337b9aecp-1, -0x1.6140p-46,
  0x1.e11676b197d5ep-1,  0x1.b480p-47,
  0x1.e264614f5a3e7p-1,  0x1.0ce0p-43,
  0x1.e3b333b16ee5cp-1,  0x1.c680p-47,
  0x1.e502ee78b3fb4p-1, -0x1.9300p-47,
  0x1.e653924676d68p-1, -0x1.5000p-49,
  0x1.e7a51fbc74c44p-1, -0x1.7f80p-47,
  0x1.e8f7977cdb726p-1, -0x1.3700p-48,
  0x1.ea4afa2a490e8p-1,  0x1.5d00p-49,
  0x1.eb9f4867ccae4p-1,  0x1.61a0p-46,
  0x1.ecf482d8e680dp-1,  0x1.5500p-48,
  0x1.ee4aaa2188514p-1,  0x1.6400p-51,
  0x1.efa1bee615a13p-1, -0x1.e800p-49,
  0x1.f0f9c1cb64106p-1, -0x1.a880p-48,
  0x1.f252b376bb963p-1, -0x1.c900p-45,
  0x1.f3ac948dd7275p-1,  0x1.a000p-53,
  0x1.f50765b6e4524p-1, -0x1.4f00p-48,
  0x1.f6632798844fdp-1,  0x1.a800p-51,
  0x1.f7bfdad9cbe38p-1,  0x1.abc0p-48,
  0x1.f91d802243c82p-1, -0x1.4600p-50,
  0x1.fa7c1819e908ep-1, -0x1.b0c0p-47,
  0x1.fbdba3692d511p-1, -0x1.0e00p-51,
  0x1.fd3c22b8f7194p-1, -0x1.0de8p-46,
  0x1.fe9d96b2a23eep-1,  0x1.e430p-49,
  0x1.0000000000000p+0,  0x0.0000p+0,
  0x1.00b1afa5abcbep+0, -0x1.3400p-52,
  0x1.0163da9fb3303p+0, -0x1.2170p-46,
  0x1.02168143b0282p+0,  0x1.a400p-52,
  0x1.02c9a3e77806cp+0,  0x1.f980p-49,
  0x1.037d42e11bbcap+0, -0x1.7400p-51,
  0x1.04315e86e7f89p+0,  0x1.8300p-50,
  0x1.04e5f72f65467p+0, -0x1.a3f0p-46,
  0x1.059b0d315855ap+0, -0x1.2840p-47,
  0x1.0650a0e3c1f95p+0,  0x1.1600p-48,
  0x1.0706b29ddf71ap+0,  0x1.5240p-46,
  0x1.07bd42b72a82dp+0, -0x1.9a00p-49,
  0x1.0874518759bd0p+0,  0x1.6400p-49,
  0x1.092bdf66607c8p+0, -0x1.0780p-47,
  0x1.09e3ecac6f383p+0, -0x1.8000p-54,
  0x1.0a9c79b1f3930p+0,  0x1.fa00p-48,
  0x1.0b5586cf988fcp+0, -0x1.ac80p-48,
  0x1.0c0f145e46c8ap+0,  0x1.9c00p-50,
  0x1.0cc922b724816p+0,  0x1.5200p-47,
  0x1.0d83b23395dd8p+0, -0x1.ad00p-48,
  0x1.0e3ec32d3d1f3p+0,  0x1.bac0p-46,
  0x1.0efa55fdfa9a6p+0, -0x1.4e80p-47,
  0x1.0fb66affed2f0p+0, -0x1.d300p-47,
  0x1.1073028d7234bp+0,  0x1.1500p-48,
  0x1.11301d0125b5bp+0,  0x1.c000p-49,
  0x1.11edbab5e2af9p+0,  0x1.6bc0p-46,
  0x1.12abdc06c31d5p+0,  0x1.8400p-49,
  0x1.136a814f2047dp+0, -0x1.ed00p-47,
  0x1.1429aaea92de9p+0,  0x1.8e00p-49,
  0x1.14e95934f3138p+0,  0x1.b400p-49,
  0x1.15a98c8a58e71p+0,  0x1.5300p-47,
  0x1.166a45471c3dfp+0,  0x1.3380p-47,
  0x1.172b83c7d5211p+0,  0x1.8d40p-45,
  0x1.17ed48695bb9fp+0, -0x1.5d00p-47,
  0x1.18af9388c8d93p+0, -0x1.c880p-46,
  0x1.1972658375d66p+0,  0x1.1f00p-46,
  0x1.1a35beb6fcba7p+0,  0x1.0480p-46,
  0x1.1af99f81387e3p+0, -0x1.7390p-43,
  0x1.1bbe084045d54p+0,  0x1.4e40p-45,
  0x1.1c82f95281c43p+0, -0x1.a200p-47,
  0x1.1d4873168b9b2p+0,  0x1.3800p-49,
  0x1.1e0e75eb44031p+0,  0x1.ac00p-49,
  0x1.1ed5022fcd938p+0,  0x1.1900p-47,
  0x1.1f9c18438cdf7p+0, -0x1.b780p-46,
  0x1.2063b88628d8fp+0,  0x1.d940p-45,
  0x1.212be3578a81ep+0,  0x1.8000p-50,
  0x1.21f49917ddd41p+0,  0x1.b340p-45,
  0x1.22bdda2791323p+0,  0x1.9f80p-46,
  0x1.2387a6e7561e7p+0, -0x1.9c80p-46,
  0x1.2451ffb821427p+0,  0x1.2300p-47,
  0x1.251ce4fb2a602p+0, -0x1.3480p-46,
  0x1.25e85711eceb0p+0,  0x1.2700p-46,
  0x1.26b4565e27d16p+0,  0x1.1d00p-46,
  0x1.2780e341de00fp+0,  0x1.1ee0p-44,
  0x1.284dfe1f5633ep+0, -0x1.4c00p-46,
  0x1.291ba7591bb30p+0, -0x1.3d80p-46,
  0x1.29e9df51fdf09p+0,  0x1.8b00p-47,
  0x1.2ab8a66d10e9bp+0, -0x1.27c0p-45,
  0x1.2b87fd0dada3ap+0,  0x1.a340p-45,
  0x1.2c57e39771af9p+0, -0x1.0800p-46,
  0x1.2d285a6e402d9p+0, -0x1.ed00p-47,
  0x1.2df961f641579p+0, -0x1.4200p-48,
  0x1.2ecafa93e2ecfp+0, -0x1.4980p-45,
  0x1.2f9d24abd8822p+0, -0x1.6300p-46,
  0x1.306fe0a31b625p+0, -0x1.2360p-44,
  0x1.31432edeea50bp+0, -0x1.0df8p-40,
  0x1.32170fc4cd7b8p+0, -0x1.2480p-45,
  0x1.32eb83ba8e9a2p+0, -0x1.5980p-45,
  0x1.33c08b2641766p+0,  0x1.ed00p-46,
  0x1.3496266e3fa27p+0, -0x1.c000p-50,
  0x1.356c55f929f0fp+0, -0x1.0d80p-44,
  0x1.36431a2de88b9p+0,  0x1.2c80p-45,
  0x1.371a7373aaa39p+0,  0x1.0600p-45,
  0x1.37f26231e74fep+0, -0x1.6600p-46,
  0x1.38cae6d05d838p+0, -0x1.ae00p-47,
  0x1.39a401b713ec3p+0, -0x1.4720p-43,
  0x1.3a7db34e5a020p+0,  0x1.8200p-47,
  0x1.3b57fbfec6e95p+0,  0x1.e800p-44,
  0x1.3c32dc313a8f2p+0,  0x1.f800p-49,
  0x1.3d0e544ede122p+0, -0x1.7a00p-46,
  0x1.3dea64c1234bbp+0,  0x1.6300p-45,
  0x1.3ec70df1c4eccp+0, -0x1.8a60p-43,
  0x1.3fa4504ac7e8cp+0, -0x1.cdc0p-44,
  0x1.40822c367a0bbp+0,  0x1.5b80p-45,
  0x1.4160a21f72e95p+0,  0x1.ec00p-46,
  0x1.423fb27094646p+0, -0x1.3600p-46,
  0x1.431f5d950a920p+0,  0x1.3980p-45,
  0x1.43ffa3f84b9ebp+0,  0x1.a000p-48,
  0x1.44e0860618919p+0, -0x1.6c00p-48,
  0x1.45c2042a7d201p+0, -0x1.bc00p-47,
  0x1.46a41ed1d0016p+0, -0x1.2800p-46,
  0x1.4786d668b3326p+0,  0x1.0e00p-44,
  0x1.486a2b5c13c00p+0, -0x1.d400p-45,
  0x1.494e1e192af04p+0,  0x1.c200p-47,
  0x1.4a32af0d7d372p+0, -0x1.e500p-46,
  0x1.4b17dea6db801p+0,  0x1.7800p-47,
  0x1.4bfdad53629e1p+0, -0x1.3800p-46,
  0x1.4ce41b817c132p+0,  0x1.0800p-47,
  0x1.4dcb299fddddbp+0,  0x1.c700p-45,
  0x1.4eb2d81d8ab96p+0, -0x1.ce00p-46,
  0x1.4f9b2769d2d02p+0,  0x1.9200p-46,
  0x1.508417f4531c1p+0, -0x1.8c00p-47,
  0x1.516daa2cf662ap+0, -0x1.a000p-48,
  0x1.5257de83f51eap+0,  0x1.a080p-43,
  0x1.5342b569d4edap+0, -0x1.6d80p-45,
  0x1.542e2f4f6ac1ap+0, -0x1.2440p-44,
  0x1.551a4ca5d94dbp+0,  0x1.83c0p-43,
  0x1.56070dde9116bp+0,  0x1.4b00p-45,
  0x1.56f4736b529dep+0,  0x1.15a0p-43,
  0x1.57e27dbe2c40ep+0, -0x1.9e00p-45,
  0x1.58d12d497c76fp+0, -0x1.3080p-45,
  0x1.59c0827ff0b4cp+0,  0x1.dec0p-43,
  0x1.5ab07dd485427p+0, -0x1.4000p-51,
  0x1.5ba11fba87af4p+0,  0x1.0080p-44,
  0x1.5c9268a59460bp+0, -0x1.6c80p-45,
  0x1.5d84590998e3fp+0,  0x1.69a0p-43,
  0x1.5e76f15ad20e1p+0, -0x1.b400p-46,
  0x1.5f6a320dcebcap+0,  0x1.7700p-46,
  0x1.605e1b976dcb8p+0,  0x1.6f80p-45,
  0x1.6152ae6cdf715p+0,  0x1.1000p-47,
  0x1.6247eb03a5531p+0, -0x1.5d00p-46,
  0x1.633dd1d1929b5p+0, -0x1.2d00p-46,
  0x1.6434634ccc313p+0, -0x1.a800p-49,
  0x1.652b9febc8efap+0, -0x1.8600p-45,
  0x1.6623882553397p+0,  0x1.1fe0p-40,
  0x1.671c1c708328ep+0, -0x1.7200p-44,
  0x1.68155d44ca97ep+0,  0x1.6800p-49,
  0x1.690f4b19e9471p+0, -0x1.9780p-45,
};

/*
 * exp2(x): compute the base 2 exponential of x
 *
 * Accuracy: Peak error < 0.503 ulp for normalized results.
 *
 * Method: (accurate tables)
 *
 *   Reduce x:
 *     x = k + y, for integer k and |y| <= 1/2.
 *     Thus we have exp2(x) = 2**k * exp2(y).
 *
 *   Reduce y:
 *     y = i/TBLSIZE + z - eps[i] for integer i near y * TBLSIZE.
 *     Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z - eps[i]),
 *     with |z - eps[i]| <= 2**-9 + 2**-39 for the table used.
 *
 *   We compute exp2(i/TBLSIZE) via table lookup and exp2(z - eps[i]) via
 *   a degree-5 minimax polynomial with maximum error under 1.3 * 2**-61.
 *   The values in exp2t[] and eps[] are chosen such that
 *   exp2t[i] = exp2(i/TBLSIZE + eps[i]), and eps[i] is a small offset such
 *   that exp2t[i] is accurate to 2**-64.
 *
 *   Note that the range of i is +-TBLSIZE/2, so we actually index the tables
 *   by i0 = i + TBLSIZE/2.  For cache efficiency, exp2t[] and eps[] are
 *   virtual tables, interleaved in the real table tbl[].
 *
 *   This method is due to Gal, with many details due to Gal and Bachelis:
 *
 *      Gal, S. and Bachelis, B.  An Accurate Elementary Mathematical Library
 *      for the IEEE Floating Point Standard.  TOMS 17(1), 26-46 (1991).
 */
double exp2(double x)
{
	double_t r, t, z;
	uint32_t ix, i0;
	union {double f; uint64_t i;} u = {x};
	union {uint32_t u; int32_t i;} k;

	/* Filter out exceptional cases. */
	ix = u.i>>32 & 0x7fffffff;
	if (ix >= 0x408ff000) {  /* |x| >= 1022 or nan */
		if (ix >= 0x40900000 && u.i>>63 == 0) {  /* x >= 1024 or nan */
			/* overflow */
			x *= 0x1p1023;
			return x;
		}
		if (ix >= 0x7ff00000)  /* -inf or -nan */
			return -1/x;
		if (u.i>>63) {  /* x <= -1022 */
			/* underflow */
			if (x <= -1075 || x - 0x1p52 + 0x1p52 != x)
				FORCE_EVAL((float)(-0x1p-149/x));
			if (x <= -1075)
				return 0;
		}
	} else if (ix < 0x3c900000) {  /* |x| < 0x1p-54 */
		return 1.0 + x;
	}

	/* Reduce x, computing z, i0, and k. */
	u.f = x + redux;
	i0 = u.i;
	i0 += TBLSIZE / 2;
	k.u = i0 / TBLSIZE * TBLSIZE;
	k.i /= TBLSIZE;
	i0 %= TBLSIZE;
	u.f -= redux;
	z = x - u.f;

	/* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
	t = tbl[2*i0];       /* exp2t[i0] */
	z -= tbl[2*i0 + 1];  /* eps[i0]   */
	r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5))));

	return scalbn(r, k.i);
}


text:
/* origin: OpenBSD /usr/src/lib/libm/src/polevll.c */
/*
 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
/*
 *      Evaluate polynomial
 *
 *
 * SYNOPSIS:
 *
 * int N;
 * long double x, y, coef[N+1], polevl[];
 *
 * y = polevll( x, coef, N );
 *
 *
 * DESCRIPTION:
 *
 * Evaluates polynomial of degree N:
 *
 *                     2          N
 * y  =  C  + C x + C x  +...+ C x
 *        0    1     2          N
 *
 * Coefficients are stored in reverse order:
 *
 * coef[0] = C  , ..., coef[N] = C  .
 *            N                   0
 *
 *  The function p1evll() assumes that coef[N] = 1.0 and is
 * omitted from the array.  Its calling arguments are
 * otherwise the same as polevll().
 *
 *
 * SPEED:
 *
 * In the interest of speed, there are no checks for out
 * of bounds arithmetic.  This routine is used by most of
 * the functions in the library.  Depending on available
 * equipment features, the user may wish to rewrite the
 * program in microcode or assembly language.
 *
 */

#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
#else
/*
 * Polynomial evaluator:
 *  P[0] x^n  +  P[1] x^(n-1)  +  ...  +  P[n]
 */
long double __polevll(long double x, const long double *P, int n)
{
	long double y;

	y = *P++;
	do {
		y = y * x + *P++;
	} while (--n);

	return y;
}

/*
 * Polynomial evaluator:
 *  x^n  +  P[0] x^(n-1)  +  P[1] x^(n-2)  +  ...  +  P[n]
 */
long double __p1evll(long double x, const long double *P, int n)
{
	long double y;

	n -= 1;
	y = x + *P++;
	do {
		y = y * x + *P++;
	} while (--n);

	return y;
}
#endif


text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_lgamma_r.c */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 *
 */
/* lgamma_r(x, signgamp)
 * Reentrant version of the logarithm of the Gamma function
 * with user provide pointer for the sign of Gamma(x).
 *
 * Method:
 *   1. Argument Reduction for 0 < x <= 8
 *      Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
 *      reduce x to a number in [1.5,2.5] by
 *              lgamma(1+s) = log(s) + lgamma(s)
 *      for example,
 *              lgamma(7.3) = log(6.3) + lgamma(6.3)
 *                          = log(6.3*5.3) + lgamma(5.3)
 *                          = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
 *   2. Polynomial approximation of lgamma around its
 *      minimun ymin=1.461632144968362245 to maintain monotonicity.
 *      On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
 *              Let z = x-ymin;
 *              lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
 *      where
 *              poly(z) is a 14 degree polynomial.
 *   2. Rational approximation in the primary interval [2,3]
 *      We use the following approximation:
 *              s = x-2.0;
 *              lgamma(x) = 0.5*s + s*P(s)/Q(s)
 *      with accuracy
 *              |P/Q - (lgamma(x)-0.5s)| < 2**-61.71
 *      Our algorithms are based on the following observation
 *
 *                             zeta(2)-1    2    zeta(3)-1    3
 * lgamma(2+s) = s*(1-Euler) + --------- * s  -  --------- * s  + ...
 *                                 2                 3
 *
 *      where Euler = 0.5771... is the Euler constant, which is very
 *      close to 0.5.
 *
 *   3. For x>=8, we have
 *      lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
 *      (better formula:
 *         lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
 *      Let z = 1/x, then we approximation
 *              f(z) = lgamma(x) - (x-0.5)(log(x)-1)
 *      by
 *                                  3       5             11
 *              w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
 *      where
 *              |w - f(z)| < 2**-58.74
 *
 *   4. For negative x, since (G is gamma function)
 *              -x*G(-x)*G(x) = pi/sin(pi*x),
 *      we have
 *              G(x) = pi/(sin(pi*x)*(-x)*G(-x))
 *      since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
 *      Hence, for x<0, signgam = sign(sin(pi*x)) and
 *              lgamma(x) = log(|Gamma(x)|)
 *                        = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
 *      Note: one should avoid compute pi*(-x) directly in the
 *            computation of sin(pi*(-x)).
 *
 *   5. Special Cases
 *              lgamma(2+s) ~ s*(1-Euler) for tiny s
 *              lgamma(1) = lgamma(2) = 0
 *              lgamma(x) ~ -log(|x|) for tiny x
 *              lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero
 *              lgamma(inf) = inf
 *              lgamma(-inf) = inf (bug for bug compatible with C99!?)
 *
 */

#include "libm.h"
#include "libc.h"

static const double
pi  =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
a0  =  7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
a1  =  3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
a2  =  6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */
a3  =  2.05808084325167332806e-02, /* 0x3F951322, 0xAC92547B */
a4  =  7.38555086081402883957e-03, /* 0x3F7E404F, 0xB68FEFE8 */
a5  =  2.89051383673415629091e-03, /* 0x3F67ADD8, 0xCCB7926B */
a6  =  1.19270763183362067845e-03, /* 0x3F538A94, 0x116F3F5D */
a7  =  5.10069792153511336608e-04, /* 0x3F40B6C6, 0x89B99C00 */
a8  =  2.20862790713908385557e-04, /* 0x3F2CF2EC, 0xED10E54D */
a9  =  1.08011567247583939954e-04, /* 0x3F1C5088, 0x987DFB07 */
a10 =  2.52144565451257326939e-05, /* 0x3EFA7074, 0x428CFA52 */
a11 =  4.48640949618915160150e-05, /* 0x3F07858E, 0x90A45837 */
tc  =  1.46163214496836224576e+00, /* 0x3FF762D8, 0x6356BE3F */
tf  = -1.21486290535849611461e-01, /* 0xBFBF19B9, 0xBCC38A42 */
/* tt = -(tail of tf) */
tt  = -3.63867699703950536541e-18, /* 0xBC50C7CA, 0xA48A971F */
t0  =  4.83836122723810047042e-01, /* 0x3FDEF72B, 0xC8EE38A2 */
t1  = -1.47587722994593911752e-01, /* 0xBFC2E427, 0x8DC6C509 */
t2  =  6.46249402391333854778e-02, /* 0x3FB08B42, 0x94D5419B */
t3  = -3.27885410759859649565e-02, /* 0xBFA0C9A8, 0xDF35B713 */
t4  =  1.79706750811820387126e-02, /* 0x3F9266E7, 0x970AF9EC */
t5  = -1.03142241298341437450e-02, /* 0xBF851F9F, 0xBA91EC6A */
t6  =  6.10053870246291332635e-03, /* 0x3F78FCE0, 0xE370E344 */
t7  = -3.68452016781138256760e-03, /* 0xBF6E2EFF, 0xB3E914D7 */
t8  =  2.25964780900612472250e-03, /* 0x3F6282D3, 0x2E15C915 */
t9  = -1.40346469989232843813e-03, /* 0xBF56FE8E, 0xBF2D1AF1 */
t10 =  8.81081882437654011382e-04, /* 0x3F4CDF0C, 0xEF61A8E9 */
t11 = -5.38595305356740546715e-04, /* 0xBF41A610, 0x9C73E0EC */
t12 =  3.15632070903625950361e-04, /* 0x3F34AF6D, 0x6C0EBBF7 */
t13 = -3.12754168375120860518e-04, /* 0xBF347F24, 0xECC38C38 */
t14 =  3.35529192635519073543e-04, /* 0x3F35FD3E, 0xE8C2D3F4 */
u0  = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
u1  =  6.32827064025093366517e-01, /* 0x3FE4401E, 0x8B005DFF */
u2  =  1.45492250137234768737e+00, /* 0x3FF7475C, 0xD119BD6F */
u3  =  9.77717527963372745603e-01, /* 0x3FEF4976, 0x44EA8450 */
u4  =  2.28963728064692451092e-01, /* 0x3FCD4EAE, 0xF6010924 */
u5  =  1.33810918536787660377e-02, /* 0x3F8B678B, 0xBF2BAB09 */
v1  =  2.45597793713041134822e+00, /* 0x4003A5D7, 0xC2BD619C */
v2  =  2.12848976379893395361e+00, /* 0x40010725, 0xA42B18F5 */
v3  =  7.69285150456672783825e-01, /* 0x3FE89DFB, 0xE45050AF */
v4  =  1.04222645593369134254e-01, /* 0x3FBAAE55, 0xD6537C88 */
v5  =  3.21709242282423911810e-03, /* 0x3F6A5ABB, 0x57D0CF61 */
s0  = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
s1  =  2.14982415960608852501e-01, /* 0x3FCB848B, 0x36E20878 */
s2  =  3.25778796408930981787e-01, /* 0x3FD4D98F, 0x4F139F59 */
s3  =  1.46350472652464452805e-01, /* 0x3FC2BB9C, 0xBEE5F2F7 */
s4  =  2.66422703033638609560e-02, /* 0x3F9B481C, 0x7E939961 */
s5  =  1.84028451407337715652e-03, /* 0x3F5E26B6, 0x7368F239 */
s6  =  3.19475326584100867617e-05, /* 0x3F00BFEC, 0xDD17E945 */
r1  =  1.39200533467621045958e+00, /* 0x3FF645A7, 0x62C4AB74 */
r2  =  7.21935547567138069525e-01, /* 0x3FE71A18, 0x93D3DCDC */
r3  =  1.71933865632803078993e-01, /* 0x3FC601ED, 0xCCFBDF27 */
r4  =  1.86459191715652901344e-02, /* 0x3F9317EA, 0x742ED475 */
r5  =  7.77942496381893596434e-04, /* 0x3F497DDA, 0xCA41A95B */
r6  =  7.32668430744625636189e-06, /* 0x3EDEBAF7, 0xA5B38140 */
w0  =  4.18938533204672725052e-01, /* 0x3FDACFE3, 0x90C97D69 */
w1  =  8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */
w2  = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */
w3  =  7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */
w4  = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
w5  =  8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
w6  = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */

/* sin(pi*x) assuming x > 2^-100, if sin(pi*x)==0 the sign is arbitrary */
static double sin_pi(double x)
{
	int n;

	/* spurious inexact if odd int */
	x = 2.0*(x*0.5 - floor(x*0.5));  /* x mod 2.0 */

	n = (int)(x*4.0);
	n = (n+1)/2;
	x -= n*0.5f;
	x *= pi;

	switch (n) {
	default: /* case 4: */
	case 0: return __sin(x, 0.0, 0);
	case 1: return __cos(x, 0.0);
	case 2: return __sin(-x, 0.0, 0);
	case 3: return -__cos(x, 0.0);
	}
}

double __lgamma_r(double x, int *signgamp)
{
	union {double f; uint64_t i;} u = {x};
	double_t t,y,z,nadj,p,p1,p2,p3,q,r,w;
	uint32_t ix;
	int sign,i;

	/* purge off +-inf, NaN, +-0, tiny and negative arguments */
	*signgamp = 1;
	sign = u.i>>63;
	ix = u.i>>32 & 0x7fffffff;
	if (ix >= 0x7ff00000)
		return x*x;
	if (ix < (0x3ff-70)<<20) {  /* |x|<2**-70, return -log(|x|) */
		if(sign) {
			x = -x;
			*signgamp = -1;
		}
		return -log(x);
	}
	if (sign) {
		x = -x;
		t = sin_pi(x);
		if (t == 0.0) /* -integer */
			return 1.0/(x-x);
		if (t > 0.0)
			*signgamp = -1;
		else
			t = -t;
		nadj = log(pi/(t*x));
	}

	/* purge off 1 and 2 */
	if ((ix == 0x3ff00000 || ix == 0x40000000) && (uint32_t)u.i == 0)
		r = 0;
	/* for x < 2.0 */
	else if (ix < 0x40000000) {
		if (ix <= 0x3feccccc) {   /* lgamma(x) = lgamma(x+1)-log(x) */
			r = -log(x);
			if (ix >= 0x3FE76944) {
				y = 1.0 - x;
				i = 0;
			} else if (ix >= 0x3FCDA661) {
				y = x - (tc-1.0);
				i = 1;
			} else {
				y = x;
				i = 2;
			}
		} else {
			r = 0.0;
			if (ix >= 0x3FFBB4C3) {  /* [1.7316,2] */
				y = 2.0 - x;
				i = 0;
			} else if(ix >= 0x3FF3B4C4) {  /* [1.23,1.73] */
				y = x - tc;
				i = 1;
			} else {
				y = x - 1.0;
				i = 2;
			}
		}
		switch (i) {
		case 0:
			z = y*y;
			p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
			p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
			p = y*p1+p2;
			r += (p-0.5*y);
			break;
		case 1:
			z = y*y;
			w = z*y;
			p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12)));    /* parallel comp */
			p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
			p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
			p = z*p1-(tt-w*(p2+y*p3));
			r += tf + p;
			break;
		case 2:
			p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
			p2 = 1.0+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
			r += -0.5*y + p1/p2;
		}
	} else if (ix < 0x40200000) {  /* x < 8.0 */
		i = (int)x;
		y = x - (double)i;
		p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
		q = 1.0+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
		r = 0.5*y+p/q;
		z = 1.0;    /* lgamma(1+s) = log(s) + lgamma(s) */
		switch (i) {
		case 7: z *= y + 6.0;  /* FALLTHRU */
		case 6: z *= y + 5.0;  /* FALLTHRU */
		case 5: z *= y + 4.0;  /* FALLTHRU */
		case 4: z *= y + 3.0;  /* FALLTHRU */
		case 3: z *= y + 2.0;  /* FALLTHRU */
			r += log(z);
			break;
		}
	} else if (ix < 0x43900000) {  /* 8.0 <= x < 2**58 */
		t = log(x);
		z = 1.0/x;
		y = z*z;
		w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
		r = (x-0.5)*(t-1.0)+w;
	} else                         /* 2**58 <= x <= inf */
		r =  x*(log(x)-1.0);
	if (sign)
		r = nadj - r;
	return r;
}

weak_alias(__lgamma_r, lgamma_r);


text:
#include <math.h>

long lroundf(float x)
{
	return roundf(x);
}


text:
#include <limits.h>
#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
int ilogbl(long double x)
{
	return ilogb(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
int ilogbl(long double x)
{
	#pragma STDC FENV_ACCESS ON
	union ldshape u = {x};
	uint64_t m = u.i.m;
	int e = u.i.se & 0x7fff;

	if (!e) {
		if (m == 0) {
			FORCE_EVAL(0/0.0f);
			return FP_ILOGB0;
		}
		/* subnormal x */
		for (e = -0x3fff+1; m>>63 == 0; e--, m<<=1);
		return e;
	}
	if (e == 0x7fff) {
		FORCE_EVAL(0/0.0f);
		return m<<1 ? FP_ILOGBNAN : INT_MAX;
	}
	return e - 0x3fff;
}
#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
int ilogbl(long double x)
{
	#pragma STDC FENV_ACCESS ON
	union ldshape u = {x};
	int e = u.i.se & 0x7fff;

	if (!e) {
		if (x == 0) {
			FORCE_EVAL(0/0.0f);
			return FP_ILOGB0;
		}
		/* subnormal x */
		x *= 0x1p120;
		return ilogbl(x) - 120;
	}
	if (e == 0x7fff) {
		FORCE_EVAL(0/0.0f);
		u.i.se = 0;
		return u.f ? FP_ILOGBNAN : INT_MAX;
	}
	return e - 0x3fff;
}
#endif


text:
#include <math.h>
#include <stdint.h>

float rintf(float x)
{
	union {float f; uint32_t i;} u = {x};
	int e = u.i>>23 & 0xff;
	int s = u.i>>31;
	float_t y;

	if (e >= 0x7f+23)
		return x;
	if (s)
		y = (float)(x - 0x1p23f) + 0x1p23f;
	else
		y = (float)(x + 0x1p23f) - 0x1p23f;
	if (y == 0)
		return s ? -0.0f : 0.0f;
	return y;
}


text:
/* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtf.c */
/*
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
 * Debugged and optimized by Bruce D. Evans.
 */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */
/* cbrtf(x)
 * Return cube root of x
 */

#include <math.h>
#include <stdint.h>

static const unsigned
B1 = 709958130, /* B1 = (127-127.0/3-0.03306235651)*2**23 */
B2 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */

float cbrtf(float x)
{
	double_t r,T;
	union {float f; uint32_t i;} u = {x};
	uint32_t hx = u.i & 0x7fffffff;

	if (hx >= 0x7f800000)  /* cbrt(NaN,INF) is itself */
		return x + x;

	/* rough cbrt to 5 bits */
	if (hx < 0x00800000) {  /* zero or subnormal? */
		if (hx == 0)
			return x;  /* cbrt(+-0) is itself */
		u.f = x*0x1p24f;
		hx = u.i & 0x7fffffff;
		hx = hx/3 + B2;
	} else
		hx = hx/3 + B1;
	u.i &= 0x80000000;
	u.i |= hx;

	/*
	 * First step Newton iteration (solving t*t-x/t == 0) to 16 bits.  In
	 * double precision so that its terms can be arranged for efficiency
	 * without causing overflow or underflow.
	 */
	T = u.f;
	r = T*T*T;
	T = T*((double_t)x+x+r)/(x+r+r);

	/*
	 * Second step Newton iteration to 47 bits.  In double precision for
	 * efficiency and accuracy.
	 */
	r = T*T*T;
	T = T*((double_t)x+x+r)/(x+r+r);

	/* rounding to 24 bits is perfect in round-to-nearest mode */
	return T;
}


text:
#define _GNU_SOURCE
#include <math.h>

float significandf(float x)
{
	return scalbnf(x, -ilogbf(x));
}


text:
#include <math.h>
#include <stdint.h>

double fmod(double x, double y)
{
	union {double f; uint64_t i;} ux = {x}, uy = {y};
	int ex = ux.i>>52 & 0x7ff;
	int ey = uy.i>>52 & 0x7ff;
	int sx = ux.i>>63;
	uint64_t i;

	/* in the followings uxi should be ux.i, but then gcc wrongly adds */
	/* float load/store to inner loops ruining performance and code size */
	uint64_t uxi = ux.i;

	if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff)
		return (x*y)/(x*y);
	if (uxi<<1 <= uy.i<<1) {
		if (uxi<<1 == uy.i<<1)
			return 0*x;
		return x;
	}

	/* normalize x and y */
	if (!ex) {
		for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1);
		uxi <<= -ex + 1;
	} else {
		uxi &= -1ULL >> 12;
		uxi |= 1ULL << 52;
	}
	if (!ey) {
		for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1);
		uy.i <<= -ey + 1;
	} else {
		uy.i &= -1ULL >> 12;
		uy.i |= 1ULL << 52;
	}

	/* x mod y */
	for (; ex > ey; ex--) {
		i = uxi - uy.i;
		if (i >> 63 == 0) {
			if (i == 0)
				return 0*x;
			uxi = i;
		}
		uxi <<= 1;
	}
	i = uxi - uy.i;
	if (i >> 63 == 0) {
		if (i == 0)
			return 0*x;
		uxi = i;
	}
	for (; uxi>>52 == 0; uxi <<= 1, ex--);

	/* scale result */
	if (ex > 0) {
		uxi -= 1ULL << 52;
		uxi |= (uint64_t)ex << 52;
	} else {
		uxi >>= -ex + 1;
	}
	uxi |= (uint64_t)sx << 63;
	ux.i = uxi;
	return ux.f;
}


text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_exp.c */
/*
 * ====================================================
 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
 *
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */
/* exp(x)
 * Returns the exponential of x.
 *
 * Method
 *   1. Argument reduction:
 *      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
 *      Given x, find r and integer k such that
 *
 *               x = k*ln2 + r,  |r| <= 0.5*ln2.
 *
 *      Here r will be represented as r = hi-lo for better
 *      accuracy.
 *
 *   2. Approximation of exp(r) by a special rational function on
 *      the interval [0,0.34658]:
 *      Write
 *          R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
 *      We use a special Remez algorithm on [0,0.34658] to generate
 *      a polynomial of degree 5 to approximate R. The maximum error
 *      of this polynomial approximation is bounded by 2**-59. In
 *      other words,
 *          R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
 *      (where z=r*r, and the values of P1 to P5 are listed below)
 *      and
 *          |                  5          |     -59
 *          | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
 *          |                             |
 *      The computation of exp(r) thus becomes
 *                              2*r
 *              exp(r) = 1 + ----------
 *                            R(r) - r
 *                                 r*c(r)
 *                     = 1 + r + ----------- (for better accuracy)
 *                                2 - c(r)
 *      where
 *                              2       4             10
 *              c(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
 *
 *   3. Scale back to obtain exp(x):
 *      From step 1, we have
 *         exp(x) = 2^k * exp(r)
 *
 * Special cases:
 *      exp(INF) is INF, exp(NaN) is NaN;
 *      exp(-INF) is 0, and
 *      for finite argument, only exp(0)=1 is exact.
 *
 * Accuracy:
 *      according to an error analysis, the error is always less than
 *      1 ulp (unit in the last place).
 *
 * Misc. info.
 *      For IEEE double
 *          if x >  709.782712893383973096 then exp(x) overflows
 *          if x < -745.133219101941108420 then exp(x) underflows
 */

#include "libm.h"

static const double
half[2] = {0.5,-0.5},
ln2hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
ln2lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */

double exp(double x)
{
	double_t hi, lo, c, xx, y;
	int k, sign;
	uint32_t hx;

	GET_HIGH_WORD(hx, x);
	sign = hx>>31;
	hx &= 0x7fffffff;  /* high word of |x| */

	/* special cases */
	if (hx >= 0x4086232b) {  /* if |x| >= 708.39... */
		if (isnan(x))
			return x;
		if (x > 709.782712893383973096) {
			/* overflow if x!=inf */
			x *= 0x1p1023;
			return x;
		}
		if (x < -708.39641853226410622) {
			/* underflow if x!=-inf */
			FORCE_EVAL((float)(-0x1p-149/x));
			if (x < -745.13321910194110842)
				return 0;
		}
	}

	/* argument reduction */
	if (hx > 0x3fd62e42) {  /* if |x| > 0.5 ln2 */
		if (hx >= 0x3ff0a2b2)  /* if |x| >= 1.5 ln2 */
			k = (int)(invln2*x + half[sign]);
		else
			k = 1 - sign - sign;
		hi = x - k*ln2hi;  /* k*ln2hi is exact here */
		lo = k*ln2lo;
		x = hi - lo;
	} else if (hx > 0x3e300000)  {  /* if |x| > 2**-28 */
		k = 0;
		hi = x;
		lo = 0;
	} else {
		/* inexact if x!=0 */
		FORCE_EVAL(0x1p1023 + x);
		return 1 + x;
	}

	/* x is now in primary range */
	xx = x*x;
	c = x - xx*(P1+xx*(P2+xx*(P3+xx*(P4+xx*P5))));
	y = 1 + (x*c/(2-c) - lo + hi);
	if (k == 0)
		return y;
	return scalbn(y, k);
}


text:
#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double asinhl(long double x)
{
	return asinh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */
long double asinhl(long double x)
{
	union ldshape u = {x};
	unsigned e = u.i.se & 0x7fff;
	unsigned s = u.i.se >> 15;

	/* |x| */
	u.i.se = e;
	x = u.f;

	if (e >= 0x3fff + 32) {
		/* |x| >= 0x1p32 or inf or nan */
		x = logl(x) + 0.693147180559945309417232121458176568L;
	} else if (e >= 0x3fff + 1) {
		/* |x| >= 2 */
		x = logl(2*x + 1/(sqrtl(x*x+1)+x));
	} else if (e >= 0x3fff - 32) {
		/* |x| >= 0x1p-32 */
		x = log1pl(x + x*x/(sqrtl(x*x+1)+1));
	} else {
		/* |x| < 0x1p-32, raise inexact if x!=0 */
		FORCE_EVAL(x + 0x1p120f);
	}
	return s ? -x : x;
}
#endif


text:
#include <math.h>

float fmaxf(float x, float y)
{
	if (isnan(x))
		return y;
	if (isnan(y))
		return x;
	/* handle signed zeroes, see C99 Annex F.9.9.2 */
	if (signbit(x) != signbit(y))
		return signbit(x) ? y : x;
	return x < y ? y : x;
}


text:
/* origin: FreeBSD /usr/src/lib/msun/ld80/k_tanl.c */
/* origin: FreeBSD /usr/src/lib/msun/ld128/k_tanl.c */
/*
 * ====================================================
 * Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
 * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
 *
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */

#include "libm.h"

#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#if LDBL_MANT_DIG == 64
/*
 * ld80 version of __tan.c.  See __tan.c for most comments.
 */
/*
 * Domain [-0.67434, 0.67434], range ~[-2.25e-22, 1.921e-22]
 * |tan(x)/x - t(x)| < 2**-71.9
 *
 * See __cosl.c for more details about the polynomial.
 */
static const long double
T3 =  0.333333333333333333180L,         /*  0xaaaaaaaaaaaaaaa5.0p-65 */
T5 =  0.133333333333333372290L,         /*  0x88888888888893c3.0p-66 */
T7 =  0.0539682539682504975744L,        /*  0xdd0dd0dd0dc13ba2.0p-68 */
pio4   =  0.785398163397448309628L,     /*  0xc90fdaa22168c235.0p-64 */
pio4lo = -1.25413940316708300586e-20L;  /* -0xece675d1fc8f8cbb.0p-130 */
static const double
T9  =  0.021869488536312216,            /*  0x1664f4882cc1c2.0p-58 */
T11 =  0.0088632355256619590,           /*  0x1226e355c17612.0p-59 */
T13 =  0.0035921281113786528,           /*  0x1d6d3d185d7ff8.0p-61 */
T15 =  0.0014558334756312418,           /*  0x17da354aa3f96b.0p-62 */
T17 =  0.00059003538700862256,          /*  0x13559358685b83.0p-63 */
T19 =  0.00023907843576635544,          /*  0x1f56242026b5be.0p-65 */
T21 =  0.000097154625656538905,         /*  0x1977efc26806f4.0p-66 */
T23 =  0.000038440165747303162,         /*  0x14275a09b3ceac.0p-67 */
T25 =  0.000018082171885432524,         /*  0x12f5e563e5487e.0p-68 */
T27 =  0.0000024196006108814377,        /*  0x144c0d80cc6896.0p-71 */
T29 =  0.0000078293456938132840,        /*  0x106b59141a6cb3.0p-69 */
T31 = -0.0000032609076735050182,        /* -0x1b5abef3ba4b59.0p-71 */
T33 =  0.0000023261313142559411;        /*  0x13835436c0c87f.0p-71 */
#define RPOLY(w) (T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + \
	w * (T25 + w * (T29 + w * T33)))))))
#define VPOLY(w) (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + \
	w * (T27 + w * T31))))))
#elif LDBL_MANT_DIG == 113
/*
 * ld128 version of __tan.c.  See __tan.c for most comments.
 */
/*
 * Domain [-0.67434, 0.67434], range ~[-3.37e-36, 1.982e-37]
 * |tan(x)/x - t(x)| < 2**-117.8 (XXX should be ~1e-37)
 *
 * See __cosl.c for more details about the polynomial.
 */
static const long double
T3 = 0x1.5555555555555555555555555553p-2L,
T5 = 0x1.1111111111111111111111111eb5p-3L,
T7 = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5L,
T9 = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6L,
T11 = 0x1.226e355e6c23c8f5b4f5762322eep-7L,
T13 = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9L,
T15 = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10L,
T17 = 0x1.355824803674477dfcf726649efep-11L,
T19 = 0x1.f57d7734d1656e0aceb716f614c2p-13L,
T21 = 0x1.967e18afcb180ed942dfdc518d6cp-14L,
T23 = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15L,
T25 = 0x1.0b132d39f055c81be49eff7afd50p-16L,
T27 = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18L,
T29 = 0x1.5ef2daf21d1113df38d0fbc00267p-19L,
T31 = 0x1.1c77d6eac0234988cdaa04c96626p-20L,
T33 = 0x1.cd2a5a292b180e0bdd701057dfe3p-22L,
T35 = 0x1.75c7357d0298c01a31d0a6f7d518p-23L,
T37 = 0x1.2f3190f4718a9a520f98f50081fcp-24L,
pio4 = 0x1.921fb54442d18469898cc51701b8p-1L,
pio4lo = 0x1.cd129024e088a67cc74020bbea60p-116L;
static const double
T39 =  0.000000028443389121318352,	/*  0x1e8a7592977938.0p-78 */
T41 =  0.000000011981013102001973,	/*  0x19baa1b1223219.0p-79 */
T43 =  0.0000000038303578044958070,	/*  0x107385dfb24529.0p-80 */
T45 =  0.0000000034664378216909893,	/*  0x1dc6c702a05262.0p-81 */
T47 = -0.0000000015090641701997785,	/* -0x19ecef3569ebb6.0p-82 */
T49 =  0.0000000029449552300483952,	/*  0x194c0668da786a.0p-81 */
T51 = -0.0000000022006995706097711,	/* -0x12e763b8845268.0p-81 */
T53 =  0.0000000015468200913196612,	/*  0x1a92fc98c29554.0p-82 */
T55 = -0.00000000061311613386849674,	/* -0x151106cbc779a9.0p-83 */
T57 =  1.4912469681508012e-10;		/*  0x147edbdba6f43a.0p-85 */
#define RPOLY(w) (T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + \
	w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 + \
	w * (T45 + w * (T49 + w * (T53 + w * T57)))))))))))))
#define VPOLY(w) (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + \
	w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 + \
	w * (T47 + w * (T51 + w * T55))))))))))))
#endif

long double __tanl(long double x, long double y, int odd) {
	long double z, r, v, w, s, a, t;
	int big, sign;

	big = fabsl(x) >= 0.67434;
	if (big) {
		sign = 0;
		if (x < 0) {
			sign = 1;
			x = -x;
			y = -y;
		}
		x = (pio4 - x) + (pio4lo - y);
		y = 0.0;
	}
	z = x * x;
	w = z * z;
	r = RPOLY(w);
	v = z * VPOLY(w);
	s = z * x;
	r = y + z * (s * (r + v) + y) + T3 * s;
	w = x + r;
	if (big) {
		s = 1 - 2*odd;
		v = s - 2.0 * (x + (r - w * w / (w + s)));
		return sign ? -v : v;
	}
	if (!odd)
		return w;
	/*
	 * if allow error up to 2 ulp, simply return
	 * -1.0 / (x+r) here
	 */
	/* compute -1.0 / (x+r) accurately */
	z = w;
	z = z + 0x1p32 - 0x1p32;
	v = r - (z - x);        /* z+v = r+x */
	t = a = -1.0 / w;       /* a = -1.0/w */
	t = t + 0x1p32 - 0x1p32;
	s = 1.0 + t * z;
	return t + a * (s + t * v);
}
#endif


text:
#include <math.h>

double ldexp(double x, int n)
{
	return scalbn(x, n);
}


text:
/* origin: FreeBSD /usr/src/lib/msun/ld80/s_exp2l.c */
/*-
 * Copyright (c) 2005-2008 David Schultz <das@FreeBSD.ORG>
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double exp2l(long double x)
{
	return exp2(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
#define TBLBITS 7
#define TBLSIZE (1 << TBLBITS)

static const double
redux = 0x1.8p63 / TBLSIZE,
P1    = 0x1.62e42fefa39efp-1,
P2    = 0x1.ebfbdff82c58fp-3,
P3    = 0x1.c6b08d7049fap-5,
P4    = 0x1.3b2ab6fba4da5p-7,
P5    = 0x1.5d8804780a736p-10,
P6    = 0x1.430918835e33dp-13;

static const double tbl[TBLSIZE * 2] = {
	0x1.6a09e667f3bcdp-1,   -0x1.bdd3413b2648p-55,
	0x1.6c012750bdabfp-1,   -0x1.2895667ff0cp-57,
	0x1.6dfb23c651a2fp-1,   -0x1.bbe3a683c88p-58,
	0x1.6ff7df9519484p-1,   -0x1.83c0f25860fp-56,
	0x1.71f75e8ec5f74p-1,   -0x1.16e4786887bp-56,
	0x1.73f9a48a58174p-1,   -0x1.0a8d96c65d5p-55,
	0x1.75feb564267c9p-1,   -0x1.0245957316ep-55,
	0x1.780694fde5d3fp-1,    0x1.866b80a0216p-55,
	0x1.7a11473eb0187p-1,   -0x1.41577ee0499p-56,
	0x1.7c1ed0130c132p-1,    0x1.f124cd1164ep-55,
	0x1.7e2f336cf4e62p-1,    0x1.05d02ba157ap-57,
	0x1.80427543e1a12p-1,   -0x1.27c86626d97p-55,
	0x1.82589994cce13p-1,   -0x1.d4c1dd41533p-55,
	0x1.8471a4623c7adp-1,   -0x1.8d684a341cep-56,
	0x1.868d99b4492edp-1,   -0x1.fc6f89bd4f68p-55,
	0x1.88ac7d98a6699p-1,    0x1.994c2f37cb5p-55,
	0x1.8ace5422aa0dbp-1,    0x1.6e9f156864bp-55,
	0x1.8cf3216b5448cp-1,   -0x1.0d55e32e9e4p-57,
	0x1.8f1ae99157736p-1,    0x1.5cc13a2e397p-56,
	0x1.9145b0b91ffc6p-1,   -0x1.dd6792e5825p-55,
	0x1.93737b0cdc5e5p-1,   -0x1.75fc781b58p-58,
	0x1.95a44cbc8520fp-1,   -0x1.64b7c96a5fp-57,
	0x1.97d829fde4e5p-1,    -0x1.d185b7c1b86p-55,
	0x1.9a0f170ca07bap-1,   -0x1.173bd91cee6p-55,
	0x1.9c49182a3f09p-1,     0x1.c7c46b071f2p-57,
	0x1.9e86319e32323p-1,    0x1.824ca78e64cp-57,
	0x1.a0c667b5de565p-1,   -0x1.359495d1cd5p-55,
	0x1.a309bec4a2d33p-1,    0x1.6305c7ddc368p-55,
	0x1.a5503b23e255dp-1,   -0x1.d2f6edb8d42p-55,
	0x1.a799e1330b358p-1,    0x1.bcb7ecac564p-55,
	0x1.a9e6b5579fdbfp-1,    0x1.0fac90ef7fdp-55,
	0x1.ac36bbfd3f37ap-1,   -0x1.f9234cae76dp-56,
	0x1.ae89f995ad3adp-1,    0x1.7a1cd345dcc8p-55,
	0x1.b0e07298db666p-1,   -0x1.bdef54c80e4p-55,
	0x1.b33a2b84f15fbp-1,   -0x1.2805e3084d8p-58,
	0x1.b59728de5593ap-1,   -0x1.c71dfbbba6ep-55,
	0x1.b7f76f2fb5e47p-1,   -0x1.5584f7e54acp-57,
	0x1.ba5b030a1064ap-1,   -0x1.efcd30e5429p-55,
	0x1.bcc1e904bc1d2p-1,    0x1.23dd07a2d9fp-56,
	0x1.bf2c25bd71e09p-1,   -0x1.efdca3f6b9c8p-55,
	0x1.c199bdd85529cp-1,    0x1.11065895049p-56,
	0x1.c40ab5fffd07ap-1,    0x1.b4537e083c6p-55,
	0x1.c67f12e57d14bp-1,    0x1.2884dff483c8p-55,
	0x1.c8f6d9406e7b5p-1,    0x1.1acbc48805cp-57,
	0x1.cb720dcef9069p-1,    0x1.503cbd1e94ap-57,
	0x1.cdf0b555dc3fap-1,   -0x1.dd83b53829dp-56,
	0x1.d072d4a07897cp-1,   -0x1.cbc3743797a8p-55,
	0x1.d2f87080d89f2p-1,   -0x1.d487b719d858p-55,
	0x1.d5818dcfba487p-1,    0x1.2ed02d75b37p-56,
	0x1.d80e316c98398p-1,   -0x1.11ec18bedep-55,
	0x1.da9e603db3285p-1,    0x1.c2300696db5p-55,
	0x1.dd321f301b46p-1,     0x1.2da5778f019p-55,
	0x1.dfc97337b9b5fp-1,   -0x1.1a5cd4f184b8p-55,
	0x1.e264614f5a129p-1,   -0x1.7b627817a148p-55,
	0x1.e502ee78b3ff6p-1,    0x1.39e8980a9cdp-56,
	0x1.e7a51fbc74c83p-1,    0x1.2d522ca0c8ep-55,
	0x1.ea4afa2a490dap-1,   -0x1.e9c23179c288p-55,
	0x1.ecf482d8e67f1p-1,   -0x1.c93f3b411ad8p-55,
	0x1.efa1bee615a27p-1,    0x1.dc7f486a4b68p-55,
	0x1.f252b376bba97p-1,    0x1.3a1a5bf0d8e8p-55,
	0x1.f50765b6e454p-1,     0x1.9d3e12dd8a18p-55,
	0x1.f7bfdad9cbe14p-1,   -0x1.dbb12d00635p-55,
	0x1.fa7c1819e90d8p-1,    0x1.74853f3a593p-56,
	0x1.fd3c22b8f71f1p-1,    0x1.2eb74966578p-58,
	0x1p+0,                  0x0p+0,
	0x1.0163da9fb3335p+0,    0x1.b61299ab8cd8p-54,
	0x1.02c9a3e778061p+0,   -0x1.19083535b08p-56,
	0x1.04315e86e7f85p+0,   -0x1.0a31c1977c98p-54,
	0x1.059b0d3158574p+0,    0x1.d73e2a475b4p-55,
	0x1.0706b29ddf6dep+0,   -0x1.c91dfe2b13cp-55,
	0x1.0874518759bc8p+0,    0x1.186be4bb284p-57,
	0x1.09e3ecac6f383p+0,    0x1.14878183161p-54,
	0x1.0b5586cf9890fp+0,    0x1.8a62e4adc61p-54,
	0x1.0cc922b7247f7p+0,    0x1.01edc16e24f8p-54,
	0x1.0e3ec32d3d1a2p+0,    0x1.03a1727c58p-59,
	0x1.0fb66affed31bp+0,   -0x1.b9bedc44ebcp-57,
	0x1.11301d0125b51p+0,   -0x1.6c51039449bp-54,
	0x1.12abdc06c31ccp+0,   -0x1.1b514b36ca8p-58,
	0x1.1429aaea92dep+0,    -0x1.32fbf9af1368p-54,
	0x1.15a98c8a58e51p+0,    0x1.2406ab9eeabp-55,
	0x1.172b83c7d517bp+0,   -0x1.19041b9d78ap-55,
	0x1.18af9388c8deap+0,   -0x1.11023d1970f8p-54,
	0x1.1a35beb6fcb75p+0,    0x1.e5b4c7b4969p-55,
	0x1.1bbe084045cd4p+0,   -0x1.95386352ef6p-54,
	0x1.1d4873168b9aap+0,    0x1.e016e00a264p-54,
	0x1.1ed5022fcd91dp+0,   -0x1.1df98027bb78p-54,
	0x1.2063b88628cd6p+0,    0x1.dc775814a85p-55,
	0x1.21f49917ddc96p+0,    0x1.2a97e9494a6p-55,
	0x1.2387a6e756238p+0,    0x1.9b07eb6c7058p-54,
	0x1.251ce4fb2a63fp+0,    0x1.ac155bef4f5p-55,
	0x1.26b4565e27cddp+0,    0x1.2bd339940eap-55,
	0x1.284dfe1f56381p+0,   -0x1.a4c3a8c3f0d8p-54,
	0x1.29e9df51fdee1p+0,    0x1.612e8afad12p-55,
	0x1.2b87fd0dad99p+0,    -0x1.10adcd6382p-59,
	0x1.2d285a6e4030bp+0,    0x1.0024754db42p-54,
	0x1.2ecafa93e2f56p+0,    0x1.1ca0f45d524p-56,
	0x1.306fe0a31b715p+0,    0x1.6f46ad23183p-55,
	0x1.32170fc4cd831p+0,    0x1.a9ce78e1804p-55,
	0x1.33c08b26416ffp+0,    0x1.327218436598p-54,
	0x1.356c55f929ff1p+0,   -0x1.b5cee5c4e46p-55,
	0x1.371a7373aa9cbp+0,   -0x1.63aeabf42ebp-54,
	0x1.38cae6d05d866p+0,   -0x1.e958d3c99048p-54,
	0x1.3a7db34e59ff7p+0,   -0x1.5e436d661f6p-56,
	0x1.3c32dc313a8e5p+0,   -0x1.efff8375d2ap-54,
	0x1.3dea64c123422p+0,    0x1.ada0911f09fp-55,
	0x1.3fa4504ac801cp+0,   -0x1.7d023f956fap-54,
	0x1.4160a21f72e2ap+0,   -0x1.ef3691c309p-58,
	0x1.431f5d950a897p+0,   -0x1.1c7dde35f7ap-55,
	0x1.44e086061892dp+0,    0x1.89b7a04ef8p-59,
	0x1.46a41ed1d0057p+0,    0x1.c944bd1648a8p-54,
	0x1.486a2b5c13cdp+0,     0x1.3c1a3b69062p-56,
	0x1.4a32af0d7d3dep+0,    0x1.9cb62f3d1be8p-54,
	0x1.4bfdad5362a27p+0,    0x1.d4397afec42p-56,
	0x1.4dcb299fddd0dp+0,    0x1.8ecdbbc6a78p-54,
	0x1.4f9b2769d2ca7p+0,   -0x1.4b309d25958p-54,
	0x1.516daa2cf6642p+0,   -0x1.f768569bd94p-55,
	0x1.5342b569d4f82p+0,   -0x1.07abe1db13dp-55,
	0x1.551a4ca5d920fp+0,   -0x1.d689cefede6p-55,
	0x1.56f4736b527dap+0,    0x1.9bb2c011d938p-54,
	0x1.58d12d497c7fdp+0,    0x1.295e15b9a1ep-55,
	0x1.5ab07dd485429p+0,    0x1.6324c0546478p-54,
	0x1.5c9268a5946b7p+0,    0x1.c4b1b81698p-60,
	0x1.5e76f15ad2148p+0,    0x1.ba6f93080e68p-54,
	0x1.605e1b976dc09p+0,   -0x1.3e2429b56de8p-54,
	0x1.6247eb03a5585p+0,   -0x1.383c17e40b48p-54,
	0x1.6434634ccc32p+0,    -0x1.c483c759d89p-55,
	0x1.6623882552225p+0,   -0x1.bb60987591cp-54,
	0x1.68155d44ca973p+0,    0x1.038ae44f74p-57,
};

/*
 * exp2l(x): compute the base 2 exponential of x
 *
 * Accuracy: Peak error < 0.511 ulp.
 *
 * Method: (equally-spaced tables)
 *
 *   Reduce x:
 *     x = 2**k + y, for integer k and |y| <= 1/2.
 *     Thus we have exp2l(x) = 2**k * exp2(y).
 *
 *   Reduce y:
 *     y = i/TBLSIZE + z for integer i near y * TBLSIZE.
 *     Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z),
 *     with |z| <= 2**-(TBLBITS+1).
 *
 *   We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a
 *   degree-6 minimax polynomial with maximum error under 2**-69.
 *   The table entries each have 104 bits of accuracy, encoded as
 *   a pair of double precision values.
 */
long double exp2l(long double x)
{
	union ldshape u = {x};
	int e = u.i.se & 0x7fff;
	long double r, z;
	uint32_t i0;
	union {uint32_t u; int32_t i;} k;

	/* Filter out exceptional cases. */
	if (e >= 0x3fff + 13) {  /* |x| >= 8192 or x is NaN */
		if (u.i.se >= 0x3fff + 14 && u.i.se >> 15 == 0)
			/* overflow */
			return x * 0x1p16383L;
		if (e == 0x7fff)  /* -inf or -nan */
			return -1/x;
		if (x < -16382) {
			if (x <= -16446 || x - 0x1p63 + 0x1p63 != x)
				/* underflow */
				FORCE_EVAL((float)(-0x1p-149/x));
			if (x <= -16446)
				return 0;
		}
	} else if (e < 0x3fff - 64) {
		return 1 + x;
	}

	/*
	 * Reduce x, computing z, i0, and k. The low bits of x + redux
	 * contain the 16-bit integer part of the exponent (k) followed by
	 * TBLBITS fractional bits (i0). We use bit tricks to extract these
	 * as integers, then set z to the remainder.
	 *
	 * Example: Suppose x is 0xabc.123456p0 and TBLBITS is 8.
	 * Then the low-order word of x + redux is 0x000abc12,
	 * We split this into k = 0xabc and i0 = 0x12 (adjusted to
	 * index into the table), then we compute z = 0x0.003456p0.
	 */
	u.f = x + redux;
	i0 = u.i.m + TBLSIZE / 2;
	k.u = i0 / TBLSIZE * TBLSIZE;
	k.i /= TBLSIZE;
	i0 %= TBLSIZE;
	u.f -= redux;
	z = x - u.f;

	/* Compute r = exp2l(y) = exp2lt[i0] * p(z). */
	long double t_hi = tbl[2*i0];
	long double t_lo = tbl[2*i0 + 1];
	/* XXX This gives > 1 ulp errors outside of FE_TONEAREST mode */
	r = t_lo + (t_hi + t_lo) * z * (P1 + z * (P2 + z * (P3 + z * (P4
	     + z * (P5 + z * P6))))) + t_hi;

	return scalbnl(r, k.i);
}
#endif


text:
#include <math.h>
#include <float.h>

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double fdiml(long double x, long double y)
{
	return fdim(x, y);
}
#else
long double fdiml(long double x, long double y)
{
	if (isnan(x))
		return x;
	if (isnan(y))
		return y;
	return x > y ? x - y : 0;
}
#endif


text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_expf.c */
/*
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
 */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */

#include "libm.h"

static const float
half[2] = {0.5,-0.5},
ln2hi   = 6.9314575195e-1f,  /* 0x3f317200 */
ln2lo   = 1.4286067653e-6f,  /* 0x35bfbe8e */
invln2  = 1.4426950216e+0f,  /* 0x3fb8aa3b */
/*
 * Domain [-0.34568, 0.34568], range ~[-4.278e-9, 4.447e-9]:
 * |x*(exp(x)+1)/(exp(x)-1) - p(x)| < 2**-27.74
 */
P1 =  1.6666625440e-1f, /*  0xaaaa8f.0p-26 */
P2 = -2.7667332906e-3f; /* -0xb55215.0p-32 */

float expf(float x)
{
	float_t hi, lo, c, xx, y;
	int k, sign;
	uint32_t hx;

	GET_FLOAT_WORD(hx, x);
	sign = hx >> 31;   /* sign bit of x */
	hx &= 0x7fffffff;  /* high word of |x| */

	/* special cases */
	if (hx >= 0x42aeac50) {  /* if |x| >= -87.33655f or NaN */
		if (hx >= 0x42b17218 && !sign) {  /* x >= 88.722839f */
			/* overflow */
			x *= 0x1p127f;
			return x;
		}
		if (sign) {
			/* underflow */
			FORCE_EVAL(-0x1p-149f/x);
			if (hx >= 0x42cff1b5)  /* x <= -103.972084f */
				return 0;
		}
	}

	/* argument reduction */
	if (hx > 0x3eb17218) {  /* if |x| > 0.5 ln2 */
		if (hx > 0x3f851592)  /* if |x| > 1.5 ln2 */
			k = invln2*x + half[sign];
		else
			k = 1 - sign - sign;
		hi = x - k*ln2hi;  /* k*ln2hi is exact here */
		lo = k*ln2lo;
		x = hi - lo;
	} else if (hx > 0x39000000) {  /* |x| > 2**-14 */
		k = 0;
		hi = x;
		lo = 0;
	} else {
		/* raise inexact */
		FORCE_EVAL(0x1p127f + x);
		return 1 + x;
	}

	/* x is now in primary range */
	xx = x*x;
	c = x - xx*(P1+xx*P2);
	y = 1 + (x*c/(2-c) - lo + hi);
	if (k == 0)
		return y;
	return scalbnf(y, k);
}


text:
/* origin: FreeBSD /usr/src/lib/msun/src/k_cos.c */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */
/*
 * __cos( x,  y )
 * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
 * Input x is assumed to be bounded by ~pi/4 in magnitude.
 * Input y is the tail of x.
 *
 * Algorithm
 *      1. Since cos(-x) = cos(x), we need only to consider positive x.
 *      2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
 *      3. cos(x) is approximated by a polynomial of degree 14 on
 *         [0,pi/4]
 *                                       4            14
 *              cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
 *         where the remez error is
 *
 *      |              2     4     6     8     10    12     14 |     -58
 *      |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
 *      |                                                      |
 *
 *                     4     6     8     10    12     14
 *      4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
 *             cos(x) ~ 1 - x*x/2 + r
 *         since cos(x+y) ~ cos(x) - sin(x)*y
 *                        ~ cos(x) - x*y,
 *         a correction term is necessary in cos(x) and hence
 *              cos(x+y) = 1 - (x*x/2 - (r - x*y))
 *         For better accuracy, rearrange to
 *              cos(x+y) ~ w + (tmp + (r-x*y))
 *         where w = 1 - x*x/2 and tmp is a tiny correction term
 *         (1 - x*x/2 == w + tmp exactly in infinite precision).
 *         The exactness of w + tmp in infinite precision depends on w
 *         and tmp having the same precision as x.  If they have extra
 *         precision due to compiler bugs, then the extra precision is
 *         only good provided it is retained in all terms of the final
 *         expression for cos().  Retention happens in all cases tested
 *         under FreeBSD, so don't pessimize things by forcibly clipping
 *         any extra precision in w.
 */

#include "libm.h"

static const double
C1  =  4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
C2  = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
C3  =  2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
C4  = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
C5  =  2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */

double __cos(double x, double y)
{
	double_t hz,z,r,w;

	z  = x*x;
	w  = z*z;
	r  = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
	hz = 0.5*z;
	w  = 1.0-hz;
	return w + (((1.0-w)-hz) + (z*r-x*y));
}


text:
#include "libm.h"

// FIXME: macro in math.h
int __signbit(double x)
{
	union {
		double d;
		uint64_t i;
	} y = { x };
	return y.i>>63;
}


text:
#include <math.h>

double nan(const char *s)
{
	return NAN;
}


text:
#include <math.h>

long lroundl(long double x)
{
	return roundl(x);
}


text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_asinf.c */
/*
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
 */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */
#include "libm.h"

static const double
pio2 = 1.570796326794896558e+00;

static const float
/* coefficients for R(x^2) */
pS0 =  1.6666586697e-01,
pS1 = -4.2743422091e-02,
pS2 = -8.6563630030e-03,
qS1 = -7.0662963390e-01;

static float R(float z)
{
	float_t p, q;
	p = z*(pS0+z*(pS1+z*pS2));
	q = 1.0f+z*qS1;
	return p/q;
}

float asinf(float x)
{
	double s;
	float z;
	uint32_t hx,ix;

	GET_FLOAT_WORD(hx, x);
	ix = hx & 0x7fffffff;
	if (ix >= 0x3f800000) {  /* |x| >= 1 */
		if (ix == 0x3f800000)  /* |x| == 1 */
			return x*pio2 + 0x1p-120f;  /* asin(+-1) = +-pi/2 with inexact */
		return 0/(x-x);  /* asin(|x|>1) is NaN */
	}
	if (ix < 0x3f000000) {  /* |x| < 0.5 */
		/* if 0x1p-126 <= |x| < 0x1p-12, avoid raising underflow */
		if (ix < 0x39800000 && ix >= 0x00800000)
			return x;
		return x + x*R(x*x);
	}
	/* 1 > |x| >= 0.5 */
	z = (1 - fabsf(x))*0.5f;
	s = sqrt(z);
	x = pio2 - 2*(s+s*R(z));
	if (hx >> 31)
		return -x;
	return x;
}


text:
#include <math.h>
#include <float.h>

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double nearbyintl(long double x)
{
	return nearbyint(x);
}
#else
#include <fenv.h>
long double nearbyintl(long double x)
{
#ifdef FE_INEXACT
	#pragma STDC FENV_ACCESS ON
	int e;

	e = fetestexcept(FE_INEXACT);
#endif
	x = rintl(x);
#ifdef FE_INEXACT
	if (!e)
		feclearexcept(FE_INEXACT);
#endif
	return x;
}
#endif