text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_jn.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.
 * ====================================================
 */
/*
 * jn(n, x), yn(n, x)
 * floating point Bessel's function of the 1st and 2nd kind
 * of order n
 *
 * Special cases:
 *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
 *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
 * Note 2. About jn(n,x), yn(n,x)
 *      For n=0, j0(x) is called,
 *      for n=1, j1(x) is called,
 *      for n<=x, forward recursion is used starting
 *      from values of j0(x) and j1(x).
 *      for n>x, a continued fraction approximation to
 *      j(n,x)/j(n-1,x) is evaluated and then backward
 *      recursion is used starting from a supposed value
 *      for j(n,x). The resulting value of j(0,x) is
 *      compared with the actual value to correct the
 *      supposed value of j(n,x).
 *
 *      yn(n,x) is similar in all respects, except
 *      that forward recursion is used for all
 *      values of n>1.
 */

#include "libm.h"

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

double jn(int n, double x)
{
	uint32_t ix, lx;
	int nm1, i, sign;
	double a, b, temp;

	EXTRACT_WORDS(ix, lx, x);
	sign = ix>>31;
	ix &= 0x7fffffff;

	if ((ix | (lx|-lx)>>31) > 0x7ff00000) /* nan */
		return x;

	/* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
	 * Thus, J(-n,x) = J(n,-x)
	 */
	/* nm1 = |n|-1 is used instead of |n| to handle n==INT_MIN */
	if (n == 0)
		return j0(x);
	if (n < 0) {
		nm1 = -(n+1);
		x = -x;
		sign ^= 1;
	} else
		nm1 = n-1;
	if (nm1 == 0)
		return j1(x);

	sign &= n;  /* even n: 0, odd n: signbit(x) */
	x = fabs(x);
	if ((ix|lx) == 0 || ix == 0x7ff00000)  /* if x is 0 or inf */
		b = 0.0;
	else if (nm1 < x) {
		/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
		if (ix >= 0x52d00000) { /* x > 2**302 */
			/* (x >> n**2)
			 *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
			 *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
			 *      Let s=sin(x), c=cos(x),
			 *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
			 *
			 *             n    sin(xn)*sqt2    cos(xn)*sqt2
			 *          ----------------------------------
			 *             0     s-c             c+s
			 *             1    -s-c            -c+s
			 *             2    -s+c            -c-s
			 *             3     s+c             c-s
			 */
			switch(nm1&3) {
			case 0: temp = -cos(x)+sin(x); break;
			case 1: temp = -cos(x)-sin(x); break;
			case 2: temp =  cos(x)-sin(x); break;
			default:
			case 3: temp =  cos(x)+sin(x); break;
			}
			b = invsqrtpi*temp/sqrt(x);
		} else {
			a = j0(x);
			b = j1(x);
			for (i=0; i<nm1; ) {
				i++;
				temp = b;
				b = b*(2.0*i/x) - a; /* avoid underflow */
				a = temp;
			}
		}
	} else {
		if (ix < 0x3e100000) { /* x < 2**-29 */
			/* x is tiny, return the first Taylor expansion of J(n,x)
			 * J(n,x) = 1/n!*(x/2)^n  - ...
			 */
			if (nm1 > 32)  /* underflow */
				b = 0.0;
			else {
				temp = x*0.5;
				b = temp;
				a = 1.0;
				for (i=2; i<=nm1+1; i++) {
					a *= (double)i; /* a = n! */
					b *= temp;      /* b = (x/2)^n */
				}
				b = b/a;
			}
		} else {
			/* use backward recurrence */
			/*                      x      x^2      x^2
			 *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
			 *                      2n  - 2(n+1) - 2(n+2)
			 *
			 *                      1      1        1
			 *  (for large x)   =  ----  ------   ------   .....
			 *                      2n   2(n+1)   2(n+2)
			 *                      -- - ------ - ------ -
			 *                       x     x         x
			 *
			 * Let w = 2n/x and h=2/x, then the above quotient
			 * is equal to the continued fraction:
			 *                  1
			 *      = -----------------------
			 *                     1
			 *         w - -----------------
			 *                        1
			 *              w+h - ---------
			 *                     w+2h - ...
			 *
			 * To determine how many terms needed, let
			 * Q(0) = w, Q(1) = w(w+h) - 1,
			 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
			 * When Q(k) > 1e4      good for single
			 * When Q(k) > 1e9      good for double
			 * When Q(k) > 1e17     good for quadruple
			 */
			/* determine k */
			double t,q0,q1,w,h,z,tmp,nf;
			int k;

			nf = nm1 + 1.0;
			w = 2*nf/x;
			h = 2/x;
			z = w+h;
			q0 = w;
			q1 = w*z - 1.0;
			k = 1;
			while (q1 < 1.0e9) {
				k += 1;
				z += h;
				tmp = z*q1 - q0;
				q0 = q1;
				q1 = tmp;
			}
			for (t=0.0, i=k; i>=0; i--)
				t = 1/(2*(i+nf)/x - t);
			a = t;
			b = 1.0;
			/*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
			 *  Hence, if n*(log(2n/x)) > ...
			 *  single 8.8722839355e+01
			 *  double 7.09782712893383973096e+02
			 *  long double 1.1356523406294143949491931077970765006170e+04
			 *  then recurrent value may overflow and the result is
			 *  likely underflow to zero
			 */
			tmp = nf*log(fabs(w));
			if (tmp < 7.09782712893383973096e+02) {
				for (i=nm1; i>0; i--) {
					temp = b;
					b = b*(2.0*i)/x - a;
					a = temp;
				}
			} else {
				for (i=nm1; i>0; i--) {
					temp = b;
					b = b*(2.0*i)/x - a;
					a = temp;
					/* scale b to avoid spurious overflow */
					if (b > 0x1p500) {
						a /= b;
						t /= b;
						b  = 1.0;
					}
				}
			}
			z = j0(x);
			w = j1(x);
			if (fabs(z) >= fabs(w))
				b = t*z/b;
			else
				b = t*w/a;
		}
	}
	return sign ? -b : b;
}

double yn(int n, double x)
{
	uint32_t ix, lx, ib;
	int nm1, sign, i;
	double a, b, temp;

	EXTRACT_WORDS(ix, lx, x);
	sign = ix>>31;
	ix &= 0x7fffffff;

	if ((ix | (lx|-lx)>>31) > 0x7ff00000) /* nan */
		return x;
	if (sign && (ix|lx)!=0) /* x < 0 */
		return 0/0.0;
	if (ix == 0x7ff00000)
		return 0.0;

	if (n == 0)
		return y0(x);
	if (n < 0) {
		nm1 = -(n+1);
		sign = n&1;
	} else {
		nm1 = n-1;
		sign = 0;
	}
	if (nm1 == 0)
		return sign ? -y1(x) : y1(x);

	if (ix >= 0x52d00000) { /* x > 2**302 */
		/* (x >> n**2)
		 *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
		 *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
		 *      Let s=sin(x), c=cos(x),
		 *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
		 *
		 *             n    sin(xn)*sqt2    cos(xn)*sqt2
		 *          ----------------------------------
		 *             0     s-c             c+s
		 *             1    -s-c            -c+s
		 *             2    -s+c            -c-s
		 *             3     s+c             c-s
		 */
		switch(nm1&3) {
		case 0: temp = -sin(x)-cos(x); break;
		case 1: temp = -sin(x)+cos(x); break;
		case 2: temp =  sin(x)+cos(x); break;
		default:
		case 3: temp =  sin(x)-cos(x); break;
		}
		b = invsqrtpi*temp/sqrt(x);
	} else {
		a = y0(x);
		b = y1(x);
		/* quit if b is -inf */
		GET_HIGH_WORD(ib, b);
		for (i=0; i<nm1 && ib!=0xfff00000; ){
			i++;
			temp = b;
			b = (2.0*i/x)*b - a;
			GET_HIGH_WORD(ib, b);
			a = temp;
		}
	}
	return sign ? -b : b;
}


text:
#include <math.h>

/* uses LONG_MAX > 2^24, see comments in lrint.c */

long lrintf(float x)
{
	return rintf(x);
}

repo_language: C

mime_type: text/x-c


text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_j0f.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.
 * ====================================================
 */

#define _GNU_SOURCE
#include "libm.h"

static float pzerof(float), qzerof(float);

static const float
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi       = 6.3661974669e-01; /* 0x3f22f983 */

static float common(uint32_t ix, float x, int y0)
{
	float z,s,c,ss,cc;
	/*
	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
	 */
	s = sinf(x);
	c = cosf(x);
	if (y0)
		c = -c;
	cc = s+c;
	if (ix < 0x7f000000) {
		ss = s-c;
		z = -cosf(2*x);
		if (s*c < 0)
			cc = z/ss;
		else
			ss = z/cc;
		if (ix < 0x58800000) {
			if (y0)
				ss = -ss;
			cc = pzerof(x)*cc-qzerof(x)*ss;
		}
	}
	return invsqrtpi*cc/sqrtf(x);
}

/* R0/S0 on [0, 2.00] */
static const float
R02 =  1.5625000000e-02, /* 0x3c800000 */
R03 = -1.8997929874e-04, /* 0xb947352e */
R04 =  1.8295404516e-06, /* 0x35f58e88 */
R05 = -4.6183270541e-09, /* 0xb19eaf3c */
S01 =  1.5619102865e-02, /* 0x3c7fe744 */
S02 =  1.1692678527e-04, /* 0x38f53697 */
S03 =  5.1354652442e-07, /* 0x3509daa6 */
S04 =  1.1661400734e-09; /* 0x30a045e8 */

float j0f(float x)
{
	float z,r,s;
	uint32_t ix;

	GET_FLOAT_WORD(ix, x);
	ix &= 0x7fffffff;
	if (ix >= 0x7f800000)
		return 1/(x*x);
	x = fabsf(x);

	if (ix >= 0x40000000) {  /* |x| >= 2 */
		/* large ulp error near zeros */
		return common(ix, x, 0);
	}
	if (ix >= 0x3a000000) {  /* |x| >= 2**-11 */
		/* up to 4ulp error near 2 */
		z = x*x;
		r = z*(R02+z*(R03+z*(R04+z*R05)));
		s = 1+z*(S01+z*(S02+z*(S03+z*S04)));
		return (1+x/2)*(1-x/2) + z*(r/s);
	}
	if (ix >= 0x21800000)  /* |x| >= 2**-60 */
		x = 0.25f*x*x;
	return 1 - x;
}

static const float
u00  = -7.3804296553e-02, /* 0xbd9726b5 */
u01  =  1.7666645348e-01, /* 0x3e34e80d */
u02  = -1.3818567619e-02, /* 0xbc626746 */
u03  =  3.4745343146e-04, /* 0x39b62a69 */
u04  = -3.8140706238e-06, /* 0xb67ff53c */
u05  =  1.9559013964e-08, /* 0x32a802ba */
u06  = -3.9820518410e-11, /* 0xae2f21eb */
v01  =  1.2730483897e-02, /* 0x3c509385 */
v02  =  7.6006865129e-05, /* 0x389f65e0 */
v03  =  2.5915085189e-07, /* 0x348b216c */
v04  =  4.4111031494e-10; /* 0x2ff280c2 */

float y0f(float x)
{
	float z,u,v;
	uint32_t ix;

	GET_FLOAT_WORD(ix, x);
	if ((ix & 0x7fffffff) == 0)
		return -1/0.0f;
	if (ix>>31)
		return 0/0.0f;
	if (ix >= 0x7f800000)
		return 1/x;
	if (ix >= 0x40000000) {  /* |x| >= 2.0 */
		/* large ulp error near zeros */
		return common(ix,x,1);
	}
	if (ix >= 0x39000000) {  /* x >= 2**-13 */
		/* large ulp error at x ~= 0.89 */
		z = x*x;
		u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
		v = 1+z*(v01+z*(v02+z*(v03+z*v04)));
		return u/v + tpi*(j0f(x)*logf(x));
	}
	return u00 + tpi*logf(x);
}

/* The asymptotic expansions of pzero is
 *      1 - 9/128 s^2 + 11025/98304 s^4 - ...,  where s = 1/x.
 * For x >= 2, We approximate pzero by
 *      pzero(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
 *      | pzero(x)-1-R/S | <= 2  ** ( -60.26)
 */
static const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
  0.0000000000e+00, /* 0x00000000 */
 -7.0312500000e-02, /* 0xbd900000 */
 -8.0816707611e+00, /* 0xc1014e86 */
 -2.5706311035e+02, /* 0xc3808814 */
 -2.4852163086e+03, /* 0xc51b5376 */
 -5.2530439453e+03, /* 0xc5a4285a */
};
static const float pS8[5] = {
  1.1653436279e+02, /* 0x42e91198 */
  3.8337448730e+03, /* 0x456f9beb */
  4.0597855469e+04, /* 0x471e95db */
  1.1675296875e+05, /* 0x47e4087c */
  4.7627726562e+04, /* 0x473a0bba */
};
static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
 -1.1412546255e-11, /* 0xad48c58a */
 -7.0312492549e-02, /* 0xbd8fffff */
 -4.1596107483e+00, /* 0xc0851b88 */
 -6.7674766541e+01, /* 0xc287597b */
 -3.3123129272e+02, /* 0xc3a59d9b */
 -3.4643338013e+02, /* 0xc3ad3779 */
};
static const float pS5[5] = {
  6.0753936768e+01, /* 0x42730408 */
  1.0512523193e+03, /* 0x44836813 */
  5.9789707031e+03, /* 0x45bad7c4 */
  9.6254453125e+03, /* 0x461665c8 */
  2.4060581055e+03, /* 0x451660ee */
};

static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
 -2.5470459075e-09, /* 0xb12f081b */
 -7.0311963558e-02, /* 0xbd8fffb8 */
 -2.4090321064e+00, /* 0xc01a2d95 */
 -2.1965976715e+01, /* 0xc1afba52 */
 -5.8079170227e+01, /* 0xc2685112 */
 -3.1447946548e+01, /* 0xc1fb9565 */
};
static const float pS3[5] = {
  3.5856033325e+01, /* 0x420f6c94 */
  3.6151397705e+02, /* 0x43b4c1ca */
  1.1936077881e+03, /* 0x44953373 */
  1.1279968262e+03, /* 0x448cffe6 */
  1.7358093262e+02, /* 0x432d94b8 */
};

static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
 -8.8753431271e-08, /* 0xb3be98b7 */
 -7.0303097367e-02, /* 0xbd8ffb12 */
 -1.4507384300e+00, /* 0xbfb9b1cc */
 -7.6356959343e+00, /* 0xc0f4579f */
 -1.1193166733e+01, /* 0xc1331736 */
 -3.2336456776e+00, /* 0xc04ef40d */
};
static const float pS2[5] = {
  2.2220300674e+01, /* 0x41b1c32d */
  1.3620678711e+02, /* 0x430834f0 */
  2.7047027588e+02, /* 0x43873c32 */
  1.5387539673e+02, /* 0x4319e01a */
  1.4657617569e+01, /* 0x416a859a */
};

static float pzerof(float x)
{
	const float *p,*q;
	float_t z,r,s;
	uint32_t ix;

	GET_FLOAT_WORD(ix, x);
	ix &= 0x7fffffff;
	if      (ix >= 0x41000000){p = pR8; q = pS8;}
	else if (ix >= 0x40f71c58){p = pR5; q = pS5;}
	else if (ix >= 0x4036db68){p = pR3; q = pS3;}
	else /*ix >= 0x40000000*/ {p = pR2; q = pS2;}
	z = 1.0f/(x*x);
	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
	s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
	return 1.0f + r/s;
}

/* For x >= 8, the asymptotic expansions of qzero is
 *      -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
 * We approximate pzero by
 *      qzero(x) = s*(-1.25 + (R/S))
 * where  R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
 *        S = 1 + qS0*s^2 + ... + qS5*s^12
 * and
 *      | qzero(x)/s +1.25-R/S | <= 2  ** ( -61.22)
 */
static const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
  0.0000000000e+00, /* 0x00000000 */
  7.3242187500e-02, /* 0x3d960000 */
  1.1768206596e+01, /* 0x413c4a93 */
  5.5767340088e+02, /* 0x440b6b19 */
  8.8591972656e+03, /* 0x460a6cca */
  3.7014625000e+04, /* 0x471096a0 */
};
static const float qS8[6] = {
  1.6377603149e+02, /* 0x4323c6aa */
  8.0983447266e+03, /* 0x45fd12c2 */
  1.4253829688e+05, /* 0x480b3293 */
  8.0330925000e+05, /* 0x49441ed4 */
  8.4050156250e+05, /* 0x494d3359 */
 -3.4389928125e+05, /* 0xc8a7eb69 */
};

static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
  1.8408595828e-11, /* 0x2da1ec79 */
  7.3242180049e-02, /* 0x3d95ffff */
  5.8356351852e+00, /* 0x40babd86 */
  1.3511157227e+02, /* 0x43071c90 */
  1.0272437744e+03, /* 0x448067cd */
  1.9899779053e+03, /* 0x44f8bf4b */
};
static const float qS5[6] = {
  8.2776611328e+01, /* 0x42a58da0 */
  2.0778142090e+03, /* 0x4501dd07 */
  1.8847289062e+04, /* 0x46933e94 */
  5.6751113281e+04, /* 0x475daf1d */
  3.5976753906e+04, /* 0x470c88c1 */
 -5.3543427734e+03, /* 0xc5a752be */
};

static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
  4.3774099900e-09, /* 0x3196681b */
  7.3241114616e-02, /* 0x3d95ff70 */
  3.3442313671e+00, /* 0x405607e3 */
  4.2621845245e+01, /* 0x422a7cc5 */
  1.7080809021e+02, /* 0x432acedf */
  1.6673394775e+02, /* 0x4326bbe4 */
};
static const float qS3[6] = {
  4.8758872986e+01, /* 0x42430916 */
  7.0968920898e+02, /* 0x44316c1c */
  3.7041481934e+03, /* 0x4567825f */
  6.4604252930e+03, /* 0x45c9e367 */
  2.5163337402e+03, /* 0x451d4557 */
 -1.4924745178e+02, /* 0xc3153f59 */
};

static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
  1.5044444979e-07, /* 0x342189db */
  7.3223426938e-02, /* 0x3d95f62a */
  1.9981917143e+00, /* 0x3fffc4bf */
  1.4495602608e+01, /* 0x4167edfd */
  3.1666231155e+01, /* 0x41fd5471 */
  1.6252708435e+01, /* 0x4182058c */
};
static const float qS2[6] = {
  3.0365585327e+01, /* 0x41f2ecb8 */
  2.6934811401e+02, /* 0x4386ac8f */
  8.4478375244e+02, /* 0x44533229 */
  8.8293585205e+02, /* 0x445cbbe5 */
  2.1266638184e+02, /* 0x4354aa98 */
 -5.3109550476e+00, /* 0xc0a9f358 */
};

static float qzerof(float x)
{
	const float *p,*q;
	float_t s,r,z;
	uint32_t ix;

	GET_FLOAT_WORD(ix, x);
	ix &= 0x7fffffff;
	if      (ix >= 0x41000000){p = qR8; q = qS8;}
	else if (ix >= 0x40f71c58){p = qR5; q = qS5;}
	else if (ix >= 0x4036db68){p = qR3; q = qS3;}
	else /*ix >= 0x40000000*/ {p = qR2; q = qS2;}
	z = 1.0f/(x*x);
	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
	s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
	return (-.125f + r/s)/x;
}


text:
#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double acoshl(long double x)
{
	return acosh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/* acosh(x) = log(x + sqrt(x*x-1)) */
long double acoshl(long double x)
{
	union ldshape u = {x};
	int e = u.i.se & 0x7fff;

	if (e < 0x3fff + 1)
		/* |x| < 2, invalid if x < 1 or nan */
		return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1)));
	if (e < 0x3fff + 32)
		/* |x| < 0x1p32 */
		return logl(2*x - 1/(x+sqrtl(x*x-1)));
	return logl(x) + 0.693147180559945309417232121458176568L;
}
#endif


text:
#include "libm.h"

/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */
float asinhf(float x)
{
	union {float f; uint32_t i;} u = {.f = x};
	uint32_t i = u.i & 0x7fffffff;
	unsigned s = u.i >> 31;

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

	if (i >= 0x3f800000 + (12<<23)) {
		/* |x| >= 0x1p12 or inf or nan */
		x = logf(x) + 0.693147180559945309417232121458176568f;
	} else if (i >= 0x3f800000 + (1<<23)) {
		/* |x| >= 2 */
		x = logf(2*x + 1/(sqrtf(x*x+1)+x));
	} else if (i >= 0x3f800000 - (12<<23)) {
		/* |x| >= 0x1p-12, up to 1.6ulp error in [0.125,0.5] */
		x = log1pf(x + x*x/(sqrtf(x*x+1)+1));
	} else {
		/* |x| < 0x1p-12, raise inexact if x!=0 */
		FORCE_EVAL(x + 0x1p120f);
	}
	return s ? -x : x;
}


text:
#include <math.h>

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


text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_log10f.c */
/*
 * ====================================================
 * 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.
 * ====================================================
 */
/*
 * See comments in log10.c.
 */

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

static const float
ivln10hi  =  4.3432617188e-01, /* 0x3ede6000 */
ivln10lo  = -3.1689971365e-05, /* 0xb804ead9 */
log10_2hi =  3.0102920532e-01, /* 0x3e9a2080 */
log10_2lo =  7.9034151668e-07, /* 0x355427db */
/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */

float log10f(float x)
{
	union {float f; uint32_t i;} u = {x};
	float_t hfsq,f,s,z,R,w,t1,t2,dk,hi,lo;
	uint32_t ix;
	int k;

	ix = u.i;
	k = 0;
	if (ix < 0x00800000 || ix>>31) {  /* x < 2**-126  */
		if (ix<<1 == 0)
			return -1/(x*x);  /* log(+-0)=-inf */
		if (ix>>31)
			return (x-x)/0.0f; /* log(-#) = NaN */
		/* subnormal number, scale up x */
		k -= 25;
		x *= 0x1p25f;
		u.f = x;
		ix = u.i;
	} else if (ix >= 0x7f800000) {
		return x;
	} else if (ix == 0x3f800000)
		return 0;

	/* reduce x into [sqrt(2)/2, sqrt(2)] */
	ix += 0x3f800000 - 0x3f3504f3;
	k += (int)(ix>>23) - 0x7f;
	ix = (ix&0x007fffff) + 0x3f3504f3;
	u.i = ix;
	x = u.f;

	f = x - 1.0f;
	s = f/(2.0f + f);
	z = s*s;
	w = z*z;
	t1= w*(Lg2+w*Lg4);
	t2= z*(Lg1+w*Lg3);
	R = t2 + t1;
	hfsq = 0.5f*f*f;

	hi = f - hfsq;
	u.f = hi;
	u.i &= 0xfffff000;
	hi = u.f;
	lo = f - hi - hfsq + s*(hfsq+R);
	dk = k;
	return dk*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi + hi*ivln10hi + dk*log10_2hi;
}


text:
#include <math.h>

long double sqrtl(long double x)
{
	/* FIXME: implement in C, this is for LDBL_MANT_DIG == 64 only */
	return sqrt(x);
}


text:
#include <math.h>
#include "libc.h"

int __signgam = 0;

weak_alias(__signgam, signgam);


text:
#define _GNU_SOURCE
#include <math.h>
#include <stdint.h>
#include "libc.h"

float exp10f(float x)
{
	static const float p10[] = {
		1e-7f, 1e-6f, 1e-5f, 1e-4f, 1e-3f, 1e-2f, 1e-1f,
		1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7
	};
	float n, y = modff(x, &n);
	union {float f; uint32_t i;} u = {n};
	/* fabsf(n) < 8 without raising invalid on nan */
	if ((u.i>>23 & 0xff) < 0x7f+3) {
		if (!y) return p10[(int)n+7];
		y = exp2f(3.32192809488736234787031942948939f * y);
		return y * p10[(int)n+7];
	}
	return exp2(3.32192809488736234787031942948939 * x);
}

weak_alias(exp10f, pow10f);


text:
/* origin: FreeBSD /usr/src/lib/msun/src/s_tan.c */
/*
 * ====================================================
 * 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.
 * ====================================================
 */
/* tan(x)
 * Return tangent function of x.
 *
 * kernel function:
 *      __tan           ... tangent function on [-pi/4,pi/4]
 *      __rem_pio2      ... argument reduction routine
 *
 * Method.
 *      Let S,C and T denote the sin, cos and tan respectively on
 *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
 *      in [-pi/4 , +pi/4], and let n = k mod 4.
 *      We have
 *
 *          n        sin(x)      cos(x)        tan(x)
 *     ----------------------------------------------------------
 *          0          S           C             T
 *          1          C          -S            -1/T
 *          2         -S          -C             T
 *          3         -C           S            -1/T
 *     ----------------------------------------------------------
 *
 * Special cases:
 *      Let trig be any of sin, cos, or tan.
 *      trig(+-INF)  is NaN, with signals;
 *      trig(NaN)    is that NaN;
 *
 * Accuracy:
 *      TRIG(x) returns trig(x) nearly rounded
 */

#include "libm.h"

double tan(double x)
{
	double y[2];
	uint32_t ix;
	unsigned n;

	GET_HIGH_WORD(ix, x);
	ix &= 0x7fffffff;

	/* |x| ~< pi/4 */
	if (ix <= 0x3fe921fb) {
		if (ix < 0x3e400000) { /* |x| < 2**-27 */
			/* raise inexact if x!=0 and underflow if subnormal */
			FORCE_EVAL(ix < 0x00100000 ? x/0x1p120f : x+0x1p120f);
			return x;
		}
		return __tan(x, 0.0, 0);
	}

	/* tan(Inf or NaN) is NaN */
	if (ix >= 0x7ff00000)
		return x - x;

	/* argument reduction */
	n = __rem_pio2(x, y);
	return __tan(y[0], y[1], n&1);
}


text:
/* origin: FreeBSD /usr/src/lib/msun/src/s_expm1.c */
/*
 * ====================================================
 * 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.
 * ====================================================
 */
/* expm1(x)
 * Returns exp(x)-1, the exponential of x minus 1.
 *
 * Method
 *   1. Argument reduction:
 *      Given x, find r and integer k such that
 *
 *               x = k*ln2 + r,  |r| <= 0.5*ln2 ~ 0.34658
 *
 *      Here a correction term c will be computed to compensate
 *      the error in r when rounded to a floating-point number.
 *
 *   2. Approximating expm1(r) by a special rational function on
 *      the interval [0,0.34658]:
 *      Since
 *          r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
 *      we define R1(r*r) by
 *          r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
 *      That is,
 *          R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
 *                   = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
 *                   = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
 *      We use a special Remez algorithm on [0,0.347] to generate
 *      a polynomial of degree 5 in r*r to approximate R1. The
 *      maximum error of this polynomial approximation is bounded
 *      by 2**-61. In other words,
 *          R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
 *      where   Q1  =  -1.6666666666666567384E-2,
 *              Q2  =   3.9682539681370365873E-4,
 *              Q3  =  -9.9206344733435987357E-6,
 *              Q4  =   2.5051361420808517002E-7,
 *              Q5  =  -6.2843505682382617102E-9;
 *              z   =  r*r,
 *      with error bounded by
 *          |                  5           |     -61
 *          | 1.0+Q1*z+...+Q5*z   -  R1(z) | <= 2
 *          |                              |
 *
 *      expm1(r) = exp(r)-1 is then computed by the following
 *      specific way which minimize the accumulation rounding error:
 *                             2     3
 *                            r     r    [ 3 - (R1 + R1*r/2)  ]
 *            expm1(r) = r + --- + --- * [--------------------]
 *                            2     2    [ 6 - r*(3 - R1*r/2) ]
 *
 *      To compensate the error in the argument reduction, we use
 *              expm1(r+c) = expm1(r) + c + expm1(r)*c
 *                         ~ expm1(r) + c + r*c
 *      Thus c+r*c will be added in as the correction terms for
 *      expm1(r+c). Now rearrange the term to avoid optimization
 *      screw up:
 *                      (      2                                    2 )
 *                      ({  ( r    [ R1 -  (3 - R1*r/2) ]  )  }    r  )
 *       expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
 *                      ({  ( 2    [ 6 - r*(3 - R1*r/2) ]  )  }    2  )
 *                      (                                             )
 *
 *                 = r - E
 *   3. Scale back to obtain expm1(x):
 *      From step 1, we have
 *         expm1(x) = either 2^k*[expm1(r)+1] - 1
 *                  = or     2^k*[expm1(r) + (1-2^-k)]
 *   4. Implementation notes:
 *      (A). To save one multiplication, we scale the coefficient Qi
 *           to Qi*2^i, and replace z by (x^2)/2.
 *      (B). To achieve maximum accuracy, we compute expm1(x) by
 *        (i)   if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
 *        (ii)  if k=0, return r-E
 *        (iii) if k=-1, return 0.5*(r-E)-0.5
 *        (iv)  if k=1 if r < -0.25, return 2*((r+0.5)- E)
 *                     else          return  1.0+2.0*(r-E);
 *        (v)   if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
 *        (vi)  if k <= 20, return 2^k((1-2^-k)-(E-r)), else
 *        (vii) return 2^k(1-((E+2^-k)-r))
 *
 * Special cases:
 *      expm1(INF) is INF, expm1(NaN) is NaN;
 *      expm1(-INF) is -1, and
 *      for finite argument, only expm1(0)=0 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 >  7.09782712893383973096e+02 then expm1(x) overflow
 *
 * Constants:
 * The hexadecimal values are the intended ones for the following
 * constants. The decimal values may be used, provided that the
 * compiler will convert from decimal to binary accurately enough
 * to produce the hexadecimal values shown.
 */

#include "libm.h"

static const double
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
ln2_hi      = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
ln2_lo      = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
invln2      = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
/* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */
Q1 = -3.33333333333331316428e-02, /* BFA11111 111110F4 */
Q2 =  1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */
Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
Q4 =  4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */

double expm1(double x)
{
	double_t y,hi,lo,c,t,e,hxs,hfx,r1,twopk;
	union {double f; uint64_t i;} u = {x};
	uint32_t hx = u.i>>32 & 0x7fffffff;
	int k, sign = u.i>>63;

	/* filter out huge and non-finite argument */
	if (hx >= 0x4043687A) {  /* if |x|>=56*ln2 */
		if (isnan(x))
			return x;
		if (sign)
			return -1;
		if (x > o_threshold) {
			x *= 0x1p1023;
			return x;
		}
	}

	/* argument reduction */
	if (hx > 0x3fd62e42) {  /* if  |x| > 0.5 ln2 */
		if (hx < 0x3FF0A2B2) {  /* and |x| < 1.5 ln2 */
			if (!sign) {
				hi = x - ln2_hi;
				lo = ln2_lo;
				k =  1;
			} else {
				hi = x + ln2_hi;
				lo = -ln2_lo;
				k = -1;
			}
		} else {
			k  = invln2*x + (sign ? -0.5 : 0.5);
			t  = k;
			hi = x - t*ln2_hi;  /* t*ln2_hi is exact here */
			lo = t*ln2_lo;
		}
		x = hi-lo;
		c = (hi-x)-lo;
	} else if (hx < 0x3c900000) {  /* |x| < 2**-54, return x */
		if (hx < 0x00100000)
			FORCE_EVAL((float)x);
		return x;
	} else
		k = 0;

	/* x is now in primary range */
	hfx = 0.5*x;
	hxs = x*hfx;
	r1 = 1.0+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
	t  = 3.0-r1*hfx;
	e  = hxs*((r1-t)/(6.0 - x*t));
	if (k == 0)   /* c is 0 */
		return x - (x*e-hxs);
	e  = x*(e-c) - c;
	e -= hxs;
	/* exp(x) ~ 2^k (x_reduced - e + 1) */
	if (k == -1)
		return 0.5*(x-e) - 0.5;
	if (k == 1) {
		if (x < -0.25)
			return -2.0*(e-(x+0.5));
		return 1.0+2.0*(x-e);
	}
	u.i = (uint64_t)(0x3ff + k)<<52;  /* 2^k */
	twopk = u.f;
	if (k < 0 || k > 56) {  /* suffice to return exp(x)-1 */
		y = x - e + 1.0;
		if (k == 1024)
			y = y*2.0*0x1p1023;
		else
			y = y*twopk;
		return y - 1.0;
	}
	u.i = (uint64_t)(0x3ff - k)<<52;  /* 2^-k */
	if (k < 20)
		y = (x-e+(1-u.f))*twopk;
	else
		y = (x-(e+u.f)+1)*twopk;
	return y;
}


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

float fabsf(float x)
{
	union {float f; uint32_t i;} u = {x};
	u.i &= 0x7fffffff;
	return u.f;
}


text:
#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double atanhl(long double x)
{
	return atanh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
long double atanhl(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 < 0x3ff - 1) {
		if (e < 0x3ff - LDBL_MANT_DIG/2) {
			/* handle underflow */
			if (e == 0)
				FORCE_EVAL((float)x);
		} else {
			/* |x| < 0.5, up to 1.7ulp error */
			x = 0.5*log1pl(2*x + 2*x*x/(1-x));
		}
	} else {
		/* avoid overflow */
		x = 0.5*log1pl(2*(x/(1-x)));
	}
	return s ? -x : x;
}
#endif


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

float fmodf(float x, float y)
{
	union {float f; uint32_t i;} ux = {x}, uy = {y};
	int ex = ux.i>>23 & 0xff;
	int ey = uy.i>>23 & 0xff;
	uint32_t sx = ux.i & 0x80000000;
	uint32_t i;
	uint32_t uxi = ux.i;

	if (uy.i<<1 == 0 || isnan(y) || ex == 0xff)
		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<<9; i>>31 == 0; ex--, i <<= 1);
		uxi <<= -ex + 1;
	} else {
		uxi &= -1U >> 9;
		uxi |= 1U << 23;
	}
	if (!ey) {
		for (i = uy.i<<9; i>>31 == 0; ey--, i <<= 1);
		uy.i <<= -ey + 1;
	} else {
		uy.i &= -1U >> 9;
		uy.i |= 1U << 23;
	}

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

	/* scale result up */
	if (ex > 0) {
		uxi -= 1U << 23;
		uxi |= (uint32_t)ex << 23;
	} else {
		uxi >>= -ex + 1;
	}
	uxi |= sx;
	ux.i = uxi;
	return ux.f;
}


text:
#include "libm.h"

#if FLT_EVAL_METHOD==2
#undef sqrtf
#define sqrtf sqrtl
#elif FLT_EVAL_METHOD==1
#undef sqrtf
#define sqrtf sqrt
#endif

/* acosh(x) = log(x + sqrt(x*x-1)) */
float acoshf(float x)
{
	union {float f; uint32_t i;} u = {x};
	uint32_t a = u.i & 0x7fffffff;

	if (a < 0x3f800000+(1<<23))
		/* |x| < 2, invalid if x < 1 or nan */
		/* up to 2ulp error in [1,1.125] */
		return log1pf(x-1 + sqrtf((x-1)*(x-1)+2*(x-1)));
	if (a < 0x3f800000+(12<<23))
		/* |x| < 0x1p12 */
		return logf(2*x - 1/(x+sqrtf(x*x-1)));
	/* x >= 0x1p12 */
	return logf(x) + 0.693147180559945309417232121458176568f;
}


text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.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.
 * ====================================================
 */
/* sqrt(x)
 * Return correctly rounded sqrt.
 *           ------------------------------------------
 *           |  Use the hardware sqrt if you have one |
 *           ------------------------------------------
 * Method:
 *   Bit by bit method using integer arithmetic. (Slow, but portable)
 *   1. Normalization
 *      Scale x to y in [1,4) with even powers of 2:
 *      find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
 *              sqrt(x) = 2^k * sqrt(y)
 *   2. Bit by bit computation
 *      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
 *           i                                                   0
 *                                     i+1         2
 *          s  = 2*q , and      y  =  2   * ( y - q  ).         (1)
 *           i      i            i                 i
 *
 *      To compute q    from q , one checks whether
 *                  i+1       i
 *
 *                            -(i+1) 2
 *                      (q + 2      ) <= y.                     (2)
 *                        i
 *                                                            -(i+1)
 *      If (2) is false, then q   = q ; otherwise q   = q  + 2      .
 *                             i+1   i             i+1   i
 *
 *      With some algebric manipulation, it is not difficult to see
 *      that (2) is equivalent to
 *                             -(i+1)
 *                      s  +  2       <= y                      (3)
 *                       i                i
 *
 *      The advantage of (3) is that s  and y  can be computed by
 *                                    i      i
 *      the following recurrence formula:
 *          if (3) is false
 *
 *          s     =  s  ,       y    = y   ;                    (4)
 *           i+1      i          i+1    i
 *
 *          otherwise,
 *                         -i                     -(i+1)
 *          s     =  s  + 2  ,  y    = y  -  s  - 2             (5)
 *           i+1      i          i+1    i     i
 *
 *      One may easily use induction to prove (4) and (5).
 *      Note. Since the left hand side of (3) contain only i+2 bits,
 *            it does not necessary to do a full (53-bit) comparison
 *            in (3).
 *   3. Final rounding
 *      After generating the 53 bits result, we compute one more bit.
 *      Together with the remainder, we can decide whether the
 *      result is exact, bigger than 1/2ulp, or less than 1/2ulp
 *      (it will never equal to 1/2ulp).
 *      The rounding mode can be detected by checking whether
 *      huge + tiny is equal to huge, and whether huge - tiny is
 *      equal to huge for some floating point number "huge" and "tiny".
 *
 * Special cases:
 *      sqrt(+-0) = +-0         ... exact
 *      sqrt(inf) = inf
 *      sqrt(-ve) = NaN         ... with invalid signal
 *      sqrt(NaN) = NaN         ... with invalid signal for signaling NaN
 */

#include "libm.h"

static const double tiny = 1.0e-300;

double sqrt(double x)
{
	double z;
	int32_t sign = (int)0x80000000;
	int32_t ix0,s0,q,m,t,i;
	uint32_t r,t1,s1,ix1,q1;

	EXTRACT_WORDS(ix0, ix1, x);

	/* take care of Inf and NaN */
	if ((ix0&0x7ff00000) == 0x7ff00000) {
		return x*x + x;  /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
	}
	/* take care of zero */
	if (ix0 <= 0) {
		if (((ix0&~sign)|ix1) == 0)
			return x;  /* sqrt(+-0) = +-0 */
		if (ix0 < 0)
			return (x-x)/(x-x);  /* sqrt(-ve) = sNaN */
	}
	/* normalize x */
	m = ix0>>20;
	if (m == 0) {  /* subnormal x */
		while (ix0 == 0) {
			m -= 21;
			ix0 |= (ix1>>11);
			ix1 <<= 21;
		}
		for (i=0; (ix0&0x00100000) == 0; i++)
			ix0<<=1;
		m -= i - 1;
		ix0 |= ix1>>(32-i);
		ix1 <<= i;
	}
	m -= 1023;    /* unbias exponent */
	ix0 = (ix0&0x000fffff)|0x00100000;
	if (m & 1) {  /* odd m, double x to make it even */
		ix0 += ix0 + ((ix1&sign)>>31);
		ix1 += ix1;
	}
	m >>= 1;      /* m = [m/2] */

	/* generate sqrt(x) bit by bit */
	ix0 += ix0 + ((ix1&sign)>>31);
	ix1 += ix1;
	q = q1 = s0 = s1 = 0;  /* [q,q1] = sqrt(x) */
	r = 0x00200000;        /* r = moving bit from right to left */

	while (r != 0) {
		t = s0 + r;
		if (t <= ix0) {
			s0   = t + r;
			ix0 -= t;
			q   += r;
		}
		ix0 += ix0 + ((ix1&sign)>>31);
		ix1 += ix1;
		r >>= 1;
	}

	r = sign;
	while (r != 0) {
		t1 = s1 + r;
		t  = s0;
		if (t < ix0 || (t == ix0 && t1 <= ix1)) {
			s1 = t1 + r;
			if ((t1&sign) == sign && (s1&sign) == 0)
				s0++;
			ix0 -= t;
			if (ix1 < t1)
				ix0--;
			ix1 -= t1;
			q1 += r;
		}
		ix0 += ix0 + ((ix1&sign)>>31);
		ix1 += ix1;
		r >>= 1;
	}

	/* use floating add to find out rounding direction */
	if ((ix0|ix1) != 0) {
		z = 1.0 - tiny; /* raise inexact flag */
		if (z >= 1.0) {
			z = 1.0 + tiny;
			if (q1 == (uint32_t)0xffffffff) {
				q1 = 0;
				q++;
			} else if (z > 1.0) {
				if (q1 == (uint32_t)0xfffffffe)
					q++;
				q1 += 2;
			} else
				q1 += q1 & 1;
		}
	}
	ix0 = (q>>1) + 0x3fe00000;
	ix1 = q1>>1;
	if (q&1)
		ix1 |= sign;
	ix0 += m << 20;
	INSERT_WORDS(z, ix0, ix1);
	return z;
}


text:
#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double nextafterl(long double x, long double y)
{
	return nextafter(x, y);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
long double nextafterl(long double x, long double y)
{
	union ldshape ux, uy;

	if (isnan(x) || isnan(y))
		return x + y;
	if (x == y)
		return y;
	ux.f = x;
	if (x == 0) {
		uy.f = y;
		ux.i.m = 1;
		ux.i.se = uy.i.se & 0x8000;
	} else if ((x < y) == !(ux.i.se & 0x8000)) {
		ux.i.m++;
		if (ux.i.m << 1 == 0) {
			ux.i.m = 1ULL << 63;
			ux.i.se++;
		}
	} else {
		if (ux.i.m << 1 == 0) {
			ux.i.se--;
			if (ux.i.se)
				ux.i.m = 0;
		}
		ux.i.m--;
	}
	/* raise overflow if ux is infinite and x is finite */
	if ((ux.i.se & 0x7fff) == 0x7fff)
		return x + x;
	/* raise underflow if ux is subnormal or zero */
	if ((ux.i.se & 0x7fff) == 0)
		FORCE_EVAL(x*x + ux.f*ux.f);
	return ux.f;
}
#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
long double nextafterl(long double x, long double y)
{
	union ldshape ux, uy;

	if (isnan(x) || isnan(y))
		return x + y;
	if (x == y)
		return y;
	ux.f = x;
	if (x == 0) {
		uy.f = y;
		ux.i.lo = 1;
		ux.i.se = uy.i.se & 0x8000;
	} else if ((x < y) == !(ux.i.se & 0x8000)) {
		ux.i2.lo++;
		if (ux.i2.lo == 0)
			ux.i2.hi++;
	} else {
		if (ux.i2.lo == 0)
			ux.i2.hi--;
		ux.i2.lo--;
	}
	/* raise overflow if ux is infinite and x is finite */
	if ((ux.i.se & 0x7fff) == 0x7fff)
		return x + x;
	/* raise underflow if ux is subnormal or zero */
	if ((ux.i.se & 0x7fff) == 0)
		FORCE_EVAL(x*x + ux.f*ux.f);
	return ux.f;
}
#endif


text:
#define _GNU_SOURCE
#include <float.h>
#include <math.h>
#include "libc.h"
#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double exp10l(long double x)
{
	return exp10(x);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
long double exp10l(long double x)
{
	static const long double p10[] = {
		1e-15L, 1e-14L, 1e-13L, 1e-12L, 1e-11L, 1e-10L,
		1e-9L, 1e-8L, 1e-7L, 1e-6L, 1e-5L, 1e-4L, 1e-3L, 1e-2L, 1e-1L,
		1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
		1e10, 1e11, 1e12, 1e13, 1e14, 1e15
	};
	long double n, y = modfl(x, &n);
	union ldshape u = {n};
	/* fabsl(n) < 16 without raising invalid on nan */
	if ((u.i.se & 0x7fff) < 0x3fff+4) {
		if (!y) return p10[(int)n+15];
		y = exp2l(3.32192809488736234787031942948939L * y);
		return y * p10[(int)n+15];
	}
	return powl(10.0, x);
}
#endif

weak_alias(exp10l, pow10l);


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

/* nearbyint is the same as rint, but it must not raise the inexact exception */

double nearbyint(double x)
{
#ifdef FE_INEXACT
	#pragma STDC FENV_ACCESS ON
	int e;

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


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

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long long llrintl(long double x)
{
	return llrint(x);
}
#elif defined(FE_INEXACT)
/*
see comments in lrint.c

Note that if LLONG_MAX == 0x7fffffffffffffff && LDBL_MANT_DIG == 64
then x == 2**63 - 0.5 is the only input that overflows and
raises inexact (with tonearest or upward rounding mode)
*/
long long llrintl(long double x)
{
	#pragma STDC FENV_ACCESS ON
	int e;

	e = fetestexcept(FE_INEXACT);
	x = rintl(x);
	if (!e && (x > LLONG_MAX || x < LLONG_MIN))
		feclearexcept(FE_INEXACT);
	/* conversion */
	return x;
}
#else
long long llrintl(long double x)
{
	return rintl(x);
}
#endif


text:
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_expl.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.
 */
/*
 *      Exponential function, long double precision
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, expl();
 *
 * y = expl( x );
 *
 *
 * DESCRIPTION:
 *
 * Returns e (2.71828...) raised to the x power.
 *
 * Range reduction is accomplished by separating the argument
 * into an integer k and fraction f such that
 *
 *     x    k  f
 *    e  = 2  e.
 *
 * A Pade' form of degree 5/6 is used to approximate exp(f) - 1
 * in the basic range [-0.5 ln 2, 0.5 ln 2].
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      +-10000     50000       1.12e-19    2.81e-20
 *
 *
 * Error amplification in the exponential function can be
 * a serious matter.  The error propagation involves
 * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
 * which shows that a 1 lsb error in representing X produces
 * a relative error of X times 1 lsb in the function.
 * While the routine gives an accurate result for arguments
 * that are exactly represented by a long double precision
 * computer number, the result contains amplified roundoff
 * error for large arguments not exactly represented.
 *
 *
 * ERROR MESSAGES:
 *
 *   message         condition      value returned
 * exp underflow    x < MINLOG         0.0
 * exp overflow     x > MAXLOG         MAXNUM
 *
 */

#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double expl(long double x)
{
	return exp(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384

static const long double P[3] = {
 1.2617719307481059087798E-4L,
 3.0299440770744196129956E-2L,
 9.9999999999999999991025E-1L,
};
static const long double Q[4] = {
 3.0019850513866445504159E-6L,
 2.5244834034968410419224E-3L,
 2.2726554820815502876593E-1L,
 2.0000000000000000000897E0L,
};
static const long double
LN2HI = 6.9314575195312500000000E-1L,
LN2LO = 1.4286068203094172321215E-6L,
LOG2E = 1.4426950408889634073599E0L;

long double expl(long double x)
{
	long double px, xx;
	int k;

	if (isnan(x))
		return x;
	if (x > 11356.5234062941439488L) /* x > ln(2^16384 - 0.5) */
		return x * 0x1p16383L;
	if (x < -11399.4985314888605581L) /* x < ln(2^-16446) */
		return -0x1p-16445L/x;

	/* Express e**x = e**f 2**k
	 *   = e**(f + k ln(2))
	 */
	px = floorl(LOG2E * x + 0.5);
	k = px;
	x -= px * LN2HI;
	x -= px * LN2LO;

	/* rational approximation of the fractional part:
	 * e**x =  1 + 2x P(x**2)/(Q(x**2) - x P(x**2))
	 */
	xx = x * x;
	px = x * __polevll(xx, P, 2);
	x = px/(__polevll(xx, Q, 3) - px);
	x = 1.0 + 2.0 * x;
	return scalbnl(x, k);
}
#endif


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

float copysignf(float x, float y)
{
	union {float f; uint32_t i;} ux={x}, uy={y};
	ux.i &= 0x7fffffff;
	ux.i |= uy.i & 0x80000000;
	return ux.f;
}


text:
/*
"A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964)
"Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001)
"An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004)

approximation method:

                        (x - 0.5)         S(x)
Gamma(x) = (x + g - 0.5)         *  ----------------
                                    exp(x + g - 0.5)

with
                 a1      a2      a3            aN
S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ]
               x + 1   x + 2   x + 3         x + N

with a0, a1, a2, a3,.. aN constants which depend on g.

for x < 0 the following reflection formula is used:

Gamma(x)*Gamma(-x) = -pi/(x sin(pi x))

most ideas and constants are from boost and python
*/
#include "libm.h"

static const double pi = 3.141592653589793238462643383279502884;

/* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */
static double sinpi(double x)
{
	int n;

	/* argument reduction: x = |x| mod 2 */
	/* spurious inexact when x is odd int */
	x = x * 0.5;
	x = 2 * (x - floor(x));

	/* reduce x into [-.25,.25] */
	n = 4 * x;
	n = (n+1)/2;
	x -= n * 0.5;

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

#define N 12
//static const double g = 6.024680040776729583740234375;
static const double gmhalf = 5.524680040776729583740234375;
static const double Snum[N+1] = {
	23531376880.410759688572007674451636754734846804940,
	42919803642.649098768957899047001988850926355848959,
	35711959237.355668049440185451547166705960488635843,
	17921034426.037209699919755754458931112671403265390,
	6039542586.3520280050642916443072979210699388420708,
	1439720407.3117216736632230727949123939715485786772,
	248874557.86205415651146038641322942321632125127801,
	31426415.585400194380614231628318205362874684987640,
	2876370.6289353724412254090516208496135991145378768,
	186056.26539522349504029498971604569928220784236328,
	8071.6720023658162106380029022722506138218516325024,
	210.82427775157934587250973392071336271166969580291,
	2.5066282746310002701649081771338373386264310793408,
};
static const double Sden[N+1] = {
	0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535,
	2637558, 357423, 32670, 1925, 66, 1,
};
/* n! for small integer n */
static const double fact[] = {
	1, 1, 2, 6, 24, 120, 720, 5040.0, 40320.0, 362880.0, 3628800.0, 39916800.0,
	479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0, 20922789888000.0,
	355687428096000.0, 6402373705728000.0, 121645100408832000.0,
	2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0,
};

/* S(x) rational function for positive x */
static double S(double x)
{
	double_t num = 0, den = 0;
	int i;

	/* to avoid overflow handle large x differently */
	if (x < 8)
		for (i = N; i >= 0; i--) {
			num = num * x + Snum[i];
			den = den * x + Sden[i];
		}
	else
		for (i = 0; i <= N; i++) {
			num = num / x + Snum[i];
			den = den / x + Sden[i];
		}
	return num/den;
}

double tgamma(double x)
{
	union {double f; uint64_t i;} u = {x};
	double absx, y;
	double_t dy, z, r;
	uint32_t ix = u.i>>32 & 0x7fffffff;
	int sign = u.i>>63;

	/* special cases */
	if (ix >= 0x7ff00000)
		/* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
		return x + INFINITY;
	if (ix < (0x3ff-54)<<20)
		/* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */
		return 1/x;

	/* integer arguments */
	/* raise inexact when non-integer */
	if (x == floor(x)) {
		if (sign)
			return 0/0.0;
		if (x <= sizeof fact/sizeof *fact)
			return fact[(int)x - 1];
	}

	/* x >= 172: tgamma(x)=inf with overflow */
	/* x =< -184: tgamma(x)=+-0 with underflow */
	if (ix >= 0x40670000) { /* |x| >= 184 */
		if (sign) {
			FORCE_EVAL((float)(0x1p-126/x));
			if (floor(x) * 0.5 == floor(x * 0.5))
				return 0;
			return -0.0;
		}
		x *= 0x1p1023;
		return x;
	}

	absx = sign ? -x : x;

	/* handle the error of x + g - 0.5 */
	y = absx + gmhalf;
	if (absx > gmhalf) {
		dy = y - absx;
		dy -= gmhalf;
	} else {
		dy = y - gmhalf;
		dy -= absx;
	}

	z = absx - 0.5;
	r = S(absx) * exp(-y);
	if (x < 0) {
		/* reflection formula for negative x */
		/* sinpi(absx) is not 0, integers are already handled */
		r = -pi / (sinpi(absx) * absx * r);
		dy = -dy;
		z = -z;
	}
	r += dy * (gmhalf+0.5) * r / y;
	z = pow(y, 0.5*z);
	y = r * z * z;
	return y;
}

#if 0
double __lgamma_r(double x, int *sign)
{
	double r, absx;

	*sign = 1;

	/* special cases */
	if (!isfinite(x))
		/* lgamma(nan)=nan, lgamma(+-inf)=inf */
		return x*x;

	/* integer arguments */
	if (x == floor(x) && x <= 2) {
		/* n <= 0: lgamma(n)=inf with divbyzero */
		/* n == 1,2: lgamma(n)=0 */
		if (x <= 0)
			return 1/0.0;
		return 0;
	}

	absx = fabs(x);

	/* lgamma(x) ~ -log(|x|) for tiny |x| */
	if (absx < 0x1p-54) {
		*sign = 1 - 2*!!signbit(x);
		return -log(absx);
	}

	/* use tgamma for smaller |x| */
	if (absx < 128) {
		x = tgamma(x);
		*sign = 1 - 2*!!signbit(x);
		return log(fabs(x));
	}

	/* second term (log(S)-g) could be more precise here.. */
	/* or with stirling: (|x|-0.5)*(log(|x|)-1) + poly(1/|x|) */
	r = (absx-0.5)*(log(absx+gmhalf)-1) + (log(S(absx)) - (gmhalf+0.5));
	if (x < 0) {
		/* reflection formula for negative x */
		x = sinpi(absx);
		*sign = 2*!!signbit(x) - 1;
		r = log(pi/(fabs(x)*absx)) - r;
	}
	return r;
}

weak_alias(__lgamma_r, lgamma_r);
#endif


text:
/* origin: FreeBSD /usr/src/lib/msun/src/s_tanf.c */
/*
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
 * 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.
 * ====================================================
 */

#include "libm.h"

/* Small multiples of pi/2 rounded to double precision. */
static const double
t1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
t2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
t3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
t4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */

float tanf(float x)
{
	double y;
	uint32_t ix;
	unsigned n, sign;

	GET_FLOAT_WORD(ix, x);
	sign = ix >> 31;
	ix &= 0x7fffffff;

	if (ix <= 0x3f490fda) {  /* |x| ~<= pi/4 */
		if (ix < 0x39800000) {  /* |x| < 2**-12 */
			/* raise inexact if x!=0 and underflow if subnormal */
			FORCE_EVAL(ix < 0x00800000 ? x/0x1p120f : x+0x1p120f);
			return x;
		}
		return __tandf(x, 0);
	}
	if (ix <= 0x407b53d1) {  /* |x| ~<= 5*pi/4 */
		if (ix <= 0x4016cbe3)  /* |x| ~<= 3pi/4 */
			return __tandf((sign ? x+t1pio2 : x-t1pio2), 1);
		else
			return __tandf((sign ? x+t2pio2 : x-t2pio2), 0);
	}
	if (ix <= 0x40e231d5) {  /* |x| ~<= 9*pi/4 */
		if (ix <= 0x40afeddf)  /* |x| ~<= 7*pi/4 */
			return __tandf((sign ? x+t3pio2 : x-t3pio2), 1);
		else
			return __tandf((sign ? x+t4pio2 : x-t4pio2), 0);
	}

	/* tan(Inf or NaN) is NaN */
	if (ix >= 0x7f800000)
		return x - x;

	/* argument reduction */
	n = __rem_pio2f(x, &y);
	return __tandf(y, n&1);
}


text:
#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double copysignl(long double x, long double y)
{
	return copysign(x, y);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
long double copysignl(long double x, long double y)
{
	union ldshape ux = {x}, uy = {y};
	ux.i.se &= 0x7fff;
	ux.i.se |= uy.i.se & 0x8000;
	return ux.f;
}
#endif


text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_lgammaf_r.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"
#include "libc.h"

static const float
pi  =  3.1415927410e+00, /* 0x40490fdb */
a0  =  7.7215664089e-02, /* 0x3d9e233f */
a1  =  3.2246702909e-01, /* 0x3ea51a66 */
a2  =  6.7352302372e-02, /* 0x3d89f001 */
a3  =  2.0580807701e-02, /* 0x3ca89915 */
a4  =  7.3855509982e-03, /* 0x3bf2027e */
a5  =  2.8905137442e-03, /* 0x3b3d6ec6 */
a6  =  1.1927076848e-03, /* 0x3a9c54a1 */
a7  =  5.1006977446e-04, /* 0x3a05b634 */
a8  =  2.2086278477e-04, /* 0x39679767 */
a9  =  1.0801156895e-04, /* 0x38e28445 */
a10 =  2.5214456400e-05, /* 0x37d383a2 */
a11 =  4.4864096708e-05, /* 0x383c2c75 */
tc  =  1.4616321325e+00, /* 0x3fbb16c3 */
tf  = -1.2148628384e-01, /* 0xbdf8cdcd */
/* tt = -(tail of tf) */
tt  =  6.6971006518e-09, /* 0x31e61c52 */
t0  =  4.8383611441e-01, /* 0x3ef7b95e */
t1  = -1.4758771658e-01, /* 0xbe17213c */
t2  =  6.4624942839e-02, /* 0x3d845a15 */
t3  = -3.2788541168e-02, /* 0xbd064d47 */
t4  =  1.7970675603e-02, /* 0x3c93373d */
t5  = -1.0314224288e-02, /* 0xbc28fcfe */
t6  =  6.1005386524e-03, /* 0x3bc7e707 */
t7  = -3.6845202558e-03, /* 0xbb7177fe */
t8  =  2.2596477065e-03, /* 0x3b141699 */
t9  = -1.4034647029e-03, /* 0xbab7f476 */
t10 =  8.8108185446e-04, /* 0x3a66f867 */
t11 = -5.3859531181e-04, /* 0xba0d3085 */
t12 =  3.1563205994e-04, /* 0x39a57b6b */
t13 = -3.1275415677e-04, /* 0xb9a3f927 */
t14 =  3.3552918467e-04, /* 0x39afe9f7 */
u0  = -7.7215664089e-02, /* 0xbd9e233f */
u1  =  6.3282704353e-01, /* 0x3f2200f4 */
u2  =  1.4549225569e+00, /* 0x3fba3ae7 */
u3  =  9.7771751881e-01, /* 0x3f7a4bb2 */
u4  =  2.2896373272e-01, /* 0x3e6a7578 */
u5  =  1.3381091878e-02, /* 0x3c5b3c5e */
v1  =  2.4559779167e+00, /* 0x401d2ebe */
v2  =  2.1284897327e+00, /* 0x4008392d */
v3  =  7.6928514242e-01, /* 0x3f44efdf */
v4  =  1.0422264785e-01, /* 0x3dd572af */
v5  =  3.2170924824e-03, /* 0x3b52d5db */
s0  = -7.7215664089e-02, /* 0xbd9e233f */
s1  =  2.1498242021e-01, /* 0x3e5c245a */
s2  =  3.2577878237e-01, /* 0x3ea6cc7a */
s3  =  1.4635047317e-01, /* 0x3e15dce6 */
s4  =  2.6642270386e-02, /* 0x3cda40e4 */
s5  =  1.8402845599e-03, /* 0x3af135b4 */
s6  =  3.1947532989e-05, /* 0x3805ff67 */
r1  =  1.3920053244e+00, /* 0x3fb22d3b */
r2  =  7.2193557024e-01, /* 0x3f38d0c5 */
r3  =  1.7193385959e-01, /* 0x3e300f6e */
r4  =  1.8645919859e-02, /* 0x3c98bf54 */
r5  =  7.7794247773e-04, /* 0x3a4beed6 */
r6  =  7.3266842264e-06, /* 0x36f5d7bd */
w0  =  4.1893854737e-01, /* 0x3ed67f1d */
w1  =  8.3333335817e-02, /* 0x3daaaaab */
w2  = -2.7777778450e-03, /* 0xbb360b61 */
w3  =  7.9365057172e-04, /* 0x3a500cfd */
w4  = -5.9518753551e-04, /* 0xba1c065c */
w5  =  8.3633989561e-04, /* 0x3a5b3dd2 */
w6  = -1.6309292987e-03; /* 0xbad5c4e8 */

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

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

	n = (int)(x*4);
	n = (n+1)/2;
	y = x - n*0.5f;
	y *= 3.14159265358979323846;
	switch (n) {
	default: /* case 4: */
	case 0: return __sindf(y);
	case 1: return __cosdf(y);
	case 2: return __sindf(-y);
	case 3: return -__cosdf(y);
	}
}

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

	/* purge off +-inf, NaN, +-0, tiny and negative arguments */
	*signgamp = 1;
	sign = u.i>>31;
	ix = u.i & 0x7fffffff;
	if (ix >= 0x7f800000)
		return x*x;
	if (ix < 0x35000000) {  /* |x| < 2**-21, return -log(|x|) */
		if (sign) {
			*signgamp = -1;
			x = -x;
		}
		return -logf(x);
	}
	if (sign) {
		x = -x;
		t = sin_pi(x);
		if (t == 0.0f) /* -integer */
			return 1.0f/(x-x);
		if (t > 0.0f)
			*signgamp = -1;
		else
			t = -t;
		nadj = logf(pi/(t*x));
	}

	/* purge off 1 and 2 */
	if (ix == 0x3f800000 || ix == 0x40000000)
		r = 0;
	/* for x < 2.0 */
	else if (ix < 0x40000000) {
		if (ix <= 0x3f666666) {  /* lgamma(x) = lgamma(x+1)-log(x) */
			r = -logf(x);
			if (ix >= 0x3f3b4a20) {
				y = 1.0f - x;
				i = 0;
			} else if (ix >= 0x3e6d3308) {
				y = x - (tc-1.0f);
				i = 1;
			} else {
				y = x;
				i = 2;
			}
		} else {
			r = 0.0f;
			if (ix >= 0x3fdda618) {  /* [1.7316,2] */
				y = 2.0f - x;
				i = 0;
			} else if (ix >= 0x3F9da620) {  /* [1.23,1.73] */
				y = x - tc;
				i = 1;
			} else {
				y = x - 1.0f;
				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.5f*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.0f+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
			r += -0.5f*y + p1/p2;
		}
	} else if (ix < 0x41000000) {  /* x < 8.0 */
		i = (int)x;
		y = x - (float)i;
		p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
		q = 1.0f+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
		r = 0.5f*y+p/q;
		z = 1.0f;    /* lgamma(1+s) = log(s) + lgamma(s) */
		switch (i) {
		case 7: z *= y + 6.0f;  /* FALLTHRU */
		case 6: z *= y + 5.0f;  /* FALLTHRU */
		case 5: z *= y + 4.0f;  /* FALLTHRU */
		case 4: z *= y + 3.0f;  /* FALLTHRU */
		case 3: z *= y + 2.0f;  /* FALLTHRU */
			r += logf(z);
			break;
		}
	} else if (ix < 0x5c800000) {  /* 8.0 <= x < 2**58 */
		t = logf(x);
		z = 1.0f/x;
		y = z*z;
		w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
		r = (x-0.5f)*(t-1.0f)+w;
	} else                         /* 2**58 <= x <= inf */
		r =  x*(logf(x)-1.0f);
	if (sign)
		r = nadj - r;
	return r;
}

weak_alias(__lgammaf_r, lgammaf_r);


text:
/* origin: FreeBSD /usr/src/lib/msun/src/e_atan2f.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
pi     = 3.1415927410e+00, /* 0x40490fdb */
pi_lo  = -8.7422776573e-08; /* 0xb3bbbd2e */

float atan2f(float y, float x)
{
	float z;
	uint32_t m,ix,iy;

	if (isnan(x) || isnan(y))
		return x+y;
	GET_FLOAT_WORD(ix, x);
	GET_FLOAT_WORD(iy, y);
	if (ix == 0x3f800000)  /* x=1.0 */
		return atanf(y);
	m = ((iy>>31)&1) | ((ix>>30)&2);  /* 2*sign(x)+sign(y) */
	ix &= 0x7fffffff;
	iy &= 0x7fffffff;

	/* when y = 0 */
	if (iy == 0) {
		switch (m) {
		case 0:
		case 1: return y;   /* atan(+-0,+anything)=+-0 */
		case 2: return  pi; /* atan(+0,-anything) = pi */
		case 3: return -pi; /* atan(-0,-anything) =-pi */
		}
	}
	/* when x = 0 */
	if (ix == 0)
		return m&1 ? -pi/2 : pi/2;
	/* when x is INF */
	if (ix == 0x7f800000) {
		if (iy == 0x7f800000) {
			switch (m) {
			case 0: return  pi/4; /* atan(+INF,+INF) */
			case 1: return -pi/4; /* atan(-INF,+INF) */
			case 2: return 3*pi/4;  /*atan(+INF,-INF)*/
			case 3: return -3*pi/4; /*atan(-INF,-INF)*/
			}
		} else {
			switch (m) {
			case 0: return  0.0f;    /* atan(+...,+INF) */
			case 1: return -0.0f;    /* atan(-...,+INF) */
			case 2: return  pi; /* atan(+...,-INF) */
			case 3: return -pi; /* atan(-...,-INF) */
			}
		}
	}
	/* |y/x| > 0x1p26 */
	if (ix+(26<<23) < iy || iy == 0x7f800000)
		return m&1 ? -pi/2 : pi/2;

	/* z = atan(|y/x|) with correct underflow */
	if ((m&2) && iy+(26<<23) < ix)  /*|y/x| < 0x1p-26, x < 0 */
		z = 0.0;
	else
		z = atanf(fabsf(y/x));
	switch (m) {
	case 0: return z;              /* atan(+,+) */
	case 1: return -z;             /* atan(-,+) */
	case 2: return pi - (z-pi_lo); /* atan(+,-) */
	default: /* case 3 */
		return (z-pi_lo) - pi; /* atan(-,-) */
	}
}


text:
#include "libm.h"

/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
float atanhf(float x)
{
	union {float f; uint32_t i;} u = {.f = x};
	unsigned s = u.i >> 31;
	float_t y;

	/* |x| */
	u.i &= 0x7fffffff;
	y = u.f;

	if (u.i < 0x3f800000 - (1<<23)) {
		if (u.i < 0x3f800000 - (32<<23)) {
			/* handle underflow */
			if (u.i < (1<<23))
				FORCE_EVAL((float)(y*y));
		} else {
			/* |x| < 0.5, up to 1.7ulp error */
			y = 0.5f*log1pf(2*y + 2*y*y/(1-y));
		}
	} else {
		/* avoid overflow */
		y = 0.5f*log1pf(2*(y/(1-y)));
	}
	return s ? -y : y;
}


text:
#include "libm.h"

double round(double x)
{
	union {double f; uint64_t i;} u = {x};
	int e = u.i >> 52 & 0x7ff;
	double_t y;

	if (e >= 0x3ff+52)
		return x;
	if (u.i >> 63)
		x = -x;
	if (e < 0x3ff-1) {
		/* raise inexact if x!=0 */
		FORCE_EVAL(x + 0x1p52);
		return 0*u.f;
	}
	y = (double)(x + 0x1p52) - 0x1p52 - x;
	if (y > 0.5)
		y = y + x - 1;
	else if (y <= -0.5)
		y = y + x + 1;
	else
		y = y + x;
	if (u.i >> 63)
		y = -y;
	return y;
}


text:
#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
double nexttoward(double x, long double y)
{
	return nextafter(x, y);
}
#else
double nexttoward(double x, long double y)
{
	union {double f; uint64_t i;} ux = {x};
	int e;

	if (isnan(x) || isnan(y))
		return x + y;
	if (x == y)
		return y;
	if (x == 0) {
		ux.i = 1;
		if (signbit(y))
			ux.i |= 1ULL<<63;
	} else if (x < y) {
		if (signbit(x))
			ux.i--;
		else
			ux.i++;
	} else {
		if (signbit(x))
			ux.i++;
		else
			ux.i--;
	}
	e = ux.i>>52 & 0x7ff;
	/* raise overflow if ux.f is infinite and x is finite */
	if (e == 0x7ff)
		FORCE_EVAL(x+x);
	/* raise underflow if ux.f is subnormal or zero */
	if (e == 0)
		FORCE_EVAL(x*x + ux.f*ux.f);
	return ux.f;
}
#endif


text:
/* origin: FreeBSD /usr/src/lib/msun/src/k_tanf.c */
/*
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
 * Optimized by Bruce D. Evans.
 */
/*
 * ====================================================
 * Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
 *
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */

#include "libm.h"

/* |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]). */
static const double T[] = {
  0x15554d3418c99f.0p-54, /* 0.333331395030791399758 */
  0x1112fd38999f72.0p-55, /* 0.133392002712976742718 */
  0x1b54c91d865afe.0p-57, /* 0.0533812378445670393523 */
  0x191df3908c33ce.0p-58, /* 0.0245283181166547278873 */
  0x185dadfcecf44e.0p-61, /* 0.00297435743359967304927 */
  0x1362b9bf971bcd.0p-59, /* 0.00946564784943673166728 */
};

float __tandf(double x, int odd)
{
	double_t z,r,w,s,t,u;

	z = x*x;
	/*
	 * Split up the polynomial into small independent terms to give
	 * opportunities for parallel evaluation.  The chosen splitting is
	 * micro-optimized for Athlons (XP, X64).  It costs 2 multiplications
	 * relative to Horner's method on sequential machines.
	 *
	 * We add the small terms from lowest degree up for efficiency on
	 * non-sequential machines (the lowest degree terms tend to be ready
	 * earlier).  Apart from this, we don't care about order of
	 * operations, and don't need to to care since we have precision to
	 * spare.  However, the chosen splitting is good for accuracy too,
	 * and would give results as accurate as Horner's method if the
	 * small terms were added from highest degree down.
	 */
	r = T[4] + z*T[5];
	t = T[2] + z*T[3];
	w = z*z;
	s = z*x;
	u = T[0] + z*T[1];
	r = (x + s*u) + (s*w)*(t + w*r);
	return odd ? -1.0/r : r;
}


text:
#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double sinl(long double x)
{
	return sin(x);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
long double sinl(long double x)
{
	union ldshape u = {x};
	unsigned n;
	long double y[2], hi, lo;

	u.i.se &= 0x7fff;
	if (u.i.se == 0x7fff)
		return x - x;
	if (u.f < M_PI_4) {
		if (u.i.se < 0x3fff - LDBL_MANT_DIG/2) {
			/* raise inexact if x!=0 and underflow if subnormal */
			FORCE_EVAL(u.i.se == 0 ? x*0x1p-120f : x+0x1p120f);
			return x;
		}
		return __sinl(x, 0.0, 0);
	}
	n = __rem_pio2l(x, y);
	hi = y[0];
	lo = y[1];
	switch (n & 3) {
	case 0:
		return __sinl(hi, lo, 1);
	case 1:
		return __cosl(hi, lo);
	case 2:
		return -__sinl(hi, lo, 1);
	case 3:
	default:
		return -__cosl(hi, lo);
	}
}
#endif


text:
#include <math.h>

long double ldexpl(long double x, int n)
{
	return scalbnl(x, n);
}


text:
/* origin: FreeBSD /usr/src/lib/msun/src/s_log1pf.c */
/*
 * ====================================================
 * 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
ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */

float log1pf(float x)
{
	union {float f; uint32_t i;} u = {x};
	float_t hfsq,f,c,s,z,R,w,t1,t2,dk;
	uint32_t ix,iu;
	int k;

	ix = u.i;
	k = 1;
	if (ix < 0x3ed413d0 || ix>>31) {  /* 1+x < sqrt(2)+  */
		if (ix >= 0xbf800000) {  /* x <= -1.0 */
			if (x == -1)
				return x/0.0f; /* log1p(-1)=+inf */
			return (x-x)/0.0f;     /* log1p(x<-1)=NaN */
		}
		if (ix<<1 < 0x33800000<<1) {   /* |x| < 2**-24 */
			/* underflow if subnormal */
			if ((ix&0x7f800000) == 0)
				FORCE_EVAL(x*x);
			return x;
		}
		if (ix <= 0xbe95f619) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
			k = 0;
			c = 0;
			f = x;
		}
	} else if (ix >= 0x7f800000)
		return x;
	if (k) {
		u.f = 1 + x;
		iu = u.i;
		iu += 0x3f800000 - 0x3f3504f3;
		k = (int)(iu>>23) - 0x7f;
		/* correction term ~ log(1+x)-log(u), avoid underflow in c/u */
		if (k < 25) {
			c = k >= 2 ? 1-(u.f-x) : x-(u.f-1);
			c /= u.f;
		} else
			c = 0;
		/* reduce u into [sqrt(2)/2, sqrt(2)] */
		iu = (iu&0x007fffff) + 0x3f3504f3;
		u.i = iu;
		f = u.f - 1;
	}
	s = f/(2.0f + f);
	z = s*s;
	w = z*z;
	t1= w*(Lg2+w*Lg4);
	t2= z*(Lg1+w*Lg3);
	R = t2 + t1;
	hfsq = 0.5f*f*f;
	dk = k;
	return s*(hfsq+R) + (dk*ln2_lo+c) - hfsq + f + dk*ln2_hi;
}


text:
/* origin: FreeBSD /usr/src/lib/msun/src/s_exp2f.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 16

static const float
redux = 0x1.8p23f / TBLSIZE,
P1    = 0x1.62e430p-1f,
P2    = 0x1.ebfbe0p-3f,
P3    = 0x1.c6b348p-5f,
P4    = 0x1.3b2c9cp-7f;

static const double exp2ft[TBLSIZE] = {
  0x1.6a09e667f3bcdp-1,
  0x1.7a11473eb0187p-1,
  0x1.8ace5422aa0dbp-1,
  0x1.9c49182a3f090p-1,
  0x1.ae89f995ad3adp-1,
  0x1.c199bdd85529cp-1,
  0x1.d5818dcfba487p-1,
  0x1.ea4afa2a490dap-1,
  0x1.0000000000000p+0,
  0x1.0b5586cf9890fp+0,
  0x1.172b83c7d517bp+0,
  0x1.2387a6e756238p+0,
  0x1.306fe0a31b715p+0,
  0x1.3dea64c123422p+0,
  0x1.4bfdad5362a27p+0,
  0x1.5ab07dd485429p+0,
};

/*
 * exp2f(x): compute the base 2 exponential of x
 *
 * Accuracy: Peak error < 0.501 ulp; location of peak: -0.030110927.
 *
 * Method: (equally-spaced tables)
 *
 *   Reduce x:
 *     x = k + y, for integer k and |y| <= 1/2.
 *     Thus we have exp2f(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**-(TBLSIZE+1).
 *
 *   We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a
 *   degree-4 minimax polynomial with maximum error under 1.4 * 2**-33.
 *   Using double precision for everything except the reduction makes
 *   roundoff error insignificant and simplifies the scaling step.
 *
 *   This method is due to Tang, but I do not use his suggested parameters:
 *
 *      Tang, P.  Table-driven Implementation of the Exponential Function
 *      in IEEE Floating-Point Arithmetic.  TOMS 15(2), 144-157 (1989).
 */
float exp2f(float x)
{
	double_t t, r, z;
	union {float f; uint32_t i;} u = {x};
	union {double f; uint64_t i;} uk;
	uint32_t ix, i0, k;

	/* Filter out exceptional cases. */
	ix = u.i & 0x7fffffff;
	if (ix > 0x42fc0000) {  /* |x| > 126 */
		if (u.i >= 0x43000000 && u.i < 0x80000000) {  /* x >= 128 */
			x *= 0x1p127f;
			return x;
		}
		if (u.i >= 0x80000000) {  /* x < -126 */
			if (u.i >= 0xc3160000 || (u.i & 0x0000ffff))
				FORCE_EVAL(-0x1p-149f/x);
			if (u.i >= 0xc3160000)  /* x <= -150 */
				return 0;
		}
	} else if (ix <= 0x33000000) {  /* |x| <= 0x1p-25 */
		return 1.0f + x;
	}

	/* Reduce x, computing z, i0, and k. */
	u.f = x + redux;
	i0 = u.i;
	i0 += TBLSIZE / 2;
	k = i0 / TBLSIZE;
	uk.i = (uint64_t)(0x3ff + k)<<52;
	i0 &= TBLSIZE - 1;
	u.f -= redux;
	z = x - u.f;
	/* Compute r = exp2(y) = exp2ft[i0] * p(z). */
	r = exp2ft[i0];
	t = r * z;
	r = r + t * (P1 + z * P2) + t * (z * z) * (P3 + z * P4);

	/* Scale by 2**k */
	return r * uk.f;
}


text:
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_lgammal.c */
/*
 * ====================================================
 * 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.
 * ====================================================
 */
/*
 * 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.
 */
/* lgammal(x)
 * 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)
 *   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)
 *      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
 *
 *   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(inf) = inf
 *              lgamma(-integer) = +-inf
 *
 */

#define _GNU_SOURCE
#include "libm.h"
#include "libc.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
double __lgamma_r(double x, int *sg);

long double __lgammal_r(long double x, int *sg)
{
	return __lgamma_r(x, sg);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
pi = 3.14159265358979323846264L,

/* lgam(1+x) = 0.5 x + x a(x)/b(x)
    -0.268402099609375 <= x <= 0
    peak relative error 6.6e-22 */
a0 = -6.343246574721079391729402781192128239938E2L,
a1 =  1.856560238672465796768677717168371401378E3L,
a2 =  2.404733102163746263689288466865843408429E3L,
a3 =  8.804188795790383497379532868917517596322E2L,
a4 =  1.135361354097447729740103745999661157426E2L,
a5 =  3.766956539107615557608581581190400021285E0L,

b0 =  8.214973713960928795704317259806842490498E3L,
b1 =  1.026343508841367384879065363925870888012E4L,
b2 =  4.553337477045763320522762343132210919277E3L,
b3 =  8.506975785032585797446253359230031874803E2L,
b4 =  6.042447899703295436820744186992189445813E1L,
/* b5 =  1.000000000000000000000000000000000000000E0 */

tc =  1.4616321449683623412626595423257213284682E0L,
tf = -1.2148629053584961146050602565082954242826E-1, /* double precision */
/* tt = (tail of tf), i.e. tf + tt has extended precision. */
tt = 3.3649914684731379602768989080467587736363E-18L,
/* lgam ( 1.4616321449683623412626595423257213284682E0 ) =
-1.2148629053584960809551455717769158215135617312999903886372437313313530E-1 */

/* lgam (x + tc) = tf + tt + x g(x)/h(x)
    -0.230003726999612341262659542325721328468 <= x
       <= 0.2699962730003876587373404576742786715318
     peak relative error 2.1e-21 */
g0 = 3.645529916721223331888305293534095553827E-18L,
g1 = 5.126654642791082497002594216163574795690E3L,
g2 = 8.828603575854624811911631336122070070327E3L,
g3 = 5.464186426932117031234820886525701595203E3L,
g4 = 1.455427403530884193180776558102868592293E3L,
g5 = 1.541735456969245924860307497029155838446E2L,
g6 = 4.335498275274822298341872707453445815118E0L,

h0 = 1.059584930106085509696730443974495979641E4L,
h1 = 2.147921653490043010629481226937850618860E4L,
h2 = 1.643014770044524804175197151958100656728E4L,
h3 = 5.869021995186925517228323497501767586078E3L,
h4 = 9.764244777714344488787381271643502742293E2L,
h5 = 6.442485441570592541741092969581997002349E1L,
/* h6 = 1.000000000000000000000000000000000000000E0 */

/* lgam (x+1) = -0.5 x + x u(x)/v(x)
    -0.100006103515625 <= x <= 0.231639862060546875
    peak relative error 1.3e-21 */
u0 = -8.886217500092090678492242071879342025627E1L,
u1 =  6.840109978129177639438792958320783599310E2L,
u2 =  2.042626104514127267855588786511809932433E3L,
u3 =  1.911723903442667422201651063009856064275E3L,
u4 =  7.447065275665887457628865263491667767695E2L,
u5 =  1.132256494121790736268471016493103952637E2L,
u6 =  4.484398885516614191003094714505960972894E0L,

v0 =  1.150830924194461522996462401210374632929E3L,
v1 =  3.399692260848747447377972081399737098610E3L,
v2 =  3.786631705644460255229513563657226008015E3L,
v3 =  1.966450123004478374557778781564114347876E3L,
v4 =  4.741359068914069299837355438370682773122E2L,
v5 =  4.508989649747184050907206782117647852364E1L,
/* v6 =  1.000000000000000000000000000000000000000E0 */

/* lgam (x+2) = .5 x + x s(x)/r(x)
     0 <= x <= 1
     peak relative error 7.2e-22 */
s0 =  1.454726263410661942989109455292824853344E6L,
s1 = -3.901428390086348447890408306153378922752E6L,
s2 = -6.573568698209374121847873064292963089438E6L,
s3 = -3.319055881485044417245964508099095984643E6L,
s4 = -7.094891568758439227560184618114707107977E5L,
s5 = -6.263426646464505837422314539808112478303E4L,
s6 = -1.684926520999477529949915657519454051529E3L,

r0 = -1.883978160734303518163008696712983134698E7L,
r1 = -2.815206082812062064902202753264922306830E7L,
r2 = -1.600245495251915899081846093343626358398E7L,
r3 = -4.310526301881305003489257052083370058799E6L,
r4 = -5.563807682263923279438235987186184968542E5L,
r5 = -3.027734654434169996032905158145259713083E4L,
r6 = -4.501995652861105629217250715790764371267E2L,
/* r6 =  1.000000000000000000000000000000000000000E0 */

/* lgam(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x w(1/x^2)
    x >= 8
    Peak relative error 1.51e-21
w0 = LS2PI - 0.5 */
w0 =  4.189385332046727417803e-1L,
w1 =  8.333333333333331447505E-2L,
w2 = -2.777777777750349603440E-3L,
w3 =  7.936507795855070755671E-4L,
w4 = -5.952345851765688514613E-4L,
w5 =  8.412723297322498080632E-4L,
w6 = -1.880801938119376907179E-3L,
w7 =  4.885026142432270781165E-3L;

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

	/* spurious inexact if odd int */
	x *= 0.5;
	x = 2.0*(x - floorl(x));  /* 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 __sinl(x, 0.0, 0);
	case 1: return __cosl(x, 0.0);
	case 2: return __sinl(-x, 0.0, 0);
	case 3: return -__cosl(x, 0.0);
	}
}

long double __lgammal_r(long double x, int *sg) {
	long double t, y, z, nadj, p, p1, p2, q, r, w;
	union ldshape u = {x};
	uint32_t ix = (u.i.se & 0x7fffU)<<16 | u.i.m>>48;
	int sign = u.i.se >> 15;
	int i;

	*sg = 1;

	/* purge off +-inf, NaN, +-0, tiny and negative arguments */
	if (ix >= 0x7fff0000)
		return x * x;
	if (ix < 0x3fc08000) {  /* |x|<2**-63, return -log(|x|) */
		if (sign) {
			*sg = -1;
			x = -x;
		}
		return -logl(x);
	}
	if (sign) {
		x = -x;
		t = sin_pi(x);
		if (t == 0.0)
			return 1.0 / (x-x); /* -integer */
		if (t > 0.0)
			*sg = -1;
		else
			t = -t;
		nadj = logl(pi / (t * x));
	}

	/* purge off 1 and 2 (so the sign is ok with downward rounding) */
	if ((ix == 0x3fff8000 || ix == 0x40008000) && u.i.m == 0) {
		r = 0;
	} else if (ix < 0x40008000) {  /* x < 2.0 */
		if (ix <= 0x3ffee666) {  /* 8.99993896484375e-1 */
			/* lgamma(x) = lgamma(x+1) - log(x) */
			r = -logl(x);
			if (ix >= 0x3ffebb4a) {  /* 7.31597900390625e-1 */
				y = x - 1.0;
				i = 0;
			} else if (ix >= 0x3ffced33) {  /* 2.31639862060546875e-1 */
				y = x - (tc - 1.0);
				i = 1;
			} else { /* x < 0.23 */
				y = x;
				i = 2;
			}
		} else {
			r = 0.0;
			if (ix >= 0x3fffdda6) {  /* 1.73162841796875 */
				/* [1.7316,2] */
				y = x - 2.0;
				i = 0;
			} else if (ix >= 0x3fff9da6) {  /* 1.23162841796875 */
				/* [1.23,1.73] */
				y = x - tc;
				i = 1;
			} else {
				/* [0.9, 1.23] */
				y = x - 1.0;
				i = 2;
			}
		}
		switch (i) {
		case 0:
			p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))));
			p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y))));
			r += 0.5 * y + y * p1/p2;
			break;
		case 1:
			p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6)))));
			p2 = h0 + y * (h1 + y * (h2 + y * (h3 + y * (h4 + y * (h5 + y)))));
			p = tt + y * p1/p2;
			r += (tf + p);
			break;
		case 2:
			p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6))))));
			p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y)))));
			r += (-0.5 * y + p1 / p2);
		}
	} else if (ix < 0x40028000) {  /* 8.0 */
		/* 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 = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
		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 += logl(z);
			break;
		}
	} else if (ix < 0x40418000) {  /* 2^66 */
		/* 8.0 <= x < 2**66 */
		t = logl(x);
		z = 1.0 / x;
		y = z * z;
		w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7))))));
		r = (x - 0.5) * (t - 1.0) + w;
	} else /* 2**66 <= x <= inf */
		r = x * (logl(x) - 1.0);
	if (sign)
		r = nadj - r;
	return r;
}
#endif

#if (LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024) || (LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384)
extern int __signgam;

long double lgammal(long double x)
{
	return __lgammal_r(x, &__signgam);
}

weak_alias(__lgammal_r, lgammal_r);
#endif


text:
#include <math.h>

long long llround(double x)
{
	return round(x);
}


text: sqrt.s


text:
.text
.global sqrtf
.type   sqrtf,%function
sqrtf:
	vsqrt.f32 s0, s0
	bx lr


text:
.text
.global fabsf
.type   fabsf,%function
fabsf:
	vabs.f32 s0, s0
	bx lr


text: fabs.s


text: sqrtf.s


text:
.text
.global fabs
.type   fabs,%function
fabs:
	vabs.f64 d0, d0
	bx lr


text:
.text
.global sqrt
.type   sqrt,%function
sqrt:
	vsqrt.f64 d0, d0
	bx lr


text: fabsf.s


text:
.global rintl
.type rintl,@function
rintl:
	fldt 4(%esp)
	frndint
	ret


text:
.global remainder
.type remainder,@function
remainder:
.weak drem
.type drem,@function
drem:
	fldl 12(%esp)
	fldl 4(%esp)
1:	fprem1
	fstsw %ax
	sahf
	jp 1b
	fstp %st(1)
	ret


text:
.global floorf
.type floorf,@function
floorf:
	flds 4(%esp)
	jmp 1f

.global floorl
.type floorl,@function
floorl:
	fldt 4(%esp)
	jmp 1f

.global floor
.type floor,@function
floor:
	fldl 4(%esp)
1:	mov $0x7,%al
1:	fstcw 4(%esp)
	mov 5(%esp),%ah
	mov %al,5(%esp)
	fldcw 4(%esp)
	frndint
	mov %ah,5(%esp)
	fldcw 4(%esp)
	ret

.global ceil
.type ceil,@function
ceil:
	fldl 4(%esp)
	mov $0xb,%al
	jmp 1b

.global ceilf
.type ceilf,@function
ceilf:
	flds 4(%esp)
	mov $0xb,%al
	jmp 1b

.global ceill
.type ceill,@function
ceill:
	fldt 4(%esp)
	mov $0xb,%al
	jmp 1b

.global trunc
.type trunc,@function
trunc:
	fldl 4(%esp)
	mov $0xf,%al
	jmp 1b

.global truncf
.type truncf,@function
truncf:
	flds 4(%esp)
	mov $0xf,%al
	jmp 1b

.global truncl
.type truncl,@function
truncl:
	fldt 4(%esp)
	mov $0xf,%al
	jmp 1b


text: # see floor.s


text: # see exp.s


text: # see scalbnf.s


text:
.global log10f
.type log10f,@function
log10f:
	fldlg2
	flds 4(%esp)
	fyl2x
	ret


text: # see acos.s


text:
.global log
.type log,@function
log:
	fldln2
	fldl 4(%esp)
	fyl2x
	ret


text:
.global sqrtf
.type sqrtf,@function
sqrtf:	flds 4(%esp)
	fsqrt
	fstps 4(%esp)
	flds 4(%esp)
	ret


text:
.global remquof
.type remquof,@function
remquof:
	mov 12(%esp),%ecx
	flds 8(%esp)
	flds 4(%esp)
	mov 11(%esp),%dh
	xor 7(%esp),%dh
	jmp 1f

.global remquol
.type remquol,@function
remquol:
	mov 28(%esp),%ecx
	fldt 16(%esp)
	fldt 4(%esp)
	mov 25(%esp),%dh
	xor 13(%esp),%dh
	jmp 1f

.global remquo
.type remquo,@function
remquo:
	mov 20(%esp),%ecx
	fldl 12(%esp)
	fldl 4(%esp)
	mov 19(%esp),%dh
	xor 11(%esp),%dh
1:	fprem1
	fnstsw %ax
	sahf
	jp 1b
	fstp %st(1)
	mov %ah,%dl
	shr %dl
	and $1,%dl
	mov %ah,%al
	shr $5,%al
	and $2,%al
	or %al,%dl
	mov %ah,%al
	shl $2,%al
	and $4,%al
	or %al,%dl
	test %dh,%dh
	jns 1f
	neg %dl
1:	movsbl %dl,%edx
	mov %edx,(%ecx)
	ret


text:
.global log1pl
.type log1pl,@function
log1pl:
	mov 10(%esp),%eax
	fldln2
	and $0x7fffffff,%eax
	fldt 4(%esp)
	cmp $0x3ffd9400,%eax
	ja 1f
	fyl2xp1
	ret
1:	fld1
	faddp
	fyl2x
	ret


text:
.global fmod
.type fmod,@function
fmod:
	fldl 12(%esp)
	fldl 4(%esp)
1:	fprem
	fstsw %ax
	sahf
	jp 1b
	fstp %st(1)
	ret


text:
.global fmodf
.type fmodf,@function
fmodf:
	flds 8(%esp)
	flds 4(%esp)
1:	fprem
	fstsw %ax
	sahf
	jp 1b
	fstp %st(1)
	ret


text:
.global sqrtl
.type sqrtl,@function
sqrtl:	fldt 4(%esp)
	fsqrt
	ret


text:
.global lrintl
.type lrintl,@function
lrintl:
	fldt 4(%esp)
	fistpl 4(%esp)
	mov 4(%esp),%eax
	ret


text:
.global logf
.type logf,@function
logf:
	fldln2
	flds 4(%esp)
	fyl2x
	ret


text:
.global ldexpl
.type ldexpl,@function
ldexpl:
	nop

.global scalblnl
.type scalblnl,@function
scalblnl:
	nop

.global scalbnl
.type scalbnl,@function
scalbnl:
	mov 16(%esp),%eax
	add $0x3ffe,%eax
	cmp $0x7ffd,%eax
	jae 1f
	inc %eax
	fldt 4(%esp)
	mov %eax,12(%esp)
	mov $0x80000000,%eax
	mov %eax,8(%esp)
	xor %eax,%eax
	mov %eax,4(%esp)
	fldt 4(%esp)
	fmulp
	ret
1:	fildl 16(%esp)
	fldt 4(%esp)
	fscale
	fstp %st(1)
	ret


text:
.global rint
.type rint,@function
rint:
	fldl 4(%esp)
	frndint
	ret


text:
.global fabsf
.type fabsf,@function
fabsf:
	flds 4(%esp)
	fabs
	ret


text:
.global log10l
.type log10l,@function
log10l:
	fldlg2
	fldt 4(%esp)
	fyl2x
	ret


text:
.global logl
.type logl,@function
logl:
	fldln2
	fldt 4(%esp)
	fyl2x
	ret


text:
.global log2
.type log2,@function
log2:
	fld1
	fldl 4(%esp)
	fyl2x
	ret


text: # see floor.s


text:
.global remainderf
.type remainderf,@function
remainderf:
.weak dremf
.type dremf,@function
dremf:
	flds 8(%esp)
	flds 4(%esp)
1:	fprem1
	fstsw %ax
	sahf
	jp 1b
	fstp %st(1)
	ret


text:
.global atan
.type atan,@function
atan:
	fldl 4(%esp)
	mov 8(%esp),%eax
	add %eax,%eax
	cmp $0x00200000,%eax
	jb 1f
	fld1
	fpatan
	ret
		# subnormal x, return x with underflow
1:	fnstsw %ax
	and $16,%ax
	jnz 2f
	fsts 4(%esp)
2:	ret


text:
.global lrintf
.type lrintf,@function
lrintf:
	flds 4(%esp)
	fistpl 4(%esp)
	mov 4(%esp),%eax
	ret


text:
.global log10
.type log10,@function
log10:
	fldlg2
	fldl 4(%esp)
	fyl2x
	ret


text:
.global llrint
.type llrint,@function
llrint:
	fldl 4(%esp)
	fistpll 4(%esp)
	mov 4(%esp),%eax
	mov 8(%esp),%edx
	ret


text: # see exp.s


text:
.global log2l
.type log2l,@function
log2l:
	fld1
	fldt 4(%esp)
	fyl2x
	ret


text:
.global hypotf
.type hypotf,@function
hypotf:
	mov 4(%esp),%eax
	mov 8(%esp),%ecx
	add %eax,%eax
	add %ecx,%ecx
	and %eax,%ecx
	cmp $0xff000000,%ecx
	jae 2f
	test %eax,%eax
	jnz 1f
	flds 8(%esp)
	fabs
	ret
1:	mov 8(%esp),%eax
	add %eax,%eax
	jnz 1f
	flds 4(%esp)
	fabs
	ret
1:	flds 4(%esp)
	fld %st(0)
	fmulp
	flds 8(%esp)
	fld %st(0)
	fmulp
	faddp
	fsqrt
	ret
2:	cmp $0xff000000,%eax
	jnz 1f
	flds 4(%esp)
	fabs
	ret
1:	mov 8(%esp),%eax
	add %eax,%eax
	cmp $0xff000000,%eax
	flds 8(%esp)
	jnz 1f
	fabs
1:	ret


text: # see asin.s


text: # see scalbnf.s


text: # see scalbnl.s


text:
.global fabs
.type fabs,@function
fabs:
	fldl 4(%esp)
	fabs
	ret


text:
# exp(x) = 2^hi + 2^hi (2^lo - 1)
# where hi+lo = log2e*x with 128bit precision
# exact log2e*x calculation depends on nearest rounding mode
# using the exact multiplication method of Dekker and Veltkamp

.global expl
.type expl,@function
expl:
	fldt 4(%esp)

		# interesting case: 0x1p-32 <= |x| < 16384
		# check if (exponent|0x8000) is in [0xbfff-32, 0xbfff+13]
	mov 12(%esp), %ax
	or $0x8000, %ax
	sub $0xbfdf, %ax
	cmp $45, %ax
	jbe 2f
	test %ax, %ax
	fld1
	js 1f
		# if |x|>=0x1p14 or nan return 2^trunc(x)
	fscale
	fstp %st(1)
	ret
		# if |x|<0x1p-32 return 1+x
1:	faddp
	ret

		# should be 0x1.71547652b82fe178p0L == 0x3fff b8aa3b29 5c17f0bc
		# it will be wrong on non-nearest rounding mode
2:	fldl2e
	subl $44, %esp
		# hi = log2e_hi*x
		# 2^hi = exp2l(hi)
	fmul %st(1),%st
	fld %st(0)
	fstpt (%esp)
	fstpt 16(%esp)
	fstpt 32(%esp)
	call exp2l
		# if 2^hi == inf return 2^hi
	fld %st(0)
	fstpt (%esp)
	cmpw $0x7fff, 8(%esp)
	je 1f
	fldt 32(%esp)
	fldt 16(%esp)
		# fpu stack: 2^hi x hi
		# exact mult: x*log2e
	fld %st(1)
		# c = 0x1p32+1
	pushl $0x41f00000
	pushl $0x00100000
	fldl (%esp)
		# xh = x - c*x + c*x
		# xl = x - xh
	fmulp
	fld %st(2)
	fsub %st(1), %st
	faddp
	fld %st(2)
	fsub %st(1), %st
		# yh = log2e_hi - c*log2e_hi + c*log2e_hi
	pushl $0x3ff71547
	pushl $0x65200000
	fldl (%esp)
		# fpu stack: 2^hi x hi xh xl yh
		# lo = hi - xh*yh + xl*yh
	fld %st(2)
	fmul %st(1), %st
	fsubp %st, %st(4)
	fmul %st(1), %st
	faddp %st, %st(3)
		# yl = log2e_hi - yh
	pushl $0x3de705fc
	pushl $0x2f000000
	fldl (%esp)
		# fpu stack: 2^hi x lo xh xl yl
		# lo += xh*yl + xl*yl
	fmul %st, %st(2)
	fmulp %st, %st(1)
	fxch %st(2)
	faddp
	faddp
		# log2e_lo
	pushl $0xbfbe
	pushl $0x82f0025f
	pushl $0x2dc582ee
	fldt (%esp)
	addl $36,%esp
		# fpu stack: 2^hi x lo log2e_lo
		# lo += log2e_lo*x
		# return 2^hi + 2^hi (2^lo - 1)
	fmulp %st, %st(2)
	faddp
	f2xm1
	fmul %st(1), %st
	faddp
1:	addl $44, %esp
	ret


text:
.global lrint
.type lrint,@function
lrint:
	fldl 4(%esp)
	fistpl 4(%esp)
	mov 4(%esp),%eax
	ret


text: # see remquo.s


text:
.global atanf
.type atanf,@function
atanf:
	flds 4(%esp)
	mov 4(%esp),%eax
	add %eax,%eax
	cmp $0x01000000,%eax
	jb 1f
	fld1
	fpatan
	ret
		# subnormal x, return x with underflow
1:	fnstsw %ax
	and $16,%ax
	jnz 2f
	fld %st(0)
	fmul %st(1)
	fstps 4(%esp)
2:	ret


text:
.global rintf
.type rintf,@function
rintf:
	flds 4(%esp)
	frndint
	ret


text:
.global hypot
.type hypot,@function
hypot:
	mov 8(%esp),%eax
	mov 16(%esp),%ecx
	add %eax,%eax
	add %ecx,%ecx
	and %eax,%ecx
	cmp $0xffe00000,%ecx
	jae 2f
	or 4(%esp),%eax
	jnz 1f
	fldl 12(%esp)
	fabs
	ret
1:	mov 16(%esp),%eax
	add %eax,%eax
	or 12(%esp),%eax
	jnz 1f
	fldl 4(%esp)
	fabs
	ret
1:	fldl 4(%esp)
	fld %st(0)
	fmulp
	fldl 12(%esp)
	fld %st(0)
	fmulp
	faddp
	fsqrt
	ret
2:	sub $0xffe00000,%eax
	or 4(%esp),%eax
	jnz 1f
	fldl 4(%esp)
	fabs
	ret
1:	mov 16(%esp),%eax
	add %eax,%eax
	sub $0xffe00000,%eax
	or 12(%esp),%eax
	fldl 12(%esp)
	jnz 1f
	fabs
1:	ret


text:
.global log1pf
.type log1pf,@function
log1pf:
	mov 4(%esp),%eax
	fldln2
	and $0x7fffffff,%eax
	flds 4(%esp)
	cmp $0x3e940000,%eax
	ja 1f
	cmp $0x00800000,%eax
	jb 2f
	fyl2xp1
	ret
1:	fld1
	faddp
	fyl2x
	ret
		# subnormal x, return x with underflow
2:	fnstsw %ax
	and $16,%ax
	jnz 1f
	fxch
	fmul %st(1)
	fstps 4(%esp)
1:	ret


text:
.global log1p
.type log1p,@function
log1p:
	mov 8(%esp),%eax
	fldln2
	and $0x7fffffff,%eax
	fldl 4(%esp)
	cmp $0x3fd28f00,%eax
	ja 1f
	cmp $0x00100000,%eax
	jb 2f
	fyl2xp1
	ret
1:	fld1
	faddp
	fyl2x
	ret
		# subnormal x, return x with underflow
2:	fnstsw %ax
	and $16,%ax
	jnz 1f
	fsts 4(%esp)
	fstp %st(1)
1:	ret


text:
.global ldexp
.type ldexp,@function
ldexp:
	nop

.global scalbln
.type scalbln,@function
scalbln:
	nop

.global scalbn
.type scalbn,@function
scalbn:
	mov 12(%esp),%eax
	add $0x3ffe,%eax
	cmp $0x7ffd,%eax
	jb 1f
	sub $0x3ffe,%eax
	sar $31,%eax
	xor $0xfff,%eax
	add $0x3ffe,%eax
1:	inc %eax
	fldl 4(%esp)
	mov %eax,12(%esp)
	mov $0x80000000,%eax
	mov %eax,8(%esp)
	xor %eax,%eax
	mov %eax,4(%esp)
	fldt 4(%esp)
	fmulp
	fstpl 4(%esp)
	fldl 4(%esp)
	ret


text:
.global asinf
.type asinf,@function
asinf:
	flds 4(%esp)
	mov 4(%esp),%eax
	add %eax,%eax
	cmp $0x01000000,%eax
	jae 1f
		# subnormal x, return x with underflow
	fnstsw %ax
	and $16,%ax
	jnz 2f
	fld %st(0)
	fmul %st(1)
	fstps 4(%esp)
2:	ret

.global asinl
.type asinl,@function
asinl:
	fldt 4(%esp)
	jmp 1f

.global asin
.type asin,@function
asin:
	fldl 4(%esp)
	mov 8(%esp),%eax
	add %eax,%eax
	cmp $0x00200000,%eax
	jae 1f
		# subnormal x, return x with underflow
	fnstsw %ax
	and $16,%ax
	jnz 2f
	fsts 4(%esp)
2:	ret
1:	fld %st(0)
	fld1
	fsub %st(0),%st(1)
	fadd %st(2)
	fmulp
	fsqrt
	fpatan
	ret


text: # see scalbn.s


text:
.global fmodl
.type fmodl,@function
fmodl:
	fldt 16(%esp)
	fldt 4(%esp)
1:	fprem
	fstsw %ax
	sahf
	jp 1b
	fstp %st(1)
	ret


text:
.global atanl
.type atanl,@function
atanl:
	fldt 4(%esp)
	fld1
	fpatan
	ret


text:
.global log2f
.type log2f,@function
log2f:
	fld1
	flds 4(%esp)
	fyl2x
	ret


text:
.global atan2f
.type atan2f,@function
atan2f:
	flds 4(%esp)
	flds 8(%esp)
	fpatan
	fsts 4(%esp)
	mov 4(%esp),%eax
	add %eax,%eax
	cmp $0x01000000,%eax
	jae 1f
		# subnormal x, return x with underflow
	fnstsw %ax
	and $16,%ax
	jnz 1f
	fld %st(0)
	fmul %st(1)
	fstps 4(%esp)
1:	ret


text:
.global sqrt
.type sqrt,@function
sqrt:	fldl 4(%esp)
	fsqrt
	fstsw %ax
	sub $12,%esp
	fld %st(0)
	fstpt (%esp)
	mov (%esp),%ecx
	and $0x7ff,%ecx
	cmp $0x400,%ecx
	jnz 1f
	and $0x200,%eax
	sub $0x100,%eax
	sub %eax,(%esp)
	fstp %st(0)
	fldt (%esp)
1:	add $12,%esp
	fstpl 4(%esp)
	fldl 4(%esp)
	ret


text:
.global llrintf
.type llrintf,@function
llrintf:
	sub $8,%esp
	flds 12(%esp)
	fistpll (%esp)
	pop %eax
	pop %edx
	ret


text:
.global expm1f
.type expm1f,@function
expm1f:
	flds 4(%esp)
	mov 4(%esp),%eax
	add %eax,%eax
	cmp $0x01000000,%eax
	jae 1f
		# subnormal x, return x with underflow
	fnstsw %ax
	and $16,%ax
	jnz 2f
	fld %st(0)
	fmul %st(1)
	fstps 4(%esp)
2:	ret

.global expm1l
.type expm1l,@function
expm1l:
	fldt 4(%esp)
	jmp 1f

.global expm1
.type expm1,@function
expm1:
	fldl 4(%esp)
	mov 8(%esp),%eax
	add %eax,%eax
	cmp $0x00200000,%eax
	jae 1f
		# subnormal x, return x with underflow
	fnstsw %ax
	and $16,%ax
	jnz 2f
	fsts 4(%esp)
2:	ret
1:	fldl2e
	fmulp
	mov $0xc2820000,%eax
	push %eax
	flds (%esp)
	pop %eax
	fucomp %st(1)
	fnstsw %ax
	sahf
	fld1
	jb 1f
		# x*log2e < -65, return -1 without underflow
	fstp %st(1)
	fchs
	ret
1:	fld %st(1)
	fabs
	fucom %st(1)
	fnstsw %ax
	fstp %st(0)
	fstp %st(0)
	sahf
	ja 1f
	f2xm1
	ret
1:	call 1f
	fld1
	fsubrp
	ret

.global exp2f
.type exp2f,@function
exp2f:
	flds 4(%esp)
	jmp 1f

.global exp2l
.type exp2l,@function
exp2l:
	fldt 4(%esp)
	jmp 1f

.global expf
.type expf,@function
expf:
	flds 4(%esp)
	jmp 2f

.global exp
.type exp,@function
exp:
	fldl 4(%esp)
2:	fldl2e
	fmulp
	jmp 1f

.global exp2
.type exp2,@function
exp2:
	fldl 4(%esp)
1:	sub $12,%esp
	fld %st(0)
	fstpt (%esp)
	mov 8(%esp),%ax
	and $0x7fff,%ax
	cmp $0x3fff+13,%ax
	jb 4f             # |x| < 8192
	cmp $0x3fff+15,%ax
	jae 3f            # |x| >= 32768
	fsts (%esp)
	cmpl $0xc67ff800,(%esp)
	jb 2f             # x > -16382
	movl $0x5f000000,(%esp)
	flds (%esp)       # 0x1p63
	fld %st(1)
	fsub %st(1)
	faddp
	fucomp %st(1)
	fnstsw
	sahf
	je 2f             # x - 0x1p63 + 0x1p63 == x
	movl $1,(%esp)
	flds (%esp)       # 0x1p-149
	fdiv %st(1)
	fstps (%esp)      # raise underflow
2:	fld1
	fld %st(1)
	frndint
	fxch %st(2)
	fsub %st(2)       # st(0)=x-rint(x), st(1)=1, st(2)=rint(x)
	f2xm1
	faddp             # 2^(x-rint(x))
1:	fscale
	fstp %st(1)
	add $12,%esp
	ret
3:	xor %eax,%eax
4:	cmp $0x3fff-64,%ax
	fld1
	jb 1b             # |x| < 0x1p-64
	fstpt (%esp)
	fistl 8(%esp)
	fildl 8(%esp)
	fsubrp %st(1)
	addl $0x3fff,8(%esp)
	f2xm1
	fld1
	faddp             # 2^(x-rint(x))
	fldt (%esp)       # 2^rint(x)
	fmulp
	add $12,%esp
	ret


text:
.global atan2
.type atan2,@function
atan2:
	fldl 4(%esp)
	fldl 12(%esp)
	fpatan
	fstl 4(%esp)
	mov 8(%esp),%eax
	add %eax,%eax
	cmp $0x00200000,%eax
	jae 1f
		# subnormal x, return x with underflow
	fnstsw %ax
	and $16,%ax
	jnz 1f
	fsts 4(%esp)
1:	ret


text:
.global atan2l
.type atan2l,@function
atan2l:
	fldt 4(%esp)
	fldt 16(%esp)
	fpatan
	ret


text: # see acos.s


text:
.global ldexpf
.type ldexpf,@function
ldexpf:
	nop

.global scalblnf
.type scalblnf,@function
scalblnf:
	nop

.global scalbnf
.type scalbnf,@function
scalbnf:
	mov 8(%esp),%eax
	add $0x3fe,%eax
	cmp $0x7fd,%eax
	jb 1f
	sub $0x3fe,%eax
	sar $31,%eax
	xor $0x1ff,%eax
	add $0x3fe,%eax
1:	inc %eax
	shl $20,%eax
	flds 4(%esp)
	mov %eax,8(%esp)
	xor %eax,%eax
	mov %eax,4(%esp)
	fldl 4(%esp)
	fmulp
	fstps 4(%esp)
	flds 4(%esp)
	ret


text:
.global fabsl
.type fabsl,@function
fabsl:
	fldt 4(%esp)
	fabs
	ret


text: # see exp.s


text:
.global llrintl
.type llrintl,@function
llrintl:
	fldt 4(%esp)
	fistpll 4(%esp)
	mov 4(%esp),%eax
	mov 8(%esp),%edx
	ret


text: # see floor.s


text:
# use acos(x) = atan2(fabs(sqrt((1-x)*(1+x))), x)

.global acosf
.type acosf,@function
acosf:
	flds 4(%esp)
	jmp 1f

.global acosl
.type acosl,@function
acosl:
	fldt 4(%esp)
	jmp 1f

.global acos
.type acos,@function
acos:
	fldl 4(%esp)
1:	fld %st(0)
	fld1
	fsub %st(0),%st(1)
	fadd %st(2)
	fmulp
	fsqrt
	fabs         # fix sign of zero (matters in downward rounding mode)
	fxch %st(1)
	fpatan
	ret


text:
.global remainderl
.type remainderl,@function
remainderl:
	fldt 16(%esp)
	fldt 4(%esp)
1:	fprem1
	fstsw %ax
	sahf
	jp 1b
	fstp %st(1)
	ret


text: # see asin.s


text:
.global rintl
.type rintl,@function
rintl:
	fldt 8(%esp)
	frndint
	ret


text:
# see ../i386/acos.s

.global acosl
.type acosl,@function
acosl:
	fldt 8(%esp)
1:	fld %st(0)
	fld1
	fsub %st(0),%st(1)
	fadd %st(2)
	fmulp
	fsqrt
	fabs
	fxch %st(1)
	fpatan
	ret


text:
.global sqrtf
.type sqrtf,@function
sqrtf:  sqrtss %xmm0, %xmm0
	ret


text:
.global log1pl
.type log1pl,@function
log1pl:
	mov 14(%esp),%eax
	fldln2
	and $0x7fffffff,%eax
	fldt 8(%esp)
	cmp $0x3ffd9400,%eax
	ja 1f
	fyl2xp1
	ret
1:	fld1
	faddp
	fyl2x
	ret


text:
.global sqrtl
.type sqrtl,@function
sqrtl:	fldt 8(%esp)
	fsqrt
	ret


text:
.global lrintl
.type lrintl,@function
lrintl:
	fldt 8(%esp)
	fistpll 8(%esp)
	mov 8(%esp),%rax
	ret