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.

2 files changed, 323 insertions, 0 deletions

diff --git a/src/internal/libm.h b/src/internal/libm.h
new file mode 100644
index 00000000..021c4e2a
--- /dev/null
+++ b/src/internal/libm.h
@@ -0,0 +1,186 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/math_private.h */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#ifndef _LIBM_H
+#define _LIBM_H
+
+#include <stdint.h>
+#include <float.h>
+#include <math.h>
+#include <complex.h>
+
+#include "longdbl.h"
+
+union fshape {
+	float value;
+	uint32_t bits;
+};
+
+union dshape {
+	double value;
+	uint64_t bits;
+};
+
+/* Get two 32 bit ints from a double.  */
+#define EXTRACT_WORDS(hi,lo,d)                                  \
+do {                                                            \
+  union dshape __u;                                             \
+  __u.value = (d);                                              \
+  (hi) = __u.bits >> 32;                                        \
+  (lo) = (uint32_t)__u.bits;                                    \
+} while (0)
+
+/* Get a 64 bit int from a double.  */
+#define EXTRACT_WORD64(i,d)                                     \
+do {                                                            \
+  union dshape __u;                                             \
+  __u.value = (d);                                              \
+  (i) = __u.bits;                                               \
+} while (0)
+
+/* Get the more significant 32 bit int from a double.  */
+#define GET_HIGH_WORD(i,d)                                      \
+do {                                                            \
+  union dshape __u;                                             \
+  __u.value = (d);                                              \
+  (i) = __u.bits >> 32;                                         \
+} while (0)
+
+/* Get the less significant 32 bit int from a double.  */
+#define GET_LOW_WORD(i,d)                                       \
+do {                                                            \
+  union dshape __u;                                             \
+  __u.value = (d);                                              \
+  (i) = (uint32_t)__u.bits;                                     \
+} while (0)
+
+/* Set a double from two 32 bit ints.  */
+#define INSERT_WORDS(d,hi,lo)                                   \
+do {                                                            \
+  union dshape __u;                                             \
+  __u.bits = ((uint64_t)(hi) << 32) | (uint32_t)(lo);           \
+  (d) = __u.value;                                              \
+} while (0)
+
+/* Set a double from a 64 bit int.  */
+#define INSERT_WORD64(d,i)                                      \
+do {                                                            \
+  union dshape __u;                                             \
+  __u.bits = (i);                                               \
+  (d) = __u.value;                                              \
+} while (0)
+
+/* Set the more significant 32 bits of a double from an int.  */
+#define SET_HIGH_WORD(d,hi)                                     \
+do {                                                            \
+  union dshape __u;                                             \
+  __u.value = (d);                                              \
+  __u.bits &= 0xffffffff;                                       \
+  __u.bits |= (uint64_t)(hi) << 32;                             \
+  (d) = __u.value;                                              \
+} while (0)
+
+/* Set the less significant 32 bits of a double from an int.  */
+#define SET_LOW_WORD(d,lo)                                      \
+do {                                                            \
+  union dshape __u;                                             \
+  __u.value = (d);                                              \
+  __u.bits &= 0xffffffff00000000ull;                            \
+  __u.bits |= (uint32_t)(lo);                                   \
+  (d) = __u.value;                                              \
+} while (0)
+
+/* Get a 32 bit int from a float.  */
+#define GET_FLOAT_WORD(i,d)                                     \
+do {                                                            \
+  union fshape __u;                                             \
+  __u.value = (d);                                              \
+  (i) = __u.bits;                                               \
+} while (0)
+
+/* Set a float from a 32 bit int.  */
+#define SET_FLOAT_WORD(d,i)                                     \
+do {                                                            \
+  union fshape __u;                                             \
+  __u.bits = (i);                                               \
+  (d) = __u.value;                                              \
+} while (0)
+
+/* fdlibm kernel functions */
+
+int    __rem_pio2_large(double*,double*,int,int,int);
+
+int    __rem_pio2(double,double*);
+double __sin(double,double,int);
+double __cos(double,double);
+double __tan(double,double,int);
+double __expo2(double);
+double complex __ldexp_cexp(double complex,int);
+
+int    __rem_pio2f(float,double*);
+float  __sindf(double);
+float  __cosdf(double);
+float  __tandf(double,int);
+float  __expo2f(float);
+float complex __ldexp_cexpf(float complex,int);
+
+long double __sinl(long double, long double, int);
+long double __cosl(long double, long double);
+long double __tanl(long double, long double, int);
+
+/* polynomial evaluation */
+long double __polevll(long double, long double *, int);
+long double __p1evll(long double, long double *, int);
+
+// FIXME: not needed when -fexcess-precision=standard is supported (>=gcc4.5)
+/*
+ * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
+ */
+#if 1
+#define STRICT_ASSIGN(type, lval, rval) do {    \
+        volatile type __v = (rval);             \
+        (lval) = __v;                           \
+} while (0)
+#else
+#define STRICT_ASSIGN(type, lval, rval) ((lval) = (type)(rval))
+#endif
+
+
+/* complex */
+
+union dcomplex {
+	double complex z;
+	double a[2];
+};
+union fcomplex {
+	float complex z;
+	float a[2];
+};
+union lcomplex {
+	long double complex z;
+	long double a[2];
+};
+
+// FIXME: move to complex.h ?
+#define creal(z) ((double)(z))
+#define crealf(z) ((float)(z))
+#define creall(z) ((long double)(z))
+#define cimag(z) ((union dcomplex){(z)}.a[1])
+#define cimagf(z) ((union fcomplex){(z)}.a[1])
+#define cimagl(z) ((union lcomplex){(z)}.a[1])
+
+/* x + y*I is not supported properly by gcc */
+#define cpack(x,y) ((union dcomplex){.a={(x),(y)}}.z)
+#define cpackf(x,y) ((union fcomplex){.a={(x),(y)}}.z)
+#define cpackl(x,y) ((union lcomplex){.a={(x),(y)}}.z)
+
+#endif
diff --git a/src/internal/longdbl.h b/src/internal/longdbl.h
new file mode 100644
index 00000000..25ec8021
--- /dev/null
+++ b/src/internal/longdbl.h
@@ -0,0 +1,137 @@
+#ifndef _LDHACK_H
+#define _LDHACK_H
+
+#include <float.h>
+#include <stdint.h>
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+union ldshape {
+	long double value;
+	struct {
+		uint64_t m;
+		uint16_t exp:15;
+		uint16_t sign:1;
+		uint16_t pad;
+	} bits;
+};
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
+union ldshape {
+	long double value;
+	struct {
+		uint64_t mlo;
+		uint64_t mhi:48;
+		uint16_t exp:15;
+		uint16_t sign:1;
+	} bits;
+};
+#else
+#error Unsupported long double representation
+#endif
+
+
+// FIXME: hacks to make freebsd+openbsd long double code happy
+
+// union and macros for freebsd
+
+#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+
+union IEEEl2bits {
+	long double e;
+	struct {
+		uint32_t manl:32;
+		uint32_t manh:32;
+		uint32_t exp:15;
+		uint32_t sign:1;
+		uint32_t pad:16;
+	} bits;
+	struct {
+		uint64_t man:64;
+		uint32_t expsign:16;
+		uint32_t pad:16;
+	} xbits;
+};
+
+#define LDBL_MANL_SIZE 32
+#define LDBL_MANH_SIZE 32
+#define LDBL_NBIT (1ull << LDBL_MANH_SIZE-1)
+#undef LDBL_IMPLICIT_NBIT
+#define mask_nbit_l(u) ((u).bits.manh &= ~LDBL_NBIT)
+
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
+/*
+// ld128 float.h
+//#define LDBL_MAX 1.189731495357231765085759326628007016E+4932L
+#define LDBL_MAX 0x1.ffffffffffffffffffffffffffffp+16383
+#define LDBL_MAX_EXP 16384
+#define LDBL_HAS_INFINITY 1
+//#define LDBL_MIN 3.362103143112093506262677817321752603E-4932L
+#define LDBL_MIN 0x1p-16382
+#define LDBL_HAS_QUIET_NAN 1
+#define LDBL_HAS_DENORM 1
+//#define LDBL_EPSILON 1.925929944387235853055977942584927319E-34L
+#define LDBL_EPSILON 0x1p-112
+#define LDBL_MANT_DIG 113
+#define LDBL_MIN_EXP (-16381)
+#define LDBL_MAX_10_EXP 4932
+#define LDBL_DENORM_MIN 0x0.0000000000000000000000000001p-16381
+#define LDBL_MIN_10_EXP (-4931)
+#define LDBL_DIG 33
+*/
+
+union IEEEl2bits {
+	long double e;
+	struct {
+		uint64_t manl:64;
+		uint64_t manh:48;
+		uint32_t exp:15;
+		uint32_t sign:1;
+	} bits;
+	struct {
+		uint64_t unused0:64;
+		uint64_t unused1:48;
+		uint32_t expsign:16;
+	} xbits;
+};
+
+#define LDBL_MANL_SIZE 64
+#define LDBL_MANH_SIZE 48
+#define LDBL_NBIT (1ull << LDBL_MANH_SIZE)
+#define LDBL_IMPLICIT_NBIT 1
+#define mask_nbit_l(u)
+
+#endif
+
+
+// macros for openbsd
+
+#define GET_LDOUBLE_WORDS(se,mh,ml, f) do{ \
+	union IEEEl2bits u; \
+	u.e = (f); \
+	(se) = u.xbits.expsign; \
+	(mh) = u.bits.manh; \
+	(ml) = u.bits.manl; \
+}while(0)
+
+#define SET_LDOUBLE_WORDS(f,  se,mh,ml) do{ \
+	union IEEEl2bits u; \
+	u.xbits.expsign = (se); \
+	u.bits.manh = (mh); \
+	u.bits.manl = (ml); \
+	(f) = u.e; \
+}while(0)
+
+#define GET_LDOUBLE_EXP(se, f) do{ \
+	union IEEEl2bits u; \
+	u.e = (f); \
+	(se) = u.xbits.expsign; \
+}while(0)
+
+#define SET_LDOUBLE_EXP(f, se) do{ \
+	union IEEEl2bits u; \
+	u.e = (f); \
+	u.xbits.expsign = (se); \
+	(f) = u.e; \
+}while(0)
+
+#endif