Make lgamma and friends work, add test.

We can't use weak aliases on variables in emscripten, so get rid of __signgam
and just use signgam.

8 files changed, 141 insertions, 11 deletions

diff --git a/system/lib/libc/musl/src/math/lgamma.c b/system/lib/libc/musl/src/math/lgamma.c
index e25ec8e6..36006682 100644
--- a/system/lib/libc/musl/src/math/lgamma.c
+++ b/system/lib/libc/musl/src/math/lgamma.c
@@ -1,9 +1,9 @@
 #include <math.h>
 
-extern int __signgam;
+extern int signgam;
 double __lgamma_r(double, int *);
 
 double lgamma(double x)
 {
-	return __lgamma_r(x, &__signgam);
+	return __lgamma_r(x, &signgam);
 }
diff --git a/system/lib/libc/musl/src/math/lgammaf.c b/system/lib/libc/musl/src/math/lgammaf.c
index badb6dfe..628e6ade 100644
--- a/system/lib/libc/musl/src/math/lgammaf.c
+++ b/system/lib/libc/musl/src/math/lgammaf.c
@@ -1,9 +1,9 @@
 #include <math.h>
 
-extern int __signgam;
+extern int signgam;
 float __lgammaf_r(float, int *);
 
 float lgammaf(float x)
 {
-	return __lgammaf_r(x, &__signgam);
+	return __lgammaf_r(x, &signgam);
 }
diff --git a/system/lib/libc/musl/src/math/lgammal.c b/system/lib/libc/musl/src/math/lgammal.c
index cc4895eb..9686d9ee 100644
--- a/system/lib/libc/musl/src/math/lgammal.c
+++ b/system/lib/libc/musl/src/math/lgammal.c
@@ -376,11 +376,11 @@ long double __lgammal_r(long double x, int *sg) {
 }
 #endif
 
-extern int __signgam;
+extern int signgam;
 
 long double lgammal(long double x)
 {
-	return __lgammal_r(x, &__signgam);
+	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
index cd728001..ee5d70f1 100644
--- a/system/lib/libc/musl/src/math/signgam.c
+++ b/system/lib/libc/musl/src/math/signgam.c
@@ -1,6 +1,4 @@
 #include <math.h>
 #include "libc.h"
 
-int __signgam = 0;
-
-weak_alias(__signgam, signgam);
+int signgam = 0;
diff --git a/system/lib/libcextra.symbols b/system/lib/libcextra.symbols
index 0780847c..8d22ec4f 100644
--- a/system/lib/libcextra.symbols
+++ b/system/lib/libcextra.symbols
@@ -5,7 +5,6 @@
          T __lgammaf_r
          T __lgammal_r
          T __memrchr
-         D __signgam
          T __sindf
          T __strchrnul
          T __strxfrm_l
@@ -73,7 +72,7 @@
          T regerror
          T regexec
          T regfree
-         W signgam
+         D signgam
          T strcasecmp_l
          T strcasestr
          W strchrnul
diff --git a/tests/math/lgamma.in b/tests/math/lgamma.in
new file mode 100644
index 00000000..e96f5610
--- /dev/null
+++ b/tests/math/lgamma.in
@@ -0,0 +1,105 @@
+#define _BSD_SOURCE 1
+#define _XOPEN_SOURCE 700
+#include <stdint.h>
+#include <stdio.h>
+#include <fenv.h>
+#include <float.h>
+#include <math.h>
+
+#define RN 0
+#define T(...) {__FILE__, __LINE__, __VA_ARGS__},
+#define POS const char *file; int line;
+struct f_fi {POS int r; float x; float y; float dy; long long i; int e; };
+
+#define DIVBYZERO 0
+#define INEXACT 0
+#define INVALID 0
+#define OVERFLOW 0
+#define UNDERFLOW 0
+
+#define inf INFINITY
+#define nan NAN
+
+static struct f_fi t[] = {
+T(RN,   -0x1.02239f3c6a8f1p+3,   -0x1.0120f61b63d5ep+3,   0x1.89ccc4p-6,          -1, INEXACT)
+T(RN,    0x1.161868e18bc67p+2,    0x1.1ef3b263fd60bp+1,  -0x1.6d0264p-3,           1, INEXACT)
+T(RN,   -0x1.0c34b3e01e6e7p+3,   -0x1.46d73255263d9p+3,   0x1.e0ec76p-3,          -1, INEXACT)
+T(RN,   -0x1.a206f0a19dcc4p+2,   -0x1.9c91f19ac48c5p+2,   0x1.c2a38cp-2,          -1, INEXACT)
+T(RN,    0x1.288bbb0d6a1e6p+3,    0x1.65c60768fcc11p+3,   0x1.2f22c2p-2,           1, INEXACT)
+T(RN,    0x1.52efd0cd80497p-1,    0x1.3cc760be720b3p-2,   0x1.0527e2p-2,           1, INEXACT)
+T(RN,   -0x1.a05cc754481d1p-2,    0x1.4ef387fea1014p+0,  -0x1.c3b036p-2,          -1, INEXACT)
+T(RN,    0x1.1f9ef934745cbp-1,    0x1.d6f0efacc5699p-2,   0x1.c0b0a8p-2,           1, INEXACT)
+T(RN,    0x1.8c5db097f7442p-1,    0x1.6c1a14cf91533p-3,   0x1.16f4cap-5,           1, INEXACT)
+T(RN,   -0x1.5b86ea8118a0ep-1,    0x1.695b1e0a0a59ep+0,   0x1.ada69ep-2,          -1, INEXACT)
+T(RN,                  0x0p+0,                     inf,          0x0p+0,           1, DIVBYZERO)
+/* T(RN,                 -0x0p+0,                     inf,          0x0p+0,          -1, DIVBYZERO) This one fails in native as well */
+T(RN,                  0x1p+0,                  0x0p+0,          0x0p+0,           1, 0)
+T(RN,                 -0x1p+0,                     inf,          0x0p+0,           1, DIVBYZERO)
+T(RN,                  0x1p+1,                  0x0p+0,          0x0p+0,           1, 0)
+T(RN,                 -0x1p+1,                     inf,          0x0p+0,           1, DIVBYZERO)
+T(RN,                     inf,                     inf,          0x0p+0,           1, 0)
+T(RN,                    -inf,                     inf,          0x0p+0,          -1, 0)
+T(RN,                     nan,                     nan,          0x0p+0,           1, 0)
+};
+
+static int eulpf(float x)
+{
+  union { float f; uint32_t i; } u = { x };
+  int e = u.i>>23 & 0xff;
+
+  if (!e)
+    e++;
+  return e - 0x7f - 23;
+}
+
+static int checkulp(float d, int r)
+{
+  // TODO: we only care about >=1.5 ulp errors for now, should be 1.0
+  if (r == RN)
+    return fabsf(d) < 1.5;
+  return 1;
+}
+
+static float ulperrf(float got, float want, float dwant)
+{
+  if (isnan(got) && isnan(want))
+    return 0;
+  if (got == want) {
+    if (signbit(got) == signbit(want))
+      return dwant;
+    return INFINITY;
+  }
+  if (isinf(got)) {
+    got = copysignf(0x1p127, got);
+    want *= 0.5;
+  }
+  return scalbn(got - want, -eulpf(want)) + dwant;
+}
+
+int main(void)
+{
+	int yi;
+	double y;
+	float d;
+	int e, i, err = 0;
+	struct f_fi *p;
+
+	for (i = 0; i < sizeof t/sizeof *t; i++) {
+		p = t + i;
+
+		if (p->r < 0)
+			continue;
+		y = lgammaf(p->x);
+		yi = signgam;
+
+    printf("%g,%d\n", y, yi);
+
+		d = ulperrf(y, p->y, p->dy);
+		if (!checkulp(d, p->r) || (!isnan(p->x) && p->x!=-inf && !(p->e&DIVBYZERO) && yi != p->i)) {
+      /* printf("%s:%d: %d lgammaf(%g) want %g,%lld got %g,%d ulperr %.3f = %g + %g\n",
+				p->file, p->line, p->r, p->x, p->y, p->i, y, yi, d, d-p->dy, p->dy); */
+			err++;
+		}
+	}
+	return !!err;
+}
diff --git a/tests/math/lgamma.out b/tests/math/lgamma.out
new file mode 100644
index 00000000..7412ffb6
--- /dev/null
+++ b/tests/math/lgamma.out
@@ -0,0 +1,18 @@
+-8.03528,-1
+2.24181,1
+-10.2138,-1
+-6.44641,-1
+11.1804,1
+0.309354,1
+1.3084,-1
+0.459903,1
+0.177784,1
+1.41155,-1
+inf,1
+0,1
+inf,1
+0,1
+inf,1
+inf,1
+inf,1
+nan,1
diff --git a/tests/test_core.py b/tests/test_core.py
index dbc717a6..8b6d4fba 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -837,6 +837,15 @@ class T(RunnerCore): # Short name, to make it more fun to use manually on the co
       expected = open(path_from_root('tests', 'hyperbolic', 'output.txt'), 'r').read()
       self.do_run(src, expected)
 
+  def test_math_lgamma(self):
+      if self.emcc_args is None: return self.skip('requires emcc')
+      if not self.is_le32(): return self.skip('le32 needed for accurate math')
+
+      test_path = path_from_root('tests', 'math', 'lgamma')
+      src, output = (test_path + s for s in ('.in', '.out'))
+
+      self.do_run_from_file(src, output)
+
   def test_frexp(self):
       test_path = path_from_root('tests', 'core', 'test_frexp')
       src, output = (test_path + s for s in ('.in', '.out'))
@@ -1289,6 +1298,7 @@ class T(RunnerCore): # Short name, to make it more fun to use manually on the co
 
     self.do_run_from_file(src, output)
 
+
   def test_white_list_exception(self):
     if os.environ.get('EMCC_FAST_COMPILER') == '1': return self.skip('todo in fastcomp')