From ae0f0fe09b7fc9d44d072c3fd08372991d852b1d Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Sun, 18 Nov 2012 22:33:27 +0100 Subject: [PATCH] math: add fenv test and new 'add' function in gen --- src/math/fenv.c | 173 +++++++++++++++++++++++++++++++++++++++ src/math/gen/Makefile | 6 +- src/math/gen/functions.h | 2 + src/math/gen/mp.c | 2 + src/math/gen/mplibm.c | 2 + 5 files changed, 182 insertions(+), 3 deletions(-) create mode 100644 src/math/fenv.c diff --git a/src/math/fenv.c b/src/math/fenv.c new file mode 100644 index 0000000..c221d7f --- /dev/null +++ b/src/math/fenv.c @@ -0,0 +1,173 @@ +#include +#include +#include +#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 --- 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 --- 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 --- 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 --- 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); } -- 2.20.1