Add lots of math functions.

12 files changed, 241 insertions, 0 deletions

diff --git a/system/lib/libc/musl/src/math/ilogb.c b/system/lib/libc/musl/src/math/ilogb.c
new file mode 100644
index 00000000..64d40154
--- /dev/null
+++ b/system/lib/libc/musl/src/math/ilogb.c
@@ -0,0 +1,26 @@
+#include <limits.h>
+#include "libm.h"
+
+int ilogb(double x)
+{
+	#pragma STDC FENV_ACCESS ON
+	union {double f; uint64_t i;} u = {x};
+	uint64_t i = u.i;
+	int e = i>>52 & 0x7ff;
+
+	if (!e) {
+		i <<= 12;
+		if (i == 0) {
+			FORCE_EVAL(0/0.0f);
+			return FP_ILOGB0;
+		}
+		/* subnormal x */
+		for (e = -0x3ff; i>>63 == 0; e--, i<<=1);
+		return e;
+	}
+	if (e == 0x7ff) {
+		FORCE_EVAL(0/0.0f);
+		return i<<12 ? FP_ILOGBNAN : INT_MAX;
+	}
+	return e - 0x3ff;
+}
diff --git a/system/lib/libc/musl/src/math/ilogbf.c b/system/lib/libc/musl/src/math/ilogbf.c
new file mode 100644
index 00000000..e23ba209
--- /dev/null
+++ b/system/lib/libc/musl/src/math/ilogbf.c
@@ -0,0 +1,26 @@
+#include <limits.h>
+#include "libm.h"
+
+int ilogbf(float x)
+{
+	#pragma STDC FENV_ACCESS ON
+	union {float f; uint32_t i;} u = {x};
+	uint32_t i = u.i;
+	int e = i>>23 & 0xff;
+
+	if (!e) {
+		i <<= 9;
+		if (i == 0) {
+			FORCE_EVAL(0/0.0f);
+			return FP_ILOGB0;
+		}
+		/* subnormal x */
+		for (e = -0x7f; i>>31 == 0; e--, i<<=1);
+		return e;
+	}
+	if (e == 0xff) {
+		FORCE_EVAL(0/0.0f);
+		return i<<9 ? FP_ILOGBNAN : INT_MAX;
+	}
+	return e - 0x7f;
+}
diff --git a/system/lib/libc/musl/src/math/ilogbl.c b/system/lib/libc/musl/src/math/ilogbl.c
new file mode 100644
index 00000000..7b1a9cf8
--- /dev/null
+++ b/system/lib/libc/musl/src/math/ilogbl.c
@@ -0,0 +1,55 @@
+#include <limits.h>
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+int ilogbl(long double x)
+{
+	return ilogb(x);
+}
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+int ilogbl(long double x)
+{
+	#pragma STDC FENV_ACCESS ON
+	union ldshape u = {x};
+	uint64_t m = u.i.m;
+	int e = u.i.se & 0x7fff;
+
+	if (!e) {
+		if (m == 0) {
+			FORCE_EVAL(0/0.0f);
+			return FP_ILOGB0;
+		}
+		/* subnormal x */
+		for (e = -0x3fff+1; m>>63 == 0; e--, m<<=1);
+		return e;
+	}
+	if (e == 0x7fff) {
+		FORCE_EVAL(0/0.0f);
+		return m<<1 ? FP_ILOGBNAN : INT_MAX;
+	}
+	return e - 0x3fff;
+}
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
+int ilogbl(long double x)
+{
+	#pragma STDC FENV_ACCESS ON
+	union ldshape u = {x};
+	int e = u.i.se & 0x7fff;
+
+	if (!e) {
+		if (x == 0) {
+			FORCE_EVAL(0/0.0f);
+			return FP_ILOGB0;
+		}
+		/* subnormal x */
+		x *= 0x1p120;
+		return ilogbl(x) - 120;
+	}
+	if (e == 0x7fff) {
+		FORCE_EVAL(0/0.0f);
+		u.i.se = 0;
+		return u.f ? FP_ILOGBNAN : INT_MAX;
+	}
+	return e - 0x3fff;
+}
+#endif
diff --git a/system/lib/libc/musl/src/math/ldexp.c b/system/lib/libc/musl/src/math/ldexp.c
new file mode 100644
index 00000000..f4d1cd6a
--- /dev/null
+++ b/system/lib/libc/musl/src/math/ldexp.c
@@ -0,0 +1,6 @@
+#include <math.h>
+
+double ldexp(double x, int n)
+{
+	return scalbn(x, n);
+}
diff --git a/system/lib/libc/musl/src/math/ldexpf.c b/system/lib/libc/musl/src/math/ldexpf.c
new file mode 100644
index 00000000..3bad5f39
--- /dev/null
+++ b/system/lib/libc/musl/src/math/ldexpf.c
@@ -0,0 +1,6 @@
+#include <math.h>
+
+float ldexpf(float x, int n)
+{
+	return scalbnf(x, n);
+}
diff --git a/system/lib/libc/musl/src/math/ldexpl.c b/system/lib/libc/musl/src/math/ldexpl.c
new file mode 100644
index 00000000..fd145ccc
--- /dev/null
+++ b/system/lib/libc/musl/src/math/ldexpl.c
@@ -0,0 +1,6 @@
+#include <math.h>
+
+long double ldexpl(long double x, int n)
+{
+	return scalbnl(x, n);
+}
diff --git a/system/lib/libc/musl/src/math/logb.c b/system/lib/libc/musl/src/math/logb.c
new file mode 100644
index 00000000..7f8bdfae
--- /dev/null
+++ b/system/lib/libc/musl/src/math/logb.c
@@ -0,0 +1,17 @@
+#include <math.h>
+
+/*
+special cases:
+	logb(+-0) = -inf, and raise divbyzero
+	logb(+-inf) = +inf
+	logb(nan) = nan
+*/
+
+double logb(double x)
+{
+	if (!isfinite(x))
+		return x * x;
+	if (x == 0)
+		return -1/(x*x);
+	return ilogb(x);
+}
diff --git a/system/lib/libc/musl/src/math/logbf.c b/system/lib/libc/musl/src/math/logbf.c
new file mode 100644
index 00000000..a0a0b5ed
--- /dev/null
+++ b/system/lib/libc/musl/src/math/logbf.c
@@ -0,0 +1,10 @@
+#include <math.h>
+
+float logbf(float x)
+{
+	if (!isfinite(x))
+		return x * x;
+	if (x == 0)
+		return -1/(x*x);
+	return ilogbf(x);
+}
diff --git a/system/lib/libc/musl/src/math/logbl.c b/system/lib/libc/musl/src/math/logbl.c
new file mode 100644
index 00000000..962973a7
--- /dev/null
+++ b/system/lib/libc/musl/src/math/logbl.c
@@ -0,0 +1,16 @@
+#include <math.h>
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double logbl(long double x)
+{
+	return logb(x);
+}
+#else
+long double logbl(long double x)
+{
+	if (!isfinite(x))
+		return x * x;
+	if (x == 0)
+		return -1/(x*x);
+	return ilogbl(x);
+}
+#endif
diff --git a/system/lib/libc/musl/src/math/scalbn.c b/system/lib/libc/musl/src/math/scalbn.c
new file mode 100644
index 00000000..530e07c7
--- /dev/null
+++ b/system/lib/libc/musl/src/math/scalbn.c
@@ -0,0 +1,31 @@
+#include <math.h>
+#include <stdint.h>
+
+double scalbn(double x, int n)
+{
+	union {double f; uint64_t i;} u;
+	double_t y = x;
+
+	if (n > 1023) {
+		y *= 0x1p1023;
+		n -= 1023;
+		if (n > 1023) {
+			y *= 0x1p1023;
+			n -= 1023;
+			if (n > 1023)
+				n = 1023;
+		}
+	} else if (n < -1022) {
+		y *= 0x1p-1022;
+		n += 1022;
+		if (n < -1022) {
+			y *= 0x1p-1022;
+			n += 1022;
+			if (n < -1022)
+				n = -1022;
+		}
+	}
+	u.i = (uint64_t)(0x3ff+n)<<52;
+	x = y * u.f;
+	return x;
+}
diff --git a/system/lib/libc/musl/src/math/scalbnf.c b/system/lib/libc/musl/src/math/scalbnf.c
new file mode 100644
index 00000000..0b62c3c7
--- /dev/null
+++ b/system/lib/libc/musl/src/math/scalbnf.c
@@ -0,0 +1,31 @@
+#include <math.h>
+#include <stdint.h>
+
+float scalbnf(float x, int n)
+{
+	union {float f; uint32_t i;} u;
+	float_t y = x;
+
+	if (n > 127) {
+		y *= 0x1p127f;
+		n -= 127;
+		if (n > 127) {
+			y *= 0x1p127f;
+			n -= 127;
+			if (n > 127)
+				n = 127;
+		}
+	} else if (n < -126) {
+		y *= 0x1p-126f;
+		n += 126;
+		if (n < -126) {
+			y *= 0x1p-126f;
+			n += 126;
+			if (n < -126)
+				n = -126;
+		}
+	}
+	u.i = (uint32_t)(0x7f+n)<<23;
+	x = y * u.f;
+	return x;
+}
diff --git a/system/lib/libcextra.symbols b/system/lib/libcextra.symbols
index 8d22ec4f..db1d3dcd 100644
--- a/system/lib/libcextra.symbols
+++ b/system/lib/libcextra.symbols
@@ -24,6 +24,9 @@
          T iconv
          T iconv_close
          T iconv_open
+         T ilogb
+         T ilogbf
+         T ilogbl
          T iswalnum
          T iswalnum_l
          T iswalpha
@@ -50,12 +53,18 @@
          T iswupper_l
          T iswxdigit
          T iswxdigit_l
+         T ldexp
+         T ldexpf
+         T ldexpl
          T lgamma
          W lgamma_r
          T lgammaf
          W lgammaf_r
          T lgammal
          W lgammal_r
+         T logb
+         T logbf
+         T logbl
          T mblen
          T mbrlen
          T mbrtowc
@@ -72,6 +81,8 @@
          T regerror
          T regexec
          T regfree
+         T scalbn
+         T scalbnf
          D signgam
          T strcasecmp_l
          T strcasestr