Add gamma functions.

12 files changed, 1563 insertions, 0 deletions

diff --git a/system/lib/libc/musl/src/math/__cosdf.c b/system/lib/libc/musl/src/math/__cosdf.c
new file mode 100644
index 00000000..2124989b
--- /dev/null
+++ b/system/lib/libc/musl/src/math/__cosdf.c
@@ -0,0 +1,35 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/k_cosf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Debugged and optimized by Bruce D. Evans.
+ */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include "libm.h"
+
+/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
+static const double
+C0  = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */
+C1  =  0x155553e1053a42.0p-57, /*  0.0416666233237390631894 */
+C2  = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */
+C3  =  0x199342e0ee5069.0p-68; /*  0.0000243904487962774090654 */
+
+float __cosdf(double x)
+{
+	double_t r, w, z;
+
+	/* Try to optimize for parallel evaluation as in __tandf.c. */
+	z = x*x;
+	w = z*z;
+	r = C2+z*C3;
+	return ((1.0+z*C0) + w*C1) + (w*z)*r;
+}
diff --git a/system/lib/libc/musl/src/math/__sindf.c b/system/lib/libc/musl/src/math/__sindf.c
new file mode 100644
index 00000000..8fec2a3f
--- /dev/null
+++ b/system/lib/libc/musl/src/math/__sindf.c
@@ -0,0 +1,36 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/k_sinf.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"
+
+/* |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). */
+static const double
+S1 = -0x15555554cbac77.0p-55, /* -0.166666666416265235595 */
+S2 =  0x111110896efbb2.0p-59, /*  0.0083333293858894631756 */
+S3 = -0x1a00f9e2cae774.0p-65, /* -0.000198393348360966317347 */
+S4 =  0x16cd878c3b46a7.0p-71; /*  0.0000027183114939898219064 */
+
+float __sindf(double x)
+{
+	double_t r, s, w, z;
+
+	/* Try to optimize for parallel evaluation as in __tandf.c. */
+	z = x*x;
+	w = z*z;
+	r = S3 + z*S4;
+	s = z*x;
+	return (x + s*(S1 + z*S2)) + s*w*r;
+}
diff --git a/system/lib/libc/musl/src/math/lgamma.c b/system/lib/libc/musl/src/math/lgamma.c
new file mode 100644
index 00000000..e25ec8e6
--- /dev/null
+++ b/system/lib/libc/musl/src/math/lgamma.c
@@ -0,0 +1,9 @@
+#include <math.h>
+
+extern int __signgam;
+double __lgamma_r(double, int *);
+
+double lgamma(double x)
+{
+	return __lgamma_r(x, &__signgam);
+}
diff --git a/system/lib/libc/musl/src/math/lgamma_r.c b/system/lib/libc/musl/src/math/lgamma_r.c
new file mode 100644
index 00000000..82e296f5
--- /dev/null
+++ b/system/lib/libc/musl/src/math/lgamma_r.c
@@ -0,0 +1,314 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_lgamma_r.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ */
+/* lgamma_r(x, signgamp)
+ * Reentrant version of the logarithm of the Gamma function
+ * with user provide pointer for the sign of Gamma(x).
+ *
+ * Method:
+ *   1. Argument Reduction for 0 < x <= 8
+ *      Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
+ *      reduce x to a number in [1.5,2.5] by
+ *              lgamma(1+s) = log(s) + lgamma(s)
+ *      for example,
+ *              lgamma(7.3) = log(6.3) + lgamma(6.3)
+ *                          = log(6.3*5.3) + lgamma(5.3)
+ *                          = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
+ *   2. Polynomial approximation of lgamma around its
+ *      minimun ymin=1.461632144968362245 to maintain monotonicity.
+ *      On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
+ *              Let z = x-ymin;
+ *              lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
+ *      where
+ *              poly(z) is a 14 degree polynomial.
+ *   2. Rational approximation in the primary interval [2,3]
+ *      We use the following approximation:
+ *              s = x-2.0;
+ *              lgamma(x) = 0.5*s + s*P(s)/Q(s)
+ *      with accuracy
+ *              |P/Q - (lgamma(x)-0.5s)| < 2**-61.71
+ *      Our algorithms are based on the following observation
+ *
+ *                             zeta(2)-1    2    zeta(3)-1    3
+ * lgamma(2+s) = s*(1-Euler) + --------- * s  -  --------- * s  + ...
+ *                                 2                 3
+ *
+ *      where Euler = 0.5771... is the Euler constant, which is very
+ *      close to 0.5.
+ *
+ *   3. For x>=8, we have
+ *      lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
+ *      (better formula:
+ *         lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
+ *      Let z = 1/x, then we approximation
+ *              f(z) = lgamma(x) - (x-0.5)(log(x)-1)
+ *      by
+ *                                  3       5             11
+ *              w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
+ *      where
+ *              |w - f(z)| < 2**-58.74
+ *
+ *   4. For negative x, since (G is gamma function)
+ *              -x*G(-x)*G(x) = pi/sin(pi*x),
+ *      we have
+ *              G(x) = pi/(sin(pi*x)*(-x)*G(-x))
+ *      since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
+ *      Hence, for x<0, signgam = sign(sin(pi*x)) and
+ *              lgamma(x) = log(|Gamma(x)|)
+ *                        = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
+ *      Note: one should avoid compute pi*(-x) directly in the
+ *            computation of sin(pi*(-x)).
+ *
+ *   5. Special Cases
+ *              lgamma(2+s) ~ s*(1-Euler) for tiny s
+ *              lgamma(1) = lgamma(2) = 0
+ *              lgamma(x) ~ -log(|x|) for tiny x
+ *              lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero
+ *              lgamma(inf) = inf
+ *              lgamma(-inf) = inf (bug for bug compatible with C99!?)
+ *
+ */
+
+#include "libm.h"
+#include "libc.h"
+
+static const double
+two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
+pi  =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
+a0  =  7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
+a1  =  3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
+a2  =  6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */
+a3  =  2.05808084325167332806e-02, /* 0x3F951322, 0xAC92547B */
+a4  =  7.38555086081402883957e-03, /* 0x3F7E404F, 0xB68FEFE8 */
+a5  =  2.89051383673415629091e-03, /* 0x3F67ADD8, 0xCCB7926B */
+a6  =  1.19270763183362067845e-03, /* 0x3F538A94, 0x116F3F5D */
+a7  =  5.10069792153511336608e-04, /* 0x3F40B6C6, 0x89B99C00 */
+a8  =  2.20862790713908385557e-04, /* 0x3F2CF2EC, 0xED10E54D */
+a9  =  1.08011567247583939954e-04, /* 0x3F1C5088, 0x987DFB07 */
+a10 =  2.52144565451257326939e-05, /* 0x3EFA7074, 0x428CFA52 */
+a11 =  4.48640949618915160150e-05, /* 0x3F07858E, 0x90A45837 */
+tc  =  1.46163214496836224576e+00, /* 0x3FF762D8, 0x6356BE3F */
+tf  = -1.21486290535849611461e-01, /* 0xBFBF19B9, 0xBCC38A42 */
+/* tt = -(tail of tf) */
+tt  = -3.63867699703950536541e-18, /* 0xBC50C7CA, 0xA48A971F */
+t0  =  4.83836122723810047042e-01, /* 0x3FDEF72B, 0xC8EE38A2 */
+t1  = -1.47587722994593911752e-01, /* 0xBFC2E427, 0x8DC6C509 */
+t2  =  6.46249402391333854778e-02, /* 0x3FB08B42, 0x94D5419B */
+t3  = -3.27885410759859649565e-02, /* 0xBFA0C9A8, 0xDF35B713 */
+t4  =  1.79706750811820387126e-02, /* 0x3F9266E7, 0x970AF9EC */
+t5  = -1.03142241298341437450e-02, /* 0xBF851F9F, 0xBA91EC6A */
+t6  =  6.10053870246291332635e-03, /* 0x3F78FCE0, 0xE370E344 */
+t7  = -3.68452016781138256760e-03, /* 0xBF6E2EFF, 0xB3E914D7 */
+t8  =  2.25964780900612472250e-03, /* 0x3F6282D3, 0x2E15C915 */
+t9  = -1.40346469989232843813e-03, /* 0xBF56FE8E, 0xBF2D1AF1 */
+t10 =  8.81081882437654011382e-04, /* 0x3F4CDF0C, 0xEF61A8E9 */
+t11 = -5.38595305356740546715e-04, /* 0xBF41A610, 0x9C73E0EC */
+t12 =  3.15632070903625950361e-04, /* 0x3F34AF6D, 0x6C0EBBF7 */
+t13 = -3.12754168375120860518e-04, /* 0xBF347F24, 0xECC38C38 */
+t14 =  3.35529192635519073543e-04, /* 0x3F35FD3E, 0xE8C2D3F4 */
+u0  = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
+u1  =  6.32827064025093366517e-01, /* 0x3FE4401E, 0x8B005DFF */
+u2  =  1.45492250137234768737e+00, /* 0x3FF7475C, 0xD119BD6F */
+u3  =  9.77717527963372745603e-01, /* 0x3FEF4976, 0x44EA8450 */
+u4  =  2.28963728064692451092e-01, /* 0x3FCD4EAE, 0xF6010924 */
+u5  =  1.33810918536787660377e-02, /* 0x3F8B678B, 0xBF2BAB09 */
+v1  =  2.45597793713041134822e+00, /* 0x4003A5D7, 0xC2BD619C */
+v2  =  2.12848976379893395361e+00, /* 0x40010725, 0xA42B18F5 */
+v3  =  7.69285150456672783825e-01, /* 0x3FE89DFB, 0xE45050AF */
+v4  =  1.04222645593369134254e-01, /* 0x3FBAAE55, 0xD6537C88 */
+v5  =  3.21709242282423911810e-03, /* 0x3F6A5ABB, 0x57D0CF61 */
+s0  = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
+s1  =  2.14982415960608852501e-01, /* 0x3FCB848B, 0x36E20878 */
+s2  =  3.25778796408930981787e-01, /* 0x3FD4D98F, 0x4F139F59 */
+s3  =  1.46350472652464452805e-01, /* 0x3FC2BB9C, 0xBEE5F2F7 */
+s4  =  2.66422703033638609560e-02, /* 0x3F9B481C, 0x7E939961 */
+s5  =  1.84028451407337715652e-03, /* 0x3F5E26B6, 0x7368F239 */
+s6  =  3.19475326584100867617e-05, /* 0x3F00BFEC, 0xDD17E945 */
+r1  =  1.39200533467621045958e+00, /* 0x3FF645A7, 0x62C4AB74 */
+r2  =  7.21935547567138069525e-01, /* 0x3FE71A18, 0x93D3DCDC */
+r3  =  1.71933865632803078993e-01, /* 0x3FC601ED, 0xCCFBDF27 */
+r4  =  1.86459191715652901344e-02, /* 0x3F9317EA, 0x742ED475 */
+r5  =  7.77942496381893596434e-04, /* 0x3F497DDA, 0xCA41A95B */
+r6  =  7.32668430744625636189e-06, /* 0x3EDEBAF7, 0xA5B38140 */
+w0  =  4.18938533204672725052e-01, /* 0x3FDACFE3, 0x90C97D69 */
+w1  =  8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */
+w2  = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */
+w3  =  7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */
+w4  = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
+w5  =  8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
+w6  = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
+
+static double sin_pi(double x)
+{
+	double y,z;
+	int n,ix;
+
+	GET_HIGH_WORD(ix, x);
+	ix &= 0x7fffffff;
+
+	if (ix < 0x3fd00000)
+		return __sin(pi*x, 0.0, 0);
+
+	y = -x;  /* negative x is assumed */
+
+	/*
+	 * argument reduction, make sure inexact flag not raised if input
+	 * is an integer
+	 */
+	z = floor(y);
+	if (z != y) {    /* inexact anyway */
+		y *= 0.5;
+		y  = 2.0*(y - floor(y));   /* y = |x| mod 2.0 */
+		n  = (int)(y*4.0);
+	} else {
+		if (ix >= 0x43400000) {
+			y = 0.0;    /* y must be even */
+			n = 0;
+		} else {
+			if (ix < 0x43300000)
+				z = y + two52;  /* exact */
+			GET_LOW_WORD(n, z);
+			n &= 1;
+			y = n;
+			n <<= 2;
+		}
+	}
+	switch (n) {
+	case 0:  y =  __sin(pi*y, 0.0, 0); break;
+	case 1:
+	case 2:  y =  __cos(pi*(0.5-y), 0.0); break;
+	case 3:
+	case 4:  y =  __sin(pi*(1.0-y), 0.0, 0); break;
+	case 5:
+	case 6:  y = -__cos(pi*(y-1.5), 0.0); break;
+	default: y =  __sin(pi*(y-2.0), 0.0, 0); break;
+	}
+	return -y;
+}
+
+
+double __lgamma_r(double x, int *signgamp)
+{
+	double t,y,z,nadj,p,p1,p2,p3,q,r,w;
+	int32_t hx;
+	int i,lx,ix;
+
+	EXTRACT_WORDS(hx, lx, x);
+
+	/* purge off +-inf, NaN, +-0, tiny and negative arguments */
+	*signgamp = 1;
+	ix = hx & 0x7fffffff;
+	if (ix >= 0x7ff00000)
+		return x*x;
+	if ((ix|lx) == 0)
+		return 1.0/0.0;
+	if (ix < 0x3b900000) {  /* |x|<2**-70, return -log(|x|) */
+		if(hx < 0) {
+			*signgamp = -1;
+			return -log(-x);
+		}
+		return -log(x);
+	}
+	if (hx < 0) {
+		if (ix >= 0x43300000)  /* |x|>=2**52, must be -integer */
+			return 1.0/0.0;
+		t = sin_pi(x);
+		if (t == 0.0) /* -integer */
+			return 1.0/0.0;
+		nadj = log(pi/fabs(t*x));
+		if (t < 0.0)
+			*signgamp = -1;
+		x = -x;
+	}
+
+	/* purge off 1 and 2 */
+	if (((ix - 0x3ff00000)|lx) == 0 || ((ix - 0x40000000)|lx) == 0)
+		r = 0;
+	/* for x < 2.0 */
+	else if (ix < 0x40000000) {
+		if (ix <= 0x3feccccc) {   /* lgamma(x) = lgamma(x+1)-log(x) */
+			r = -log(x);
+			if (ix >= 0x3FE76944) {
+				y = 1.0 - x;
+				i = 0;
+			} else if (ix >= 0x3FCDA661) {
+				y = x - (tc-1.0);
+				i = 1;
+			} else {
+				y = x;
+				i = 2;
+			}
+		} else {
+			r = 0.0;
+			if (ix >= 0x3FFBB4C3) {  /* [1.7316,2] */
+				y = 2.0 - x;
+				i = 0;
+			} else if(ix >= 0x3FF3B4C4) {  /* [1.23,1.73] */
+				y = x - tc;
+				i = 1;
+			} else {
+				y = x - 1.0;
+				i = 2;
+			}
+		}
+		switch (i) {
+		case 0:
+			z = y*y;
+			p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
+			p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
+			p = y*p1+p2;
+			r += (p-0.5*y);
+			break;
+		case 1:
+			z = y*y;
+			w = z*y;
+			p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12)));    /* parallel comp */
+			p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
+			p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
+			p = z*p1-(tt-w*(p2+y*p3));
+			r += tf + p;
+			break;
+		case 2:
+			p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
+			p2 = 1.0+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
+			r += -0.5*y + p1/p2;
+		}
+	} else if (ix < 0x40200000) {  /* x < 8.0 */
+		i = (int)x;
+		y = x - (double)i;
+		p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
+		q = 1.0+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
+		r = 0.5*y+p/q;
+		z = 1.0;    /* lgamma(1+s) = log(s) + lgamma(s) */
+		switch (i) {
+		case 7: z *= y + 6.0;  /* FALLTHRU */
+		case 6: z *= y + 5.0;  /* FALLTHRU */
+		case 5: z *= y + 4.0;  /* FALLTHRU */
+		case 4: z *= y + 3.0;  /* FALLTHRU */
+		case 3: z *= y + 2.0;  /* FALLTHRU */
+			r += log(z);
+			break;
+		}
+	} else if (ix < 0x43900000) {  /* 8.0 <= x < 2**58 */
+		t = log(x);
+		z = 1.0/x;
+		y = z*z;
+		w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
+		r = (x-0.5)*(t-1.0)+w;
+	} else                         /* 2**58 <= x <= inf */
+		r =  x*(log(x)-1.0);
+	if (hx < 0)
+		r = nadj - r;
+	return r;
+}
+
+weak_alias(__lgamma_r, lgamma_r);
diff --git a/system/lib/libc/musl/src/math/lgammaf.c b/system/lib/libc/musl/src/math/lgammaf.c
new file mode 100644
index 00000000..badb6dfe
--- /dev/null
+++ b/system/lib/libc/musl/src/math/lgammaf.c
@@ -0,0 +1,9 @@
+#include <math.h>
+
+extern int __signgam;
+float __lgammaf_r(float, int *);
+
+float lgammaf(float x)
+{
+	return __lgammaf_r(x, &__signgam);
+}
diff --git a/system/lib/libc/musl/src/math/lgammaf_r.c b/system/lib/libc/musl/src/math/lgammaf_r.c
new file mode 100644
index 00000000..dc65bace
--- /dev/null
+++ b/system/lib/libc/musl/src/math/lgammaf_r.c
@@ -0,0 +1,249 @@
+/* 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
+two23= 8.3886080000e+06, /* 0x4b000000 */
+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 */
+
+static float sin_pif(float x)
+{
+	float y,z;
+	int n,ix;
+
+	GET_FLOAT_WORD(ix, x);
+	ix &= 0x7fffffff;
+
+	if(ix < 0x3e800000)
+		return __sindf(pi*x);
+
+	y = -x;  /* negative x is assumed */
+
+	/*
+	 * argument reduction, make sure inexact flag not raised if input
+	 * is an integer
+	 */
+	z = floorf(y);
+	if (z != y) {   /* inexact anyway */
+		y *= 0.5f;
+		y  = 2.0f*(y - floorf(y));   /* y = |x| mod 2.0 */
+		n  = (int)(y*4.0f);
+	} else {
+		if (ix >= 0x4b800000) {
+			y = 0.0f;  /* y must be even */
+			n = 0;
+		} else {
+			if (ix < 0x4b000000)
+				z = y + two23;  /* exact */
+			GET_FLOAT_WORD(n, z);
+			n &= 1;
+			y = n;
+			n <<= 2;
+		}
+	}
+	switch (n) {
+	case 0:  y =  __sindf(pi*y); break;
+	case 1:
+	case 2:  y =  __cosdf(pi*(0.5f - y)); break;
+	case 3:
+	case 4:  y =  __sindf(pi*(1.0f - y)); break;
+	case 5:
+	case 6:  y = -__cosdf(pi*(y - 1.5f)); break;
+	default: y =  __sindf(pi*(y - 2.0f)); break;
+	}
+	return -y;
+}
+
+
+float __lgammaf_r(float x, int *signgamp)
+{
+	float t,y,z,nadj,p,p1,p2,p3,q,r,w;
+	int32_t hx;
+	int i,ix;
+
+	GET_FLOAT_WORD(hx, x);
+
+	/* purge off +-inf, NaN, +-0, tiny and negative arguments */
+	*signgamp = 1;
+	ix = hx & 0x7fffffff;
+	if (ix >= 0x7f800000)
+		return x*x;
+	if (ix == 0)
+		return 1.0f/0.0f;
+	if (ix < 0x35000000) {  /* |x| < 2**-21, return -log(|x|) */
+		if (hx < 0) {
+			*signgamp = -1;
+			return -logf(-x);
+		}
+		return -logf(x);
+	}
+	if (hx < 0) {
+		if (ix >= 0x4b000000)  /* |x| >= 2**23, must be -integer */
+			return 1.0f/0.0f;
+		t = sin_pif(x);
+		if (t == 0.0f) /* -integer */
+			return 1.0f/0.0f;
+		nadj = logf(pi/fabsf(t*x));
+		if (t < 0.0f)
+			*signgamp = -1;
+		x = -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 (hx < 0)
+		r = nadj - r;
+	return r;
+}
+
+weak_alias(__lgammaf_r, lgammaf_r);
diff --git a/system/lib/libc/musl/src/math/lgammal.c b/system/lib/libc/musl/src/math/lgammal.c
new file mode 100644
index 00000000..cc4895eb
--- /dev/null
+++ b/system/lib/libc/musl/src/math/lgammal.c
@@ -0,0 +1,386 @@
+/* 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,
+two63 = 9.223372036854775808e18L,
+
+/* 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;
+
+static long double sin_pi(long double x)
+{
+	union ldshape u = {x};
+	uint32_t ix = (u.i.se & 0x7fffU)<<16 | u.i.m>>48;
+	long double y, z;
+	int n;
+
+	if (ix < 0x3ffd8000)  /* 0.25 */
+		return sinl(pi * x);
+	y = -x;  /* x is assume negative */
+
+	/*
+	 * argument reduction, make sure inexact flag not raised if input
+	 * is an integer
+	 */
+	z = floorl(y);
+	if (z != y) {  /* inexact anyway */
+		y *= 0.5;
+		y = 2.0*(y - floorl(y));/* y = |x| mod 2.0 */
+		n = (int) (y*4.0);
+	} else {
+		if (ix >= 0x403f8000) {  /* 2^64 */
+			y = 0.0;  /* y must be even */
+			n = 0;
+		} else {
+			if (ix < 0x403e8000)  /* 2^63 */
+				z = y + two63;  /* exact */
+			u.f = z;
+			n = u.i.m & 1;
+			y = n;
+			n <<= 2;
+		}
+	}
+
+	switch (n) {
+	case 0:
+		y = sinl(pi * y);
+		break;
+	case 1:
+	case 2:
+		y = cosl(pi * (0.5 - y));
+		break;
+	case 3:
+	case 4:
+		y = sinl(pi * (1.0 - y));
+		break;
+	case 5:
+	case 6:
+		y = -cosl(pi * (y - 1.5));
+		break;
+	default:
+		y = sinl(pi * (y - 2.0));
+		break;
+	}
+	return -y;
+}
+
+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, and negative arguments */
+	if (ix >= 0x7fff0000)
+		return x * x;
+	if (x == 0) {
+		*sg -= 2*sign;
+		return 1.0 / fabsl(x);
+	}
+	if (ix < 0x3fc08000) {  /* |x|<2**-63, return -log(|x|) */
+		if (sign) {
+			*sg = -1;
+			return -logl(-x);
+		}
+		return -logl(x);
+	}
+	if (sign) {
+		t = sin_pi (x);
+		if (t == 0.0)
+			return 1.0 / fabsl(t); /* -integer */
+		nadj = logl(pi / fabsl(t * x));
+		if (t < 0.0)
+			*sg = -1;
+		x = -x;
+	}
+
+	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;
+		t = 0.0;
+		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
+
+extern int __signgam;
+
+long double lgammal(long double x)
+{
+	return __lgammal_r(x, &__signgam);
+}
+
+weak_alias(__lgammal_r, lgammal_r);
diff --git a/system/lib/libc/musl/src/math/signgam.c b/system/lib/libc/musl/src/math/signgam.c
new file mode 100644
index 00000000..cd728001
--- /dev/null
+++ b/system/lib/libc/musl/src/math/signgam.c
@@ -0,0 +1,6 @@
+#include <math.h>
+#include "libc.h"
+
+int __signgam = 0;
+
+weak_alias(__signgam, signgam);
diff --git a/system/lib/libc/musl/src/math/tgamma.c b/system/lib/libc/musl/src/math/tgamma.c
new file mode 100644
index 00000000..f91af735
--- /dev/null
+++ b/system/lib/libc/musl/src/math/tgamma.c
@@ -0,0 +1,222 @@
+/*
+"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 > 0 && isnormal(x) assumption */
+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:
+		/* sin(0-x) and -sin(x) have different sign at 0 */
+		return __sin(0-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)
+{
+	double absx, y, dy, z, r;
+
+	/* special cases */
+	if (!isfinite(x))
+		/* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
+		return x + INFINITY;
+
+	/* integer arguments */
+	/* raise inexact when non-integer */
+	if (x == floor(x)) {
+		if (x == 0)
+			/* tgamma(+-0)=+-inf with divide-by-zero */
+			return 1/x;
+		if (x < 0)
+			return 0/0.0;
+		if (x <= sizeof fact/sizeof *fact)
+			return fact[(int)x - 1];
+	}
+
+	absx = fabs(x);
+
+	/* x ~ 0: tgamma(x) ~ 1/x */
+	if (absx < 0x1p-54)
+		return 1/x;
+
+	/* x >= 172: tgamma(x)=inf with overflow */
+	/* x =< -184: tgamma(x)=+-0 with underflow */
+	if (absx >= 184) {
+		if (x < 0) {
+			FORCE_EVAL((float)(0x1p-126/x));
+			if (floor(x) * 0.5 == floor(x * 0.5))
+				return 0;
+			return -0.0;
+		}
+		x *= 0x1p1023;
+		return 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 */
+		r = -pi / (sinpi(absx) * absx * r);
+		dy = -dy;
+		z = -z;
+	}
+	r += dy * (gmhalf+0.5) * r / y;
+	z = pow(y, 0.5*z);
+	r = r * z * z;
+	return r;
+}
+
+#if 0
+double __lgamma_r(double x, int *sign)
+{
+	double r, absx, z, zz, w;
+
+	*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
diff --git a/system/lib/libc/musl/src/math/tgammaf.c b/system/lib/libc/musl/src/math/tgammaf.c
new file mode 100644
index 00000000..b4ca51c9
--- /dev/null
+++ b/system/lib/libc/musl/src/math/tgammaf.c
@@ -0,0 +1,6 @@
+#include <math.h>
+
+float tgammaf(float x)
+{
+	return tgamma(x);
+}
diff --git a/system/lib/libc/musl/src/math/tgammal.c b/system/lib/libc/musl/src/math/tgammal.c
new file mode 100644
index 00000000..5c1a02a6
--- /dev/null
+++ b/system/lib/libc/musl/src/math/tgammal.c
@@ -0,0 +1,275 @@
+/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_tgammal.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.
+ */
+/*
+ *      Gamma function
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, tgammal();
+ *
+ * y = tgammal( x );
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns gamma function of the argument.  The result is
+ * correctly signed.
+ *
+ * Arguments |x| <= 13 are reduced by recurrence and the function
+ * approximated by a rational function of degree 7/8 in the
+ * interval (2,3).  Large arguments are handled by Stirling's
+ * formula. Large negative arguments are made positive using
+ * a reflection formula.
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     -40,+40      10000       3.6e-19     7.9e-20
+ *    IEEE    -1755,+1755   10000       4.8e-18     6.5e-19
+ *
+ * Accuracy for large arguments is dominated by error in powl().
+ *
+ */
+
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double tgammal(long double x)
+{
+	return tgamma(x);
+}
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+/*
+tgamma(x+2) = tgamma(x+2) P(x)/Q(x)
+0 <= x <= 1
+Relative error
+n=7, d=8
+Peak error =  1.83e-20
+Relative error spread =  8.4e-23
+*/
+static const long double P[8] = {
+ 4.212760487471622013093E-5L,
+ 4.542931960608009155600E-4L,
+ 4.092666828394035500949E-3L,
+ 2.385363243461108252554E-2L,
+ 1.113062816019361559013E-1L,
+ 3.629515436640239168939E-1L,
+ 8.378004301573126728826E-1L,
+ 1.000000000000000000009E0L,
+};
+static const long double Q[9] = {
+-1.397148517476170440917E-5L,
+ 2.346584059160635244282E-4L,
+-1.237799246653152231188E-3L,
+-7.955933682494738320586E-4L,
+ 2.773706565840072979165E-2L,
+-4.633887671244534213831E-2L,
+-2.243510905670329164562E-1L,
+ 4.150160950588455434583E-1L,
+ 9.999999999999999999908E-1L,
+};
+
+/*
+static const long double P[] = {
+-3.01525602666895735709e0L,
+-3.25157411956062339893e1L,
+-2.92929976820724030353e2L,
+-1.70730828800510297666e3L,
+-7.96667499622741999770e3L,
+-2.59780216007146401957e4L,
+-5.99650230220855581642e4L,
+-7.15743521530849602425e4L
+};
+static const long double Q[] = {
+ 1.00000000000000000000e0L,
+-1.67955233807178858919e1L,
+ 8.85946791747759881659e1L,
+ 5.69440799097468430177e1L,
+-1.98526250512761318471e3L,
+ 3.31667508019495079814e3L,
+ 1.60577839621734713377e4L,
+-2.97045081369399940529e4L,
+-7.15743521530849602412e4L
+};
+*/
+#define MAXGAML 1755.455L
+/*static const long double LOGPI = 1.14472988584940017414L;*/
+
+/* Stirling's formula for the gamma function
+tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
+z(x) = x
+13 <= x <= 1024
+Relative error
+n=8, d=0
+Peak error =  9.44e-21
+Relative error spread =  8.8e-4
+*/
+static const long double STIR[9] = {
+ 7.147391378143610789273E-4L,
+-2.363848809501759061727E-5L,
+-5.950237554056330156018E-4L,
+ 6.989332260623193171870E-5L,
+ 7.840334842744753003862E-4L,
+-2.294719747873185405699E-4L,
+-2.681327161876304418288E-3L,
+ 3.472222222230075327854E-3L,
+ 8.333333333333331800504E-2L,
+};
+
+#define MAXSTIR 1024.0L
+static const long double SQTPI = 2.50662827463100050242E0L;
+
+/* 1/tgamma(x) = z P(z)
+ * z(x) = 1/x
+ * 0 < x < 0.03125
+ * Peak relative error 4.2e-23
+ */
+static const long double S[9] = {
+-1.193945051381510095614E-3L,
+ 7.220599478036909672331E-3L,
+-9.622023360406271645744E-3L,
+-4.219773360705915470089E-2L,
+ 1.665386113720805206758E-1L,
+-4.200263503403344054473E-2L,
+-6.558780715202540684668E-1L,
+ 5.772156649015328608253E-1L,
+ 1.000000000000000000000E0L,
+};
+
+/* 1/tgamma(-x) = z P(z)
+ * z(x) = 1/x
+ * 0 < x < 0.03125
+ * Peak relative error 5.16e-23
+ * Relative error spread =  2.5e-24
+ */
+static const long double SN[9] = {
+ 1.133374167243894382010E-3L,
+ 7.220837261893170325704E-3L,
+ 9.621911155035976733706E-3L,
+-4.219773343731191721664E-2L,
+-1.665386113944413519335E-1L,
+-4.200263503402112910504E-2L,
+ 6.558780715202536547116E-1L,
+ 5.772156649015328608727E-1L,
+-1.000000000000000000000E0L,
+};
+
+static const long double PIL = 3.1415926535897932384626L;
+
+/* Gamma function computed by Stirling's formula.
+ */
+static long double stirf(long double x)
+{
+	long double y, w, v;
+
+	w = 1.0/x;
+	/* For large x, use rational coefficients from the analytical expansion.  */
+	if (x > 1024.0)
+		w = (((((6.97281375836585777429E-5L * w
+		 + 7.84039221720066627474E-4L) * w
+		 - 2.29472093621399176955E-4L) * w
+		 - 2.68132716049382716049E-3L) * w
+		 + 3.47222222222222222222E-3L) * w
+		 + 8.33333333333333333333E-2L) * w
+		 + 1.0;
+	else
+		w = 1.0 + w * __polevll(w, STIR, 8);
+	y = expl(x);
+	if (x > MAXSTIR) { /* Avoid overflow in pow() */
+		v = powl(x, 0.5L * x - 0.25L);
+		y = v * (v / y);
+	} else {
+		y = powl(x, x - 0.5L) / y;
+	}
+	y = SQTPI * y * w;
+	return y;
+}
+
+long double tgammal(long double x)
+{
+	long double p, q, z;
+
+	if (!isfinite(x))
+		return x + INFINITY;
+
+	q = fabsl(x);
+	if (q > 13.0) {
+		if (x < 0.0) {
+			p = floorl(q);
+			z = q - p;
+			if (z == 0)
+				return 0 / z;
+			if (q > MAXGAML) {
+				z = 0;
+			} else {
+				if (z > 0.5) {
+					p += 1.0;
+					z = q - p;
+				}
+				z = q * sinl(PIL * z);
+				z = fabsl(z) * stirf(q);
+				z = PIL/z;
+			}
+			if (0.5 * p == floorl(q * 0.5))
+				z = -z;
+		} else if (x > MAXGAML) {
+			z = x * 0x1p16383L;
+		} else {
+			z = stirf(x);
+		}
+		return z;
+	}
+
+	z = 1.0;
+	while (x >= 3.0) {
+		x -= 1.0;
+		z *= x;
+	}
+	while (x < -0.03125L) {
+		z /= x;
+		x += 1.0;
+	}
+	if (x <= 0.03125L)
+		goto small;
+	while (x < 2.0) {
+		z /= x;
+		x += 1.0;
+	}
+	if (x == 2.0)
+		return z;
+
+	x -= 2.0;
+	p = __polevll(x, P, 7);
+	q = __polevll(x, Q, 8);
+	z = z * p / q;
+	return z;
+
+small:
+	/* z==1 if x was originally +-0 */
+	if (x == 0 && z != 1)
+		return x / x;
+	if (x < 0.0) {
+		x = -x;
+		q = z / (x * __polevll(x, SN, 8));
+	} else
+		q = z / (x * __polevll(x, S, 8));
+	return q;
+}
+#endif
diff --git a/system/lib/libcextra.symbols b/system/lib/libcextra.symbols
index a32bdbdb..0780847c 100644
--- a/system/lib/libcextra.symbols
+++ b/system/lib/libcextra.symbols
@@ -1,6 +1,12 @@
+         T __cosdf
          T __intscan
          W __iswctype_l
+         T __lgamma_r
+         T __lgammaf_r
+         T __lgammal_r
          T __memrchr
+         D __signgam
+         T __sindf
          T __strchrnul
          T __strxfrm_l
          W __towlower_l
@@ -45,6 +51,12 @@
          T iswupper_l
          T iswxdigit
          T iswxdigit_l
+         T lgamma
+         W lgamma_r
+         T lgammaf
+         W lgammaf_r
+         T lgammal
+         W lgammal_r
          T mblen
          T mbrlen
          T mbrtowc
@@ -61,6 +73,7 @@
          T regerror
          T regexec
          T regfree
+         W signgam
          T strcasecmp_l
          T strcasestr
          W strchrnul
@@ -74,6 +87,9 @@
          T strxfrm
          W strxfrm_l
          T swprintf
+         T tgamma
+         T tgammaf
+         T tgammal
          T towctrans
          T towctrans_l
          T towlower