first commit of the new libm!

thanks to the hard work of Szabolcs Nagy (nsz), identifying the best
(from correctness and license standpoint) implementations from freebsd
and openbsd and cleaning them up! musl should now fully support c99
float and long double math functions, and has near-complete complex
math support. tgmath should also work (fully on gcc-compatible
compilers, and mostly on any c99 compiler).

based largely on commit 0376d44a890fea261506f1fc63833e7a686dca19 from
nsz's libm git repo, with some additions (dummy versions of a few
missing long double complex functions, etc.) by me.

various cleanups still need to be made, including re-adding (if
they're correct) some asm functions that were dropped.

3 files changed, 354 insertions, 9 deletions

diff --git a/include/complex.h b/include/complex.h
new file mode 100644
index 00000000..8ee70575
--- /dev/null
+++ b/include/complex.h
@@ -0,0 +1,103 @@
+#ifndef _COMPLEX_H
+#define _COMPLEX_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#define complex _Complex
+#define _Complex_I 1.0fi
+#define I _Complex_I
+
+double complex cacos(double complex);
+float complex cacosf(float complex);
+long double complex cacosl(long double complex);
+
+double complex casin(double complex);
+float complex casinf(float complex);
+long double complex casinl(long double complex);
+
+double complex catan(double complex);
+float complex catanf(float complex);
+long double complex catanl(long double complex);
+
+double complex ccos(double complex);
+float complex ccosf(float complex);
+long double complex ccosl(long double complex);
+
+double complex csin(double complex);
+float complex csinf(float complex);
+long double complex csinl(long double complex);
+
+double complex ctan(double complex);
+float complex ctanf(float complex);
+long double complex ctanl(long double complex);
+
+double complex cacosh(double complex);
+float complex cacoshf(float complex);
+long double complex cacoshl(long double complex);
+
+double complex casinh(double complex);
+float complex casinhf(float complex);
+long double complex casinhl(long double complex);
+
+double complex catanh(double complex);
+float complex catanhf(float complex);
+long double complex catanhl(long double complex);
+
+double complex ccosh(double complex);
+float complex ccoshf(float complex);
+long double complex ccoshl(long double complex);
+
+double complex csinh(double complex);
+float complex csinhf(float complex);
+long double complex csinhl(long double complex);
+
+double complex ctanh(double complex);
+float complex ctanhf(float complex);
+long double complex ctanhl(long double complex);
+
+double complex cexp(double complex);
+float complex cexpf(float complex);
+long double complex cexpl(long double complex);
+
+double complex clog(double complex);
+float complex clogf(float complex);
+long double complex clogl(long double complex);
+
+double cabs(double complex);
+float cabsf(float complex);
+long double cabsl(long double complex);
+
+double complex cpow(double complex, double complex);
+float complex cpowf(float complex, float complex);
+long double complex cpowl(long double complex, long double complex);
+
+double complex csqrt(double complex);
+float complex csqrtf(float complex);
+long double complex csqrtl(long double complex);
+
+double carg(double complex);
+float cargf(float complex);
+long double cargl(long double complex);
+
+double cimag(double complex);
+float cimagf(float complex);
+long double cimagl(long double complex);
+
+double complex conj(double complex);
+float complex conjf(float complex);
+long double complex conjl(long double complex);
+
+double complex cproj(double complex);
+float complex cprojf(float complex);
+long double complex cprojl(long double complex);
+
+double creal(double complex);
+float crealf(float complex);
+long double creall(long double complex);
+
+#ifdef __cplusplus
+}
+#endif
+#endif
diff --git a/include/math.h b/include/math.h
index ae84a731..f320b8e9 100644
--- a/include/math.h
+++ b/include/math.h
@@ -37,27 +37,53 @@ extern "C" {
 #define FP_SUBNORMAL 3
 #define FP_NORMAL    4
 
-int __fpclassifyf(float);
 int __fpclassify(double);
+int __fpclassifyf(float);
 int __fpclassifyl(long double);
 
+#define __FLOAT_BITS(f) (((union { float __f; __uint32_t __i; }){ (f) }).__i)
+#define __DOUBLE_BITS(f) (((union { double __f; __uint64_t __i; }){ (f) }).__i)
+
 #define fpclassify(x) ( \
 	sizeof(x) == sizeof(float) ? __fpclassifyf(x) : \
 	sizeof(x) == sizeof(double) ? __fpclassify(x) : \
 	__fpclassifyl(x) )
 
-#define isinf(x)    (fpclassify(x) == FP_INFINITE)
-#define isnan(x)    (fpclassify(x) == FP_NAN)
-#define isnormal(x) (fpclassify(x) == FP_NORMAL)
-#define isfinite(x) (fpclassify(x) > FP_INFINITE)
+#define isinf(x) ( \
+	sizeof(x) == sizeof(float) ? (__FLOAT_BITS(x) & 0x7fffffff) == 0x7f800000 : \
+	sizeof(x) == sizeof(double) ? (__DOUBLE_BITS(x) & (__uint64_t)-1>>1) == (__uint64_t)0x7ff<<52 : \
+	__fpclassifyl(x) == FP_INFINITE)
+
+#define isnan(x) ( \
+	sizeof(x) == sizeof(float) ? (__FLOAT_BITS(x) & 0x7fffffff) > 0x7f800000 : \
+	sizeof(x) == sizeof(double) ? (__DOUBLE_BITS(x) & (__uint64_t)-1>>1) > (__uint64_t)0x7ff<<52 : \
+	__fpclassifyl(x) == FP_NAN)
+
+#define isnormal(x) ( \
+	sizeof(x) == sizeof(float) ? ((__FLOAT_BITS(x)+0x00800000) & 0x7fffffff) >= 0x01000000 : \
+	sizeof(x) == sizeof(double) ? ((__DOUBLE_BITS(x)+((__uint64_t)1<<52)) & (__uint64_t)-1>>1) >= (__uint64_t)1<<53 : \
+	__fpclassifyl(x) == FP_NORMAL)
 
-#define isunordered(x,y) (isnan((x)) ? ((y),1) : isnan((y)))
+#define isfinite(x) ( \
+	sizeof(x) == sizeof(float) ? (__FLOAT_BITS(x) & 0x7fffffff) < 0x7f800000 : \
+	sizeof(x) == sizeof(double) ? (__DOUBLE_BITS(x) & (__uint64_t)-1>>1) < (__uint64_t)0x7ff<<52 : \
+	__fpclassifyl(x) > FP_INFINITE)
+
+int __signbit(double);
+int __signbitf(float);
+int __signbitl(long double);
+
+#define signbit(x) ( \
+	sizeof(x) == sizeof(float) ? !!(__FLOAT_BITS(x) & 0x80000000) : \
+	sizeof(x) == sizeof(double) ? !!(__DOUBLE_BITS(x) & (__uint64_t)1<<63) : \
+	__signbitl(x) )
+
+#define isunordered(x,y) (isnan((x)) ? ((void)(y),1) : isnan((y)))
 
-static
 #if __STDC_VERSION__ >= 199901L
 inline
 #endif
-int __isrel(long double __x, long double __y, int __rel)
+static int __isrel(long double __x, long double __y, int __rel)
 {
 	if (isunordered(__x, __y)) return 0;
 	if (__rel==-2) return __x < __y;
@@ -316,17 +342,46 @@ long double truncl(long double);
 #define M_2_SQRTPI      1.12837916709551257390  /* 2/sqrt(pi) */
 #define M_SQRT2         1.41421356237309504880  /* sqrt(2) */
 #define M_SQRT1_2       0.70710678118654752440  /* 1/sqrt(2) */
+
+extern int signgam;
+
+double      gamma(double);
+float       gammaf(float);
+long double gammal(long double);
+
+double      lgamma_r(double, int*);
+float       lgammaf_r(float, int*);
+long double lgammal_r(long double, int*);
+
 double      j0(double);
+float       j0f(float);
+long double j0l(long double);
+
 double      j1(double);
+float       j1f(float);
+long double j1l(long double);
+
 double      jn(int, double);
+float       jnf(int, float);
+long double jnl(int, long double);
+
 double      y0(double);
+float       y0f(float);
+long double y0l(long double);
+
 double      y1(double);
+float       y1f(float);
+long double y1l(long double);
+
 double      yn(int, double);
-extern int signgam;
+float       ynf(int, float);
+long double ynl(int, long double);
 #endif
 
 #ifdef _GNU_SOURCE
 double      scalb(double, double);
+float       scalbf(float, float);
+long double scalbl(long double, long double);
 #endif
 
 #ifdef __cplusplus
diff --git a/include/tgmath.h b/include/tgmath.h
new file mode 100644
index 00000000..52913913
--- /dev/null
+++ b/include/tgmath.h
@@ -0,0 +1,187 @@
+#ifndef _TGMATH_H
+#define _TGMATH_H
+
+/*
+the return types are only correct with gcc (__GNUC__)
+otherwise they are long double or long double complex
+
+the long double version of a function is never chosen when
+sizeof(double) == sizeof(long double)
+(but the return type is set correctly with gcc)
+*/
+
+#include <math.h>
+#include <complex.h>
+
+#define __IS_FP(x) !!((1?1:(x))/2)
+#define __IS_CX(x) (__IS_FP(x) && sizeof(x) == sizeof((x)+I))
+#define __IS_REAL(x) (__IS_FP(x) && 2*sizeof(x) == sizeof((x)+I))
+
+#define __FLT(x) (__IS_REAL(x) && sizeof(x) == sizeof(float))
+#define __LDBL(x) (__IS_REAL(x) && sizeof(x) == sizeof(long double) && sizeof(long double) != sizeof(double))
+
+#define __FLTCX(x) (__IS_CX(x) && sizeof(x) == sizeof(float complex))
+#define __DBLCX(x) (__IS_CX(x) && sizeof(x) == sizeof(double complex))
+#define __LDBLCX(x) (__IS_CX(x) && sizeof(x) == sizeof(long double complex) && sizeof(long double) != sizeof(double))
+
+/* return type */
+
+#ifdef __GNUC__
+/* cast to double when x is integral, otherwise use typeof(x) */
+#define __RETCAST(x) (__typeof__(*( \
+	0 ? (__typeof__(0 ? (double *)0 : (void *)__IS_FP(x)))0 : \
+	    (__typeof__(0 ? (__typeof__(x) *)0 : (void *)!__IS_FP(x)))0 )))
+/* 2 args case, consider complex types (for cpow) */
+#define __RETCAST_2(x, y) (__typeof__(*( \
+	0 ? (__typeof__(0 ? (double *)0 : \
+		(void *)!((!__IS_FP(x) || !__IS_FP(y)) && __FLT((x)+(y)+1.0f))))0 : \
+	0 ? (__typeof__(0 ? (double complex *)0 : \
+		(void *)!((!__IS_FP(x) || !__IS_FP(y)) && __FLTCX((x)+(y)))))0 : \
+	    (__typeof__(0 ? (__typeof__((x)+(y)) *)0 : \
+		(void *)((!__IS_FP(x) || !__IS_FP(y)) && (__FLT((x)+(y)+1.0f) || __FLTCX((x)+(y))))))0 )))
+/* 3 args case, don't consider complex types (fma only) */
+#define __RETCAST_3(x, y, z) (__typeof__(*( \
+	0 ? (__typeof__(0 ? (double *)0 : \
+		(void *)!((!__IS_FP(x) || !__IS_FP(y) || !__IS_FP(z)) && __FLT((x)+(y)+(z)+1.0f))))0 : \
+	    (__typeof__(0 ? (__typeof__((x)+(y)) *)0 : \
+		(void *)((!__IS_FP(x) || !__IS_FP(y) || !__IS_FP(z)) && __FLT((x)+(y)+(z)+1.0f))))0 )))
+/* drop complex from the type of x */
+#define __TO_REAL(x) *( \
+	0 ? (__typeof__(0 ? (double *)0 : (void *)!__DBLCX(x)))0 : \
+	0 ? (__typeof__(0 ? (float *)0 : (void *)!__FLTCX(x)))0 : \
+	0 ? (__typeof__(0 ? (long double *)0 : (void *)!__LDBLCX(x)))0 : \
+	    (__typeof__(0 ? (__typeof__(x) *)0 : (void *)__IS_CX(x)))0 )
+#else
+#define __RETCAST(x)
+#define __RETCAST_2(x, y)
+#define __RETCAST_3(x, y, z)
+#endif
+
+/* function selection */
+
+#define __tg_real(fun, x) (__RETCAST(x)( \
+	__FLT(x) ? fun ## f (x) : \
+	__LDBL(x) ? fun ## l (x) : \
+	fun(x) ))
+
+#define __tg_real_2_1(fun, x, y) (__RETCAST(x)( \
+	__FLT(x) ? fun ## f (x, y) : \
+	__LDBL(x) ? fun ## l (x, y) : \
+	fun(x, y) ))
+
+#define __tg_real_2(fun, x, y) (__RETCAST_2(x, y)( \
+	__FLT(x) && __FLT(y) ? fun ## f (x, y) : \
+	__LDBL((x)+(y)) ? fun ## l (x, y) : \
+	fun(x, y) ))
+
+#define __tg_complex(fun, x) (__RETCAST((x)+I)( \
+	__FLTCX((x)+I) && __IS_FP(x) ? fun ## f (x) : \
+	__LDBLCX((x)+I) ? fun ## l (x) : \
+	fun(x) ))
+
+#define __tg_complex_retreal(fun, x) (__RETCAST(__TO_REAL(x))( \
+	__FLTCX((x)+I) && __IS_FP(x) ? fun ## f (x) : \
+	__LDBLCX((x)+I) ? fun ## l (x) : \
+	fun(x) ))
+
+#define __tg_real_complex(fun, x) (__RETCAST(x)( \
+	__FLTCX(x) ? c ## fun ## f (x) : \
+	__DBLCX(x) ? c ## fun (x) : \
+	__LDBLCX(x) ? c ## fun ## l (x) : \
+	__FLT(x) ? fun ## f (x) : \
+	__LDBL(x) ? fun ## l (x) : \
+	fun(x) ))
+
+/* special cases */
+
+#define __tg_real_remquo(x, y, z) (__RETCAST_2(x, y)( \
+	__FLT(x) && __FLT(y) ? remquof(x, y, z) : \
+	__LDBL((x)+(y)) ? remquol(x, y, z) : \
+	remquo(x, y, z) ))
+
+#define __tg_real_fma(x, y, z) (__RETCAST_3(x, y, z)( \
+	__FLT(x) && __FLT(y) && __FLT(z) ? fmaf(x, y, z) : \
+	__LDBL((x)+(y)+(z)) ? fmal(x, y, z) : \
+	fma(x, y, z) ))
+
+#define __tg_real_complex_pow(x, y) (__RETCAST_2(x, y)( \
+	__FLTCX((x)+(y)) && __IS_FP(x) && __IS_FP(y) ? cpowf(x, y) : \
+	__FLTCX((x)+(y)) ? cpow(x, y) : \
+	__DBLCX((x)+(y)) ? cpow(x, y) : \
+	__LDBLCX((x)+(y)) ? cpowl(x, y) : \
+	__FLT(x) && __FLT(y) ? powf(x, y) : \
+	__LDBL((x)+(y)) ? powl(x, y) : \
+	pow(x, y) ))
+
+#define __tg_real_complex_fabs(x) (__RETCAST(__TO_REAL(x))( \
+	__FLTCX(x) ? cabsf(x) : \
+	__DBLCX(x) ? cabs(x) : \
+	__LDBLCX(x) ? cabsl(x) : \
+	__FLT(x) ? fabsf(x) : \
+	__LDBL(x) ? fabsl(x) : \
+	fabs(x) ))
+
+/* tg functions */
+
+#define acos(x)         __tg_real_complex(acos, (x))
+#define acosh(x)        __tg_real_complex(acosh, (x))
+#define asin(x)         __tg_real_complex(asin, (x))
+#define asinh(x)        __tg_real_complex(asinh, (x))
+#define atan(x)         __tg_real_complex(atan, (x))
+#define atan2(x,y)      __tg_real_2(atan2, (x), (y))
+#define atanh(x)        __tg_real_complex(atanh, (x))
+#define carg(x)         __tg_complex_retreal(carg, (x))
+#define cbrt(x)         __tg_real(cbrt, (x))
+#define ceil(x)         __tg_real(ceil, (x))
+#define cimag(x)        __tg_complex_retreal(cimag, (x))
+#define conj(x)         __tg_complex(conj, (x))
+#define copysign(x,y)   __tg_real_2(copysign, (x), (y))
+#define cos(x)          __tg_real_complex(cos, (x))
+#define cosh(x)         __tg_real_complex(cosh, (x))
+#define cproj(x)        __tg_complex(cproj, (x))
+#define creal(x)        __tg_complex_retreal(creal, (x))
+#define erf(x)          __tg_real(erf, (x))
+#define erfc(x)         __tg_real(erfc, (x))
+#define exp(x)          __tg_real_complex(exp, (x))
+#define exp2(x)         __tg_real(exp2, (x))
+#define expm1(x)        __tg_real(expm1, (x))
+#define fabs(x)         __tg_real_complex_fabs(x)
+#define fdim(x,y)       __tg_real_2(fdim, (x), (y))
+#define floor(x)        __tg_real(floor, (x))
+#define fma(x,y,z)      __tg_real_fma((x), (y), (z))
+#define fmax(x,y)       __tg_real_2(fmax, (x), (y))
+#define fmin(x,y)       __tg_real_2(fmin, (x), (y))
+#define fmod(x,y)       __tg_real_2(fmod, (x), (y))
+#define frexp(x,y)      __tg_real_2_1(frexp, (x), (y))
+#define hypot(x,y)      __tg_real_2(hypot, (x), (y))
+#define ilogb(x)        __tg_real(ilogb, (x))
+#define ldexp(x,y)      __tg_real_2_1(ldexp, (x), (y))
+#define lgamma(x)       __tg_real(lgamma, (x))
+#define llrint(x)       __tg_real(llrint, (x))
+#define llround(x)      __tg_real(llround, (x))
+#define log(x)          __tg_real_complex(log, (x))
+#define log10(x)        __tg_real(log10, (x))
+#define log1p(x)        __tg_real(log1p, (x))
+#define log2(x)         __tg_real(log2, (x))
+#define logb(x)         __tg_real(logb, (x))
+#define lrint(x)        __tg_real(lrint, (x))
+#define lround(x)       __tg_real(lround, (x))
+#define nearbyint(x)    __tg_real(nearbyint, (x))
+#define nextafter(x,y)  __tg_real_2(nextafter, (x), (y)
+#define nexttoward(x,y) __tg_real_2(nexttoward, (x), (y))
+#define pow(x,y)        __tg_real_complex_pow((x), (y))
+#define remainder(x,y)  __tg_real_2(remainder, (x), (y))
+#define remquo(x,y,z)   __tg_real_remquo((x), (y), (z))
+#define rint(x)         __tg_real(rint, (x))
+#define round(x)        __tg_real(round, (x))
+#define scalbln(x,y)    __tg_real_2_1(scalbln, (x), (y))
+#define scalbn(x,y)     __tg_real_2_1(scalbn, (x), (y))
+#define sin(x)          __tg_real_complex(sin, (x))
+#define sinh(x)         __tg_real_complex(sinh, (x))
+#define sqrt(x)         __tg_real_complex(sqrt, (x))
+#define tan(x)          __tg_real_complex(tan, (x))
+#define tanh(x)         __tg_real_complex(tanh, (x))
+#define tgamma(x)       __tg_real(tgamma, (x))
+#define trunc(x)        __tg_real(trunc, (x))
+
+#endif