math: add fenv test and new 'add' function in gen

diff --git a/src/math/fenv.c b/src/math/fenv.c
new file mode 100644 (file)
index 0000000..c221d7f
--- /dev/null
+++ b/src/math/fenv.c
@@ -0,0 +1,173 @@
+#include <stdint.h>
+#include <stdio.h>
+#include <stdarg.h>
+#include "util.h"
+
+static int test_status;
+
+#define error(...) print(__FILE__, __LINE__, __VA_ARGS__)
+static void print(char *f, int l, char *fmt, ...)
+{
+       test_status = 1;
+       va_list ap;
+       printf("%s:%d: ", f, l);
+       va_start(ap, fmt);
+       vprintf(fmt, ap);
+       va_end(ap);
+}
+
+#define F(n) {#n, n}
+
+static struct {
+       char *name;
+       int i;
+} te[] = {
+#ifdef FE_DIVBYZERO
+       F(FE_DIVBYZERO),
+#endif
+#ifdef FE_INEXACT
+       F(FE_INEXACT),
+#endif
+#ifdef FE_INVALID
+       F(FE_INVALID),
+#endif
+#ifdef FE_OVERFLOW
+       F(FE_OVERFLOW),
+#endif
+#ifdef FE_UNDERFLOW
+       F(FE_UNDERFLOW),
+#endif
+};
+
+static void test_except()
+{
+       #pragma STDC FENV_ACCESS ON
+       int i,r;
+
+       for (i=0; i < sizeof te/sizeof*te; i++) {
+               feclearexcept(FE_ALL_EXCEPT);
+
+               r = feraiseexcept(te[i].i);
+               if (r)
+                       error("feraiseexcept(%s) returned %d\n", te[i].name, r);
+               r = fetestexcept(FE_ALL_EXCEPT);
+               if (r != te[i].i)
+                       error("feraiseexcept(%s) want %d got %d\n",
+                               te[i].name, te[i].i, r);
+       }
+}
+
+static struct {
+       char *name;
+       int i;
+} tr[] = {
+       F(FE_TONEAREST),
+#ifdef FE_UPWARD
+       F(FE_UPWARD),
+#endif
+#ifdef FE_DOWNWARD
+       F(FE_DOWNWARD),
+#endif
+#ifdef FE_TOWARDZERO
+       F(FE_TOWARDZERO),
+#endif
+};
+
+static void test_round()
+{
+       #pragma STDC FENV_ACCESS ON
+       int i,r;
+
+       for (i=0; i < sizeof tr/sizeof*tr; i++) {
+               if (tr[i].i < 0)
+                       error("%s (%d) < 0\n", tr[i].name, tr[i].i);
+               for (r=0; r < i; r++)
+                       if (tr[r].i == tr[i].i)
+                               error("%s (%d) == %s (%d)\n",
+                                       tr[r].name, tr[r].i, tr[i].name, tr[i].i);
+       }
+
+       for (i=0; i < sizeof tr/sizeof*tr; i++) {
+               r = fesetround(tr[i].i);
+               if (r != 0)
+                       error("fesetround %d\n", r);
+               r = fegetround();
+               if (r != tr[i].i)
+                       error("fegetround %x wanted %x\n", r, tr[i].i);
+       }
+}
+
+/* ieee double precision add operation */
+static struct dd_d t[] = {
+T(RN,                  0x1p+0,                 0x1p-52,    0x1.0000000000001p+0,          0x0p+0, 0)
+T(RN,                  0x1p+0,                 0x1p-53,                  0x1p+0,         -0x1p-1, INEXACT)
+T(RN,                  0x1p+0,              0x1.01p-53,    0x1.0000000000001p+0,       0x1.fep-2, INEXACT)
+T(RN,                  0x1p+0,                -0x1p-54,                  0x1p+0,          0x1p-2, INEXACT)
+T(RN,                  0x1p+0,             -0x1.01p-54,    0x1.fffffffffffffp-1,      -0x1.fep-2, INEXACT)
+T(RN,                 -0x1p+0,                -0x1p-53,                 -0x1p+0,          0x1p-1, INEXACT)
+T(RN,                 -0x1p+0,             -0x1.01p-53,   -0x1.0000000000001p+0,      -0x1.fep-2, INEXACT)
+T(RN,                 -0x1p+0,                 0x1p-54,                 -0x1p+0,         -0x1p-2, INEXACT)
+T(RN,                 -0x1p+0,              0x1.01p-54,   -0x1.fffffffffffffp-1,       0x1.fep-2, INEXACT)
+
+T(RU,                  0x1p+0,                 0x1p-52,    0x1.0000000000001p+0,          0x0p+0, 0)
+T(RU,                  0x1p+0,                 0x1p-53,    0x1.0000000000001p+0,          0x1p-1, INEXACT)
+T(RU,                  0x1p+0,              0x1.01p-53,    0x1.0000000000001p+0,       0x1.fep-2, INEXACT)
+T(RU,                  0x1p+0,                -0x1p-54,                  0x1p+0,          0x1p-2, INEXACT)
+T(RU,                  0x1p+0,             -0x1.01p-54,                  0x1p+0,       0x1.01p-2, INEXACT)
+T(RU,                 -0x1p+0,                -0x1p-53,                 -0x1p+0,          0x1p-1, INEXACT)
+T(RU,                 -0x1p+0,             -0x1.01p-53,                 -0x1p+0,       0x1.01p-1, INEXACT)
+T(RU,                 -0x1p+0,                 0x1p-54,   -0x1.fffffffffffffp-1,          0x1p-1, INEXACT)
+T(RU,                 -0x1p+0,              0x1.01p-54,   -0x1.fffffffffffffp-1,       0x1.fep-2, INEXACT)
+
+T(RD,                  0x1p+0,                 0x1p-52,    0x1.0000000000001p+0,          0x0p+0, 0)
+T(RD,                  0x1p+0,                 0x1p-53,                  0x1p+0,         -0x1p-1, INEXACT)
+T(RD,                  0x1p+0,              0x1.01p-53,                  0x1p+0,      -0x1.01p-1, INEXACT)
+T(RD,                  0x1p+0,                -0x1p-54,    0x1.fffffffffffffp-1,         -0x1p-1, INEXACT)
+T(RD,                  0x1p+0,             -0x1.01p-54,    0x1.fffffffffffffp-1,      -0x1.fep-2, INEXACT)
+T(RD,                 -0x1p+0,                -0x1p-53,   -0x1.0000000000001p+0,         -0x1p-1, INEXACT)
+T(RD,                 -0x1p+0,             -0x1.01p-53,   -0x1.0000000000001p+0,      -0x1.fep-2, INEXACT)
+T(RD,                 -0x1p+0,                 0x1p-54,                 -0x1p+0,         -0x1p-2, INEXACT)
+T(RD,                 -0x1p+0,              0x1.01p-54,                 -0x1p+0,      -0x1.01p-2, INEXACT)
+
+T(RZ,                  0x1p+0,                 0x1p-52,    0x1.0000000000001p+0,          0x0p+0, 0)
+T(RZ,                  0x1p+0,                 0x1p-53,                  0x1p+0,         -0x1p-1, INEXACT)
+T(RZ,                  0x1p+0,              0x1.01p-53,                  0x1p+0,      -0x1.01p-1, INEXACT)
+T(RZ,                  0x1p+0,                -0x1p-54,    0x1.fffffffffffffp-1,         -0x1p-1, INEXACT)
+T(RZ,                  0x1p+0,             -0x1.01p-54,    0x1.fffffffffffffp-1,      -0x1.fep-2, INEXACT)
+T(RZ,                 -0x1p+0,                -0x1p-53,                 -0x1p+0,          0x1p-1, INEXACT)
+T(RZ,                 -0x1p+0,             -0x1.01p-53,                 -0x1p+0,       0x1.01p-1, INEXACT)
+T(RZ,                 -0x1p+0,                 0x1p-54,   -0x1.fffffffffffffp-1,          0x1p-1, INEXACT)
+T(RZ,                 -0x1p+0,              0x1.01p-54,   -0x1.fffffffffffffp-1,       0x1.fep-2, INEXACT)
+};
+
+static void test_round_add(void)
+{
+       #pragma STDC FENV_ACCESS ON
+       double y;
+       float d;
+       int i;
+       struct dd_d *p;
+
+       for (i = 0; i < sizeof t/sizeof *t; i++) {
+               p = t + i;
+
+               if (p->r < 0)
+                       continue;
+               fesetround(p->r);
+               y = p->x + p->x2;
+               d = ulperr(y, p->y, p->dy);
+               if (!checkcr(y, p->y, p->r)) {
+                       printf("%s:%d: %s %a+%a want %a got %a ulperr %.3f = %a + %a\n",
+                               p->file, p->line, rstr(p->r), p->x, p->x2, p->y, y, d, d-p->dy, p->dy);
+                       test_status = 1;
+               }
+       }
+}
+
+int main(void)
+{
+       test_except();
+       test_round();
+       test_round_add();
+       return test_status;
+}

diff --git a/src/math/gen/Makefile b/src/math/gen/Makefile
index 8394478..5760c23 100644 (file)
--- a/src/math/gen/Makefile
+++ b/src/math/gen/Makefile
@@ -2,7 +2,7 @@ CFLAGS=-I. -Wall -fno-builtin -ffloat-store -D_GNU_SOURCE
 U=mpfr
 T=$(wildcard t*.c)
 
-all: gen check
+all: gen check mgen
 
 gen: gen.c util.c mp.c $(U)/lib/libmpfr.a $(U)/lib/libgmp.a
        musl-gcc -o $@ $(CFLAGS) -I$(U)/include $^
@@ -10,8 +10,8 @@ gen: gen.c util.c mp.c $(U)/lib/libmpfr.a $(U)/lib/libgmp.a
 check: gen.c util.c mplibm.c
        musl-gcc -o $@ $(CFLAGS) -lm $^
 
-#genlibm: gen.c util.c mplibm.c
-#      musl-gcc -o $@ $(CFLAGS) -lm $^
+mgen: gen.c util.c mplibm.c
+       musl-gcc -o $@ $(CFLAGS) -lm $^
 
 clean:
        rm -f gen check genlibm

diff --git a/src/math/gen/functions.h b/src/math/gen/functions.h
index f23fad1..c8f2756 100644 (file)
--- a/src/math/gen/functions.h
+++ b/src/math/gen/functions.h
@@ -1,3 +1,5 @@
+T(add,         dd_d)
+
 T(acos,        d_d)
 T(acosf,       f_f)
 T(acosl,       l_l)

diff --git a/src/math/gen/mp.c b/src/math/gen/mp.c
index 172b454..be2ed39 100644 (file)
--- a/src/math/gen/mp.c
+++ b/src/math/gen/mp.c
@@ -391,6 +391,8 @@ static int wrap_pow10(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
        return mpfr_ui_pow(my, 10, mx, r);
 }
 
+int mpadd(struct t *t) { return mpd2(t, mpfr_add); }
+
 int mpacos(struct t *t) { return mpd1(t, mpfr_acos); }
 int mpacosf(struct t *t) { return mpf1(t, mpfr_acos); }
 int mpacosl(struct t *t) { return mpl1(t, mpfr_acos); }

diff --git a/src/math/gen/mplibm.c b/src/math/gen/mplibm.c
index ccb7c81..04af049 100644 (file)
--- a/src/math/gen/mplibm.c
+++ b/src/math/gen/mplibm.c
@@ -54,6 +54,8 @@ static int mpl2(struct t *s, long double (*f)(long double, long double))
        return 0;
 }
 
+static double add(double x, double y) { double z = x + y; return z; }
+int mpadd(struct t *t) { return mpd2(t, add); }
 
 int mpacos(struct t *t) { return mpd1(t, acos); }
 int mpacosf(struct t *t) { return mpf1(t, acosf); }