From 21dd45c44f745d3b85afab304362d01b3b45a5e4 Mon Sep 17 00:00:00 2001 From: nsz Date: Sat, 13 Oct 2012 19:35:16 +0200 Subject: [PATCH] math: add jn and yn --- src/math/gen/functions.h | 4 ++++ src/math/gen/gentests.sh | 4 ++-- src/math/gen/mp.c | 16 ++++++++++++++++ src/math/gen/mplibm.c | 17 +++++++++++++++++ src/math/jn.c | 41 ++++++++++++++++++++++++++++++++++++++++ src/math/jnf.c | 41 ++++++++++++++++++++++++++++++++++++++++ src/math/sanity/jn.h | 10 ++++++++++ src/math/sanity/jnf.h | 10 ++++++++++ src/math/sanity/yn.h | 10 ++++++++++ src/math/sanity/ynf.h | 10 ++++++++++ src/math/yn.c | 41 ++++++++++++++++++++++++++++++++++++++++ src/math/ynf.c | 41 ++++++++++++++++++++++++++++++++++++++++ 12 files changed, 243 insertions(+), 2 deletions(-) create mode 100644 src/math/jn.c create mode 100644 src/math/jnf.c create mode 100644 src/math/sanity/jn.h create mode 100644 src/math/sanity/jnf.h create mode 100644 src/math/sanity/yn.h create mode 100644 src/math/sanity/ynf.h create mode 100644 src/math/yn.c create mode 100644 src/math/ynf.c diff --git a/src/math/gen/functions.h b/src/math/gen/functions.h index bbc4864..f23fad1 100644 --- a/src/math/gen/functions.h +++ b/src/math/gen/functions.h @@ -196,3 +196,7 @@ T(fma, ddd_d) T(fmaf, fff_f) T(fmal, lll_l) +T(jn, di_d) +T(jnf, fi_f) +T(yn, di_d) +T(ynf, fi_f) diff --git a/src/math/gen/gentests.sh b/src/math/gen/gentests.sh index 23eee22..b47ce7c 100755 --- a/src/math/gen/gentests.sh +++ b/src/math/gen/gentests.sh @@ -5,7 +5,7 @@ sed 's/^T(//;s/,//;s/)//' functions.h | while read N T do [ "$T" ] || continue -# [ -e $D/$N.c ] || { + [ -e $D/$N.c ] || { cp template/$T.c $D/$N.c || continue ND=`echo $N |sed 's/l$//'` @@ -21,5 +21,5 @@ do done sed -i "s/___/$N/g;s,DHEADERS,$DH,;s,HEADERS,$H," $D/$N.c -# } + } done diff --git a/src/math/gen/mp.c b/src/math/gen/mp.c index cb10c41..bebc019 100644 --- a/src/math/gen/mp.c +++ b/src/math/gen/mp.c @@ -355,6 +355,15 @@ static int wrap_remquo(mpfr_t my, const mpfr_t mx, const mpfr_t mx2, mpfr_rnd_t { return mpfr_remquo(my, &mpremquo_q, mx, mx2, r); } +static int mpbessel_n; +static int wrap_jn(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r) +{ + return mpfr_jn(my, mpbessel_n, mx, r); +} +static int wrap_yn(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r) +{ + return mpfr_yn(my, mpbessel_n, mx, r); +} static int wrap_ceil(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r) { return mpfr_ceil(my, mx); @@ -961,3 +970,10 @@ int mpfmal(struct t *t) #endif } +int mpjn(struct t *t) { mpbessel_n = t->i; return mpd1(t, wrap_jn); } +int mpjnf(struct t *t) { mpbessel_n = t->i; return mpf1(t, wrap_jn); } +int mpjnl(struct t *t) { mpbessel_n = t->i; return mpl1(t, wrap_jn); } +int mpyn(struct t *t) { mpbessel_n = t->i; return mpd1(t, wrap_yn); } +int mpynf(struct t *t) { mpbessel_n = t->i; return mpf1(t, wrap_yn); } +int mpynl(struct t *t) { mpbessel_n = t->i; return mpl1(t, wrap_yn); } + diff --git a/src/math/gen/mplibm.c b/src/math/gen/mplibm.c index cf17d6d..ccb7c81 100644 --- a/src/math/gen/mplibm.c +++ b/src/math/gen/mplibm.c @@ -411,3 +411,20 @@ mp_fff_f(fma) mp_fff_f(fmaf) mp_fff_f(fmal) +#define mp_if_f(n) \ +int mp##n(struct t *t) \ +{ \ + t->dy = 0; \ + setupfenv(t->r); \ + t->y = n(t->i, t->x); \ + t->e = getexcept(); \ + return 0; \ +} + +mp_if_f(jn) +mp_if_f(jnf) +//mp_if_f(jnl) +mp_if_f(yn) +mp_if_f(ynf) +//mp_if_f(ynl) + diff --git a/src/math/jn.c b/src/math/jn.c new file mode 100644 index 0000000..96658dd --- /dev/null +++ b/src/math/jn.c @@ -0,0 +1,41 @@ +#include +#include +#include "util.h" + +static struct di_d t[] = { +#include "sanity/jn.h" +}; + +int main(void) +{ + #pragma STDC FENV_ACCESS ON + double y; + float d; + int e, i, err = 0; + struct di_d *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 = jn(p->i, p->x); + e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW); + + if (!checkexcept(e, p->e, p->r)) { + printf("%s:%d: bad fp exception: %s jn(%a, %lld)=%a, want %s", + p->file, p->line, rstr(p->r), p->x, p->i, p->y, estr(p->e)); + printf(" got %s\n", estr(e)); + err++; + } + d = ulperr(y, p->y, p->dy); + if (!checkulp(d, p->r)) { + printf("%s:%d: %s jn(%a, %lld) want %a got %a, ulperr %.3f = %a + %a\n", + p->file, p->line, rstr(p->r), p->x, p->i, p->y, y, d, d-p->dy, p->dy); + err++; + } + } + return !!err; +} diff --git a/src/math/jnf.c b/src/math/jnf.c new file mode 100644 index 0000000..7b163d7 --- /dev/null +++ b/src/math/jnf.c @@ -0,0 +1,41 @@ +#include +#include +#include "util.h" + +static struct fi_f t[] = { +#include "sanity/jnf.h" +}; + +int main(void) +{ + #pragma STDC FENV_ACCESS ON + double y; + float d; + int e, i, err = 0; + struct fi_f *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 = jnf(p->i, p->x); + e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW); + + if (!checkexcept(e, p->e, p->r)) { + printf("%s:%d: bad fp exception: %s jnf(%a, %lld)=%a, want %s", + p->file, p->line, rstr(p->r), p->x, p->i, p->y, estr(p->e)); + printf(" got %s\n", estr(e)); + err++; + } + d = ulperrf(y, p->y, p->dy); + if (!checkulp(d, p->r)) { + printf("%s:%d: %s jnf(%a, %lld) want %a got %a, ulperr %.3f = %a + %a\n", + p->file, p->line, rstr(p->r), p->x, p->i, p->y, y, d, d-p->dy, p->dy); + err++; + } + } + return !!err; +} diff --git a/src/math/sanity/jn.h b/src/math/sanity/jn.h new file mode 100644 index 0000000..ff1ba79 --- /dev/null +++ b/src/math/sanity/jn.h @@ -0,0 +1,10 @@ +T(RN, -0x1.02239f3c6a8f1p+3, -2, -0x1.8637401cdd96bp-4, -0x1.b2ea52p-3, INEXACT) +T(RN, 0x1.161868e18bc67p+2, -1, 0x1.7d48aacc4b11fp-3, -0x1.e6896ap-2, INEXACT) +T(RN, -0x1.0c34b3e01e6e7p+3, 0, 0x1.2fd575f8ad8b4p-4, -0x1.7d645cp-3, INEXACT) +T(RN, -0x1.a206f0a19dcc4p+2, 1, 0x1.288dbd86fec93p-3, -0x1.933dep-2, INEXACT) +T(RN, 0x1.288bbb0d6a1e6p+3, 2, 0x1.904ebb8f3e76p-3, 0x1.db86aep-3, INEXACT) +T(RN, 0x1.52efd0cd80497p-1, 3, 0x1.815a0471e7b9fp-8, 0x1.40aa2p-3, INEXACT) +T(RN, -0x1.a05cc754481d1p-2, 4, 0x1.2816cfe1d5454p-14, 0x1.49e706p-2, INEXACT) +T(RN, 0x1.1f9ef934745cbp-1, 5, 0x1.e274364abf2d5p-17, 0x1.fbe5c8p-2, INEXACT) +T(RN, 0x1.8c5db097f7442p-1, 6, 0x1.32d8157822faep-18, -0x1.6bc402p-6, INEXACT) +T(RN, -0x1.5b86ea8118a0ep-1, 7, -0x1.b39a9fa627d2p-24, 0x1.308cbp-3, INEXACT) diff --git a/src/math/sanity/jnf.h b/src/math/sanity/jnf.h new file mode 100644 index 0000000..9db5d86 --- /dev/null +++ b/src/math/sanity/jnf.h @@ -0,0 +1,10 @@ +T(RN, -0x1.0223ap+3, -2, -0x1.863726p-4, -0x1.28885p-7, INEXACT) +T(RN, 0x1.161868p+2, -1, 0x1.7d48a2p-3, -0x1.cd365cp-6, INEXACT) +T(RN, -0x1.0c34b4p+3, 0, 0x1.2fd572p-4, 0x1.556ac8p-3, INEXACT) +T(RN, -0x1.a206fp+2, 1, 0x1.288dc4p-3, 0x1.598192p-2, INEXACT) +T(RN, 0x1.288bbcp+3, 2, 0x1.904ec6p-3, 0x1.f826c6p-3, INEXACT) +T(RN, 0x1.52efdp-1, 3, 0x1.815a02p-8, 0x1.f13c2cp-4, INEXACT) +T(RN, -0x1.a05cc8p-2, 4, 0x1.2816d2p-14, 0x1.bdd79p-4, INEXACT) +T(RN, 0x1.1f9efap-1, 5, 0x1.e2743cp-17, -0x1.d9c372p-2, INEXACT) +T(RN, 0x1.8c5dbp-1, 6, 0x1.32d812p-18, -0x1.76a6fp-2, INEXACT) +T(RN, -0x1.5b86eap-1, 7, -0x1.b39a9cp-24, -0x1.83c12p-2, INEXACT) diff --git a/src/math/sanity/yn.h b/src/math/sanity/yn.h new file mode 100644 index 0000000..6be798b --- /dev/null +++ b/src/math/sanity/yn.h @@ -0,0 +1,10 @@ +T(RN, -0x1.02239f3c6a8f1p+3, -2, nan, 0x0p+0, INVALID) +T(RN, 0x1.161868e18bc67p+2, -1, -0x1.5ab54680a5966p-2, -0x1.498536p-2, INEXACT) +T(RN, -0x1.0c34b3e01e6e7p+3, 0, nan, 0x0p+0, INVALID) +T(RN, -0x1.a206f0a19dcc4p+2, 1, nan, 0x0p+0, INVALID) +T(RN, 0x1.288bbb0d6a1e6p+3, 2, -0x1.6e6fa7c6fbbf7p-3, -0x1.d687b8p-3, INEXACT) +T(RN, 0x1.52efd0cd80497p-1, 3, -0x1.2935d021584c1p+4, 0x1.739fecp-2, INEXACT) +T(RN, -0x1.a05cc754481d1p-2, 4, nan, 0x0p+0, INVALID) +T(RN, 0x1.1f9ef934745cbp-1, 5, -0x1.1691a15d9cdfp+12, 0x1.cdf17cp-6, INEXACT) +T(RN, 0x1.8c5db097f7442p-1, 6, -0x1.6dbc18f9d5ddep+13, -0x1.c1596ap-3, INEXACT) +T(RN, -0x1.5b86ea8118a0ep-1, 7, nan, 0x0p+0, INVALID) diff --git a/src/math/sanity/ynf.h b/src/math/sanity/ynf.h new file mode 100644 index 0000000..afe055e --- /dev/null +++ b/src/math/sanity/ynf.h @@ -0,0 +1,10 @@ +T(RN, -0x1.0223ap+3, -2, nan, 0x0p+0, INVALID) +T(RN, 0x1.161868p+2, -1, -0x1.5ab54ap-2, -0x1.6996cep-3, INEXACT) +T(RN, -0x1.0c34b4p+3, 0, nan, 0x0p+0, INVALID) +T(RN, -0x1.a206fp+2, 1, nan, 0x0p+0, INVALID) +T(RN, 0x1.288bbcp+3, 2, -0x1.6e6f9cp-3, -0x1.b964fp-3, INEXACT) +T(RN, 0x1.52efdp-1, 3, -0x1.2935d2p+4, 0x1.47d13p-4, INEXACT) +T(RN, -0x1.a05cc8p-2, 4, nan, 0x0p+0, INVALID) +T(RN, 0x1.1f9efap-1, 5, -0x1.16919ep+12, -0x1.d10586p-3, INEXACT) +T(RN, 0x1.8c5dbp-1, 6, -0x1.6dbc1cp+13, 0x1.d597eep-4, INEXACT) +T(RN, -0x1.5b86eap-1, 7, nan, 0x0p+0, INVALID) diff --git a/src/math/yn.c b/src/math/yn.c new file mode 100644 index 0000000..eb7e5ce --- /dev/null +++ b/src/math/yn.c @@ -0,0 +1,41 @@ +#include +#include +#include "util.h" + +static struct di_d t[] = { +#include "sanity/yn.h" +}; + +int main(void) +{ + #pragma STDC FENV_ACCESS ON + double y; + float d; + int e, i, err = 0; + struct di_d *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 = yn(p->i, p->x); + e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW); + + if (!checkexcept(e, p->e, p->r)) { + printf("%s:%d: bad fp exception: %s yn(%a, %lld)=%a, want %s", + p->file, p->line, rstr(p->r), p->x, p->i, p->y, estr(p->e)); + printf(" got %s\n", estr(e)); + err++; + } + d = ulperr(y, p->y, p->dy); + if (!checkulp(d, p->r)) { + printf("%s:%d: %s yn(%a, %lld) want %a got %a, ulperr %.3f = %a + %a\n", + p->file, p->line, rstr(p->r), p->x, p->i, p->y, y, d, d-p->dy, p->dy); + err++; + } + } + return !!err; +} diff --git a/src/math/ynf.c b/src/math/ynf.c new file mode 100644 index 0000000..fbd0359 --- /dev/null +++ b/src/math/ynf.c @@ -0,0 +1,41 @@ +#include +#include +#include "util.h" + +static struct fi_f t[] = { +#include "sanity/ynf.h" +}; + +int main(void) +{ + #pragma STDC FENV_ACCESS ON + double y; + float d; + int e, i, err = 0; + struct fi_f *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 = ynf(p->i, p->x); + e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW); + + if (!checkexcept(e, p->e, p->r)) { + printf("%s:%d: bad fp exception: %s ynf(%a, %lld)=%a, want %s", + p->file, p->line, rstr(p->r), p->x, p->i, p->y, estr(p->e)); + printf(" got %s\n", estr(e)); + err++; + } + d = ulperrf(y, p->y, p->dy); + if (!checkulp(d, p->r)) { + printf("%s:%d: %s ynf(%a, %lld) want %a got %a, ulperr %.3f = %a + %a\n", + p->file, p->line, rstr(p->r), p->x, p->i, p->y, y, d, d-p->dy, p->dy); + err++; + } + } + return !!err; +} -- 2.20.1