From cb9f87af9a1f917facd1140603b24f1087729751 Mon Sep 17 00:00:00 2001 From: nsz Date: Sat, 13 Oct 2012 17:07:30 +0200 Subject: [PATCH] math: add modf and sincos --- src/math/gen/functions.h | 7 +++ src/math/gen/gensanity.sh | 2 + src/math/gen/gentests.sh | 2 + src/math/gen/mp.c | 90 +++++++++++++++++++++++++++++++++++++++ src/math/gen/mplibm.c | 67 +++++++++++++++++++++++++++++ src/math/modf.c | 42 ++++++++++++++++++ src/math/modff.c | 42 ++++++++++++++++++ src/math/modfl.c | 46 ++++++++++++++++++++ src/math/sanity/modf.h | 10 +++++ src/math/sanity/modff.h | 10 +++++ src/math/sanity/modfl.h | 10 +++++ src/math/sanity/sincos.h | 10 +++++ src/math/sanity/sincosf.h | 10 +++++ src/math/sanity/sincosl.h | 10 +++++ src/math/sincos.c | 43 +++++++++++++++++++ src/math/sincosf.c | 43 +++++++++++++++++++ src/math/sincosl.c | 47 ++++++++++++++++++++ src/math/util.h | 3 ++ 18 files changed, 494 insertions(+) create mode 100644 src/math/modf.c create mode 100644 src/math/modff.c create mode 100644 src/math/modfl.c create mode 100644 src/math/sanity/modf.h create mode 100644 src/math/sanity/modff.h create mode 100644 src/math/sanity/modfl.h create mode 100644 src/math/sanity/sincos.h create mode 100644 src/math/sanity/sincosf.h create mode 100644 src/math/sanity/sincosl.h create mode 100644 src/math/sincos.c create mode 100644 src/math/sincosf.c create mode 100644 src/math/sincosl.c diff --git a/src/math/gen/functions.h b/src/math/gen/functions.h index 0612aaf..8cff966 100644 --- a/src/math/gen/functions.h +++ b/src/math/gen/functions.h @@ -176,3 +176,10 @@ T(llroundl, l_i) T(lround, d_i) T(lroundf, f_i) T(lroundl, l_i) +T(modf, d_dd) +T(modff, f_ff) +T(modfl, l_ll) +T(sincos, d_dd) +T(sincosf, f_ff) +T(sincosl, l_ll) + diff --git a/src/math/gen/gensanity.sh b/src/math/gen/gensanity.sh index 949c4b6..31bb996 100755 --- a/src/math/gen/gensanity.sh +++ b/src/math/gen/gensanity.sh @@ -5,6 +5,8 @@ D=../sanity sed 's/^T(//;s/,//;s/)//' functions.h | while read N T do case $T in + '') continue ;; + d_*|f_*|l_*) ./gen $N >$D/$N.h <y2 = t->y; + t->dy2 = t->dy; + e = t->e; + r = mpd1(t, mpfr_frac); + t->e |= e; + return r; +} + +int mpmodff(struct t *t) +{ + int e, r; + + r = mpf1(t, wrap_trunc); + if (r) + return r; + t->y2 = t->y; + t->dy2 = t->dy; + e = t->e; + r = mpf1(t, mpfr_frac); + t->e |= e; + return r; +} + +int mpmodfl(struct t *t) +{ + int e, r; + + r = mpl1(t, wrap_trunc); + if (r) + return r; + t->y2 = t->y; + t->dy2 = t->dy; + e = t->e; + r = mpl1(t, mpfr_frac); + t->e |= e; + return r; +} + +int mpsincos(struct t *t) +{ + int e, r; + + r = mpd1(t, mpfr_cos); + if (r) + return r; + t->y2 = t->y; + t->dy2 = t->dy; + e = t->e; + r = mpd1(t, mpfr_sin); + t->e |= e; + return r; +} + +int mpsincosf(struct t *t) +{ + int e, r; + + r = mpf1(t, mpfr_cos); + if (r) + return r; + t->y2 = t->y; + t->dy2 = t->dy; + e = t->e; + r = mpf1(t, mpfr_sin); + t->e |= e; + return r; +} + +int mpsincosl(struct t *t) +{ + int e, r; + + r = mpl1(t, mpfr_cos); + if (r) + return r; + t->y2 = t->y; + t->dy2 = t->dy; + e = t->e; + r = mpl1(t, mpfr_sin); + t->e |= e; + return r; +} + diff --git a/src/math/gen/mplibm.c b/src/math/gen/mplibm.c index 93225a9..62b8721 100644 --- a/src/math/gen/mplibm.c +++ b/src/math/gen/mplibm.c @@ -313,3 +313,70 @@ mp_f_i(lround) mp_f_i(lroundf) mp_f_i(lroundl) +int mpmodf(struct t *t) +{ + double y2; + + t->dy = t->dy2 = 0; + setupfenv(t->r); + t->y = modf(t->x, &y2); + t->y2 = y2; + t->e = getexcept(); + return 0; +} + +int mpmodff(struct t *t) +{ + float y2; + + t->dy = t->dy2 = 0; + setupfenv(t->r); + t->y = modff(t->x, &y2); + t->y2 = y2; + t->e = getexcept(); + return 0; +} + +int mpmodfl(struct t *t) +{ + t->dy = t->dy2 = 0; + setupfenv(t->r); + t->y = modfl(t->x, &t->y2); + t->e = getexcept(); + return 0; +} + +int mpsincos(struct t *t) +{ + double y, y2; + + t->dy = t->dy2 = 0; + setupfenv(t->r); + sincos(t->x, &y, &y2); + t->y = y; + t->y2 = y2; + t->e = getexcept(); + return 0; +} + +int mpsincosf(struct t *t) +{ + float y, y2; + + t->dy = t->dy2 = 0; + setupfenv(t->r); + sincosf(t->x, &y, &y2); + t->y = y; + t->y2 = y2; + t->e = getexcept(); + return 0; +} + +int mpsincosl(struct t *t) +{ + t->dy = t->dy2 = 0; + setupfenv(t->r); + sincosl(t->x, &t->y, &t->y2); + t->e = getexcept(); + return 0; +} diff --git a/src/math/modf.c b/src/math/modf.c new file mode 100644 index 0000000..86be04f --- /dev/null +++ b/src/math/modf.c @@ -0,0 +1,42 @@ +#include +#include +#include "util.h" + +static struct d_dd t[] = { +#include "sanity/modf.h" +}; + +int main(void) +{ + #pragma STDC FENV_ACCESS ON + double y, yi; + float d, di; + int e, i, err = 0; + struct d_dd *p; + + for (i = 0; i < sizeof t/sizeof *t; i++) { + p = t + i; + + if (p->r < 0) + continue; + fesetround(p->r); + feclearexcept(FE_ALL_EXCEPT); + y = modf(p->x, &yi); + e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW); + + if (!checkexcept(e, p->e, p->r)) { + printf("%s:%d: bad fp exception: %s modf(%a)=%a,%a, want %s", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, estr(p->e)); + printf(" got %s\n", estr(e)); + err++; + } + d = ulperr(y, p->y, p->dy); + di = ulperr(yi, p->y2, p->dy2); + if (!checkulp(d, p->r) || !checkulp(di, p->r)) { + printf("%s:%d: %s modf(%a) want %a,%a got %a,%a, ulperr %.3f = %a + %a, %.3f = %a + %a\n", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, y, yi, d, d-p->dy, p->dy, di, di-p->dy2, p->dy2); + err++; + } + } + return !!err; +} diff --git a/src/math/modff.c b/src/math/modff.c new file mode 100644 index 0000000..e6ed290 --- /dev/null +++ b/src/math/modff.c @@ -0,0 +1,42 @@ +#include +#include +#include "util.h" + +static struct f_ff t[] = { +#include "sanity/modff.h" +}; + +int main(void) +{ + #pragma STDC FENV_ACCESS ON + float y, yi; + float d, di; + int e, i, err = 0; + struct f_ff *p; + + for (i = 0; i < sizeof t/sizeof *t; i++) { + p = t + i; + + if (p->r < 0) + continue; + fesetround(p->r); + feclearexcept(FE_ALL_EXCEPT); + y = modff(p->x, &yi); + e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW); + + if (!checkexcept(e, p->e, p->r)) { + printf("%s:%d: bad fp exception: %s modff(%a)=%a,%a, want %s", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, estr(p->e)); + printf(" got %s\n", estr(e)); + err++; + } + d = ulperr(y, p->y, p->dy); + di = ulperr(yi, p->y2, p->dy2); + if (!checkulp(d, p->r) || !checkulp(di, p->r)) { + printf("%s:%d: %s modff(%a) want %a,%a got %a,%a, ulperr %.3f = %a + %a, %.3f = %a + %a\n", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, y, yi, d, d-p->dy, p->dy, di, di-p->dy2, p->dy2); + err++; + } + } + return !!err; +} diff --git a/src/math/modfl.c b/src/math/modfl.c new file mode 100644 index 0000000..72d3b1d --- /dev/null +++ b/src/math/modfl.c @@ -0,0 +1,46 @@ +#include +#include +#include "util.h" + +static struct l_ll t[] = { +#if LDBL_MANT_DIG == 53 +#include "sanity/modf.h" +#elif LDBL_MANT_DIG == 64 +#include "sanity/modfl.h" +#endif +}; + +int main(void) +{ + #pragma STDC FENV_ACCESS ON + long double y, yi; + float d, di; + int e, i, err = 0; + struct l_ll *p; + + for (i = 0; i < sizeof t/sizeof *t; i++) { + p = t + i; + + if (p->r < 0) + continue; + fesetround(p->r); + feclearexcept(FE_ALL_EXCEPT); + y = modfl(p->x, &yi); + e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW); + + if (!checkexcept(e, p->e, p->r)) { + printf("%s:%d: bad fp exception: %s modf(%La)=%La,%La, want %s", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, estr(p->e)); + printf(" got %s\n", estr(e)); + err++; + } + d = ulperr(y, p->y, p->dy); + di = ulperr(yi, p->y2, p->dy2); + if (!checkulp(d, p->r) || !checkulp(di, p->r)) { + printf("%s:%d: %s modf(%La) want %La,%La got %La,%La, ulperr %.3f = %a + %a, %.3f = %a + %a\n", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, y, yi, d, d-p->dy, p->dy, di, di-p->dy2, p->dy2); + err++; + } + } + return !!err; +} diff --git a/src/math/sanity/modf.h b/src/math/sanity/modf.h new file mode 100644 index 0000000..c8fdd24 --- /dev/null +++ b/src/math/sanity/modf.h @@ -0,0 +1,10 @@ +T(RN, -0x1.02239f3c6a8f1p+3, -0x1.11cf9e354788p-4, 0x0p+0, -0x1p+3, 0x0p+0, INEXACT) +T(RN, 0x1.161868e18bc67p+2, 0x1.61868e18bc67p-2, 0x0p+0, 0x1p+2, 0x0p+0, INEXACT) +T(RN, -0x1.0c34b3e01e6e7p+3, -0x1.86967c03cdcep-2, 0x0p+0, -0x1p+3, 0x0p+0, INEXACT) +T(RN, -0x1.a206f0a19dcc4p+2, -0x1.1037850cee62p-1, 0x0p+0, -0x1.8p+2, 0x0p+0, INEXACT) +T(RN, 0x1.288bbb0d6a1e6p+3, 0x1.117761ad43ccp-2, 0x0p+0, 0x1.2p+3, 0x0p+0, INEXACT) +T(RN, 0x1.52efd0cd80497p-1, 0x1.52efd0cd80497p-1, 0x0p+0, 0x0p+0, 0x0p+0, INEXACT) +T(RN, -0x1.a05cc754481d1p-2, -0x1.a05cc754481d1p-2, 0x0p+0, -0x0p+0, 0x0p+0, INEXACT) +T(RN, 0x1.1f9ef934745cbp-1, 0x1.1f9ef934745cbp-1, 0x0p+0, 0x0p+0, 0x0p+0, INEXACT) +T(RN, 0x1.8c5db097f7442p-1, 0x1.8c5db097f7442p-1, 0x0p+0, 0x0p+0, 0x0p+0, INEXACT) +T(RN, -0x1.5b86ea8118a0ep-1, -0x1.5b86ea8118a0ep-1, 0x0p+0, -0x0p+0, 0x0p+0, INEXACT) diff --git a/src/math/sanity/modff.h b/src/math/sanity/modff.h new file mode 100644 index 0000000..3b9f798 --- /dev/null +++ b/src/math/sanity/modff.h @@ -0,0 +1,10 @@ +T(RN, -0x1.0223ap+3, -0x1.11dp-4, 0x0p+0, -0x1p+3, 0x0p+0, INEXACT) +T(RN, 0x1.161868p+2, 0x1.61868p-2, 0x0p+0, 0x1p+2, 0x0p+0, INEXACT) +T(RN, -0x1.0c34b4p+3, -0x1.86968p-2, 0x0p+0, -0x1p+3, 0x0p+0, INEXACT) +T(RN, -0x1.a206fp+2, -0x1.10378p-1, 0x0p+0, -0x1.8p+2, 0x0p+0, INEXACT) +T(RN, 0x1.288bbcp+3, 0x1.11778p-2, 0x0p+0, 0x1.2p+3, 0x0p+0, INEXACT) +T(RN, 0x1.52efdp-1, 0x1.52efdp-1, 0x0p+0, 0x0p+0, 0x0p+0, INEXACT) +T(RN, -0x1.a05cc8p-2, -0x1.a05cc8p-2, 0x0p+0, -0x0p+0, 0x0p+0, INEXACT) +T(RN, 0x1.1f9efap-1, 0x1.1f9efap-1, 0x0p+0, 0x0p+0, 0x0p+0, INEXACT) +T(RN, 0x1.8c5dbp-1, 0x1.8c5dbp-1, 0x0p+0, 0x0p+0, 0x0p+0, INEXACT) +T(RN, -0x1.5b86eap-1, -0x1.5b86eap-1, 0x0p+0, -0x0p+0, 0x0p+0, INEXACT) diff --git a/src/math/sanity/modfl.h b/src/math/sanity/modfl.h new file mode 100644 index 0000000..aa42f90 --- /dev/null +++ b/src/math/sanity/modfl.h @@ -0,0 +1,10 @@ +T(RN, -0x1.02239f3c6a8f13dep+3L, -0x1.11cf9e354789efp-4L, 0x0p+0, -0x1p+3L, 0x0p+0, INEXACT) +T(RN, 0x1.161868e18bc67782p+2L, 0x1.61868e18bc67782p-2L, 0x0p+0, 0x1p+2L, 0x0p+0, INEXACT) +T(RN, -0x1.0c34b3e01e6e682cp+3L, -0x1.86967c03cdcd058p-2L, 0x0p+0, -0x1p+3L, 0x0p+0, INEXACT) +T(RN, -0x1.a206f0a19dcc3948p+2L, -0x1.1037850cee61ca4p-1L, 0x0p+0, -0x1.8p+2L, 0x0p+0, INEXACT) +T(RN, 0x1.288bbb0d6a1e5bdap+3L, 0x1.117761ad43cb7b4p-2L, 0x0p+0, 0x1.2p+3L, 0x0p+0, INEXACT) +T(RN, 0x1.52efd0cd80496a5ap-1L, 0x1.52efd0cd80496a5ap-1L, 0x0p+0, 0x0p+0L, 0x0p+0, INEXACT) +T(RN, -0x1.a05cc754481d0bdp-2L, -0x1.a05cc754481d0bdp-2L, 0x0p+0, -0x0p+0L, 0x0p+0, INEXACT) +T(RN, 0x1.1f9ef934745cad6p-1L, 0x1.1f9ef934745cad6p-1L, 0x0p+0, 0x0p+0L, 0x0p+0, INEXACT) +T(RN, 0x1.8c5db097f744257ep-1L, 0x1.8c5db097f744257ep-1L, 0x0p+0, 0x0p+0L, 0x0p+0, INEXACT) +T(RN, -0x1.5b86ea8118a0e2bcp-1L, -0x1.5b86ea8118a0e2bcp-1L, 0x0p+0, -0x0p+0L, 0x0p+0, INEXACT) diff --git a/src/math/sanity/sincos.h b/src/math/sanity/sincos.h new file mode 100644 index 0000000..283e961 --- /dev/null +++ b/src/math/sanity/sincos.h @@ -0,0 +1,10 @@ +T(RN, -0x1.02239f3c6a8f1p+3, -0x1.f4719cbe20bd2p-1, -0x1.2a4a16p-3, -0x1.b0aa8f2c9baf6p-3, -0x1.c105d2p-4, INEXACT) +T(RN, 0x1.161868e18bc67p+2, -0x1.dde0a33834424p-1, -0x1.6902d6p-4, -0x1.6f922aed88704p-2, -0x1.b8b8fap-4, INEXACT) +T(RN, -0x1.0c34b3e01e6e7p+3, -0x1.ba6a5410cb9ccp-1, -0x1.e1078ap-4, -0x1.01b4e00041423p-1, -0x1.5f1decp-6, INEXACT) +T(RN, -0x1.a206f0a19dcc4p+2, -0x1.f7aed6ca5f32fp-3, -0x1.040d5p-3, 0x1.f0462a6686a9cp-1, -0x1.ea474ep-2, INEXACT) +T(RN, 0x1.288bbb0d6a1e6p+3, 0x1.41acd05fae3c4p-3, -0x1.e4265ap-6, -0x1.f9a51be5829b7p-1, 0x1.f3c7cep-2, INEXACT) +T(RN, 0x1.52efd0cd80497p-1, 0x1.3ab7ecc98df9ap-1, -0x1.98a5aep-4, 0x1.93da10e89d4d1p-1, 0x1.044604p-3, INEXACT) +T(RN, -0x1.a05cc754481d1p-2, -0x1.94fbf72645bfcp-2, -0x1.77aebcp-2, 0x1.d64199a5cb117p-1, -0x1.0b79e2p-2, INEXACT) +T(RN, 0x1.1f9ef934745cbp-1, 0x1.10baf3a5f550ep-1, -0x1.6b8fcep-2, 0x1.b150bae7795b1p-1, -0x1.35d926p-2, INEXACT) +T(RN, 0x1.8c5db097f7442p-1, 0x1.65f1c5e591db2p-1, -0x1.b5efc2p-2, 0x1.6e164e427022bp-1, -0x1.5db4c2p-4, INEXACT) +T(RN, -0x1.5b86ea8118a0ep-1, -0x1.417318671b83dp-1, -0x1.87ffcp-2, 0x1.8e83d35a366cp-1, 0x1.3c524p-2, INEXACT) diff --git a/src/math/sanity/sincosf.h b/src/math/sanity/sincosf.h new file mode 100644 index 0000000..757317d --- /dev/null +++ b/src/math/sanity/sincosf.h @@ -0,0 +1,10 @@ +T(RN, -0x1.0223ap+3, -0x1.f4719ap-1, 0x1.481cf2p-4, -0x1.b0aabep-3, 0x1.eee272p-2, INEXACT) +T(RN, 0x1.161868p+2, -0x1.dde0ap-1, 0x1.6107cap-2, -0x1.6f9238p-2, 0x1.5c33e2p-5, INEXACT) +T(RN, -0x1.0c34b4p+3, -0x1.ba6a54p-1, -0x1.dfe862p-2, -0x1.01b4e2p-1, -0x1.1be494p-3, INEXACT) +T(RN, -0x1.a206fp+2, -0x1.f7aec4p-3, -0x1.95029cp-2, 0x1.f0462cp-1, 0x1.6df7bcp-3, INEXACT) +T(RN, 0x1.288bbcp+3, 0x1.41ac94p-3, -0x1.eba8d2p-3, -0x1.f9a51ep-1, 0x1.1c971cp-3, INEXACT) +T(RN, 0x1.52efdp-1, 0x1.3ab7ecp-1, -0x1.3bafcap-4, 0x1.93da12p-1, 0x1.322268p-2, INEXACT) +T(RN, -0x1.a05cc8p-2, -0x1.94fbf8p-2, -0x1.e01394p-4, 0x1.d6419ap-1, 0x1.f0a754p-3, INEXACT) +T(RN, 0x1.1f9efap-1, 0x1.10baf4p-1, -0x1.48e402p-3, 0x1.b150bap-1, -0x1.ec3366p-3, INEXACT) +T(RN, 0x1.8c5dbp-1, 0x1.65f1c6p-1, 0x1.0e2d0ap-2, 0x1.6e164ep-1, -0x1.595b9cp-2, INEXACT) +T(RN, -0x1.5b86eap-1, -0x1.417318p-1, 0x1.5010ccp-8, 0x1.8e83d4p-1, 0x1.52f278p-3, INEXACT) diff --git a/src/math/sanity/sincosl.h b/src/math/sanity/sincosl.h new file mode 100644 index 0000000..5abadc1 --- /dev/null +++ b/src/math/sanity/sincosl.h @@ -0,0 +1,10 @@ +T(RN, -0x1.02239f3c6a8f13dep+3L, -0x1.f4719cbe20bd109ap-1L, -0x1.f6addp-2, -0x1.b0aa8f2c9bb05028p-3L, 0x1.cb66d8p-2, INEXACT) +T(RN, 0x1.161868e18bc67782p+2L, -0x1.dde0a33834425426p-1L, 0x1.2fdc6p-3, -0x1.6f922aed886fce28p-2L, 0x1.857366p-2, INEXACT) +T(RN, -0x1.0c34b3e01e6e682cp+3L, -0x1.ba6a5410cb9cfd2ap-1L, -0x1.41edcap-4, -0x1.01b4e0004141c36ep-1L, -0x1.496d12p-6, INEXACT) +T(RN, -0x1.a206f0a19dcc3948p+2L, -0x1.f7aed6ca5f321d92p-3L, 0x1.ba498cp-3, 0x1.f0462a6686a9d4e2p-1L, 0x1.882674p-2, INEXACT) +T(RN, 0x1.288bbb0d6a1e5bdap+3L, 0x1.41acd05fae3d46aep-3L, 0x1.30647p-2, -0x1.f9a51be5829b6d62p-1L, 0x1.cb35d4p-4, INEXACT) +T(RN, 0x1.52efd0cd80496a5ap-1L, 0x1.3ab7ecc98df99d24p-1L, -0x1.4d10a4p-5, 0x1.93da10e89d4d117p-1L, -0x1.25ac1p-3, INEXACT) +T(RN, -0x1.a05cc754481d0bdp-2L, -0x1.94fbf72645bfb648p-2L, 0x1.5772aep-2, 0x1.d64199a5cb117502p-1L, 0x1.c94d6p-5, INEXACT) +T(RN, 0x1.1f9ef934745cad6p-1L, 0x1.10baf3a5f550e376p-1L, 0x1.eaff04p-3, 0x1.b150bae7795b163ep-1L, 0x1.4c1268p-2, INEXACT) +T(RN, 0x1.8c5db097f744257ep-1L, 0x1.65f1c5e591db2ac6p-1L, 0x1.e36f2ep-2, 0x1.6e164e427022ad86p-1L, -0x1.83a52p-2, INEXACT) +T(RN, -0x1.5b86ea8118a0e2bcp-1L, -0x1.417318671b83ccp-1L, 0x1.b27416p-2, 0x1.8e83d35a366bf958p-1L, 0x1.88bbdp-2, INEXACT) diff --git a/src/math/sincos.c b/src/math/sincos.c new file mode 100644 index 0000000..3e574f3 --- /dev/null +++ b/src/math/sincos.c @@ -0,0 +1,43 @@ +#include +#include +#include "util.h" + +static struct d_dd t[] = { +#include "sanity/sincos.h" +}; + +int main(void) +{ + #pragma STDC FENV_ACCESS ON + double ysin, ycos; + float dsin, dcos; + int e, i, err = 0; + struct d_dd *p; + + for (i = 0; i < sizeof t/sizeof *t; i++) { + p = t + i; + + if (p->r < 0) + continue; + fesetround(p->r); + feclearexcept(FE_ALL_EXCEPT); + sincos(p->x, &ysin, &ycos); + e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW); + + if (!checkexcept(e, p->e, p->r)) { + printf("%s:%d: bad fp exception: %s sincos(%a)=%a,%a, want %s", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, estr(p->e)); + printf(" got %s\n", estr(e)); + err++; + } + dsin = ulperr(ysin, p->y, p->dy); + dcos = ulperr(ycos, p->y2, p->dy2); + if (!checkulp(dsin, p->r) || !checkulp(dcos, p->r)) { + printf("%s:%d: %s sincos(%a) want %a,%a got %a,%a, ulperr %.3f = %a + %a, %.3f = %a + %a\n", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, ysin, ycos, + dsin, dsin-p->dy, p->dy, dcos, dcos-p->dy2, p->dy2); + err++; + } + } + return !!err; +} diff --git a/src/math/sincosf.c b/src/math/sincosf.c new file mode 100644 index 0000000..e21fbd9 --- /dev/null +++ b/src/math/sincosf.c @@ -0,0 +1,43 @@ +#include +#include +#include "util.h" + +static struct f_ff t[] = { +#include "sanity/sincosf.h" +}; + +int main(void) +{ + #pragma STDC FENV_ACCESS ON + float ysin, ycos; + float dsin, dcos; + int e, i, err = 0; + struct f_ff *p; + + for (i = 0; i < sizeof t/sizeof *t; i++) { + p = t + i; + + if (p->r < 0) + continue; + fesetround(p->r); + feclearexcept(FE_ALL_EXCEPT); + sincosf(p->x, &ysin, &ycos); + e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW); + + if (!checkexcept(e, p->e, p->r)) { + printf("%s:%d: bad fp exception: %s sincosf(%a)=%a,%a, want %s", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, estr(p->e)); + printf(" got %s\n", estr(e)); + err++; + } + dsin = ulperr(ysin, p->y, p->dy); + dcos = ulperr(ycos, p->y2, p->dy2); + if (!checkulp(dsin, p->r) || !checkulp(dcos, p->r)) { + printf("%s:%d: %s sincosf(%a) want %a,%a got %a,%a, ulperr %.3f = %a + %a, %.3f = %a + %a\n", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, ysin, ycos, + dsin, dsin-p->dy, p->dy, dcos, dcos-p->dy2, p->dy2); + err++; + } + } + return !!err; +} diff --git a/src/math/sincosl.c b/src/math/sincosl.c new file mode 100644 index 0000000..697eeda --- /dev/null +++ b/src/math/sincosl.c @@ -0,0 +1,47 @@ +#include +#include +#include "util.h" + +static struct l_ll t[] = { +#if LDBL_MANT_DIG == 53 +#include "sanity/sincos.h" +#elif LDBL_MANT_DIG == 64 +#include "sanity/sincosl.h" +#endif +}; + +int main(void) +{ + #pragma STDC FENV_ACCESS ON + long double ysin, ycos; + float dsin, dcos; + int e, i, err = 0; + struct l_ll *p; + + for (i = 0; i < sizeof t/sizeof *t; i++) { + p = t + i; + + if (p->r < 0) + continue; + fesetround(p->r); + feclearexcept(FE_ALL_EXCEPT); + sincosl(p->x, &ysin, &ycos); + e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW); + + if (!checkexcept(e, p->e, p->r)) { + printf("%s:%d: bad fp exception: %s sincosl(%La)=%La,%La, want %s", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, estr(p->e)); + printf(" got %s\n", estr(e)); + err++; + } + dsin = ulperr(ysin, p->y, p->dy); + dcos = ulperr(ycos, p->y2, p->dy2); + if (!checkulp(dsin, p->r) || !checkulp(dcos, p->r)) { + printf("%s:%d: %s sincosl(%La) want %La,%La got %La,%La, ulperr %.3f = %a + %a, %.3f = %a + %a\n", + p->file, p->line, rstr(p->r), p->x, p->y, p->y2, ysin, ycos, + dsin, dsin-p->dy, p->dy, dcos, dcos-p->dy2, p->dy2); + err++; + } + } + return !!err; +} diff --git a/src/math/util.h b/src/math/util.h index 50b923a..1fc6939 100644 --- a/src/math/util.h +++ b/src/math/util.h @@ -81,6 +81,9 @@ struct li_l {POS int r; long double x; long long i; long double y; float dy; int struct d_i {POS int r; double x; long long i; int e; }; struct f_i {POS int r; float x; long long i; int e; }; struct l_i {POS int r; long double x; long long i; int e; }; +struct d_dd {POS int r; double x; double y; float dy; double y2; float dy2; int e; }; +struct f_ff {POS int r; float x; float y; float dy; float y2; float dy2; int e; }; +struct l_ll {POS int r; long double x; long double y; float dy; long double y2; float dy2; int e; }; #undef POS char *estr(int); -- 2.20.1