math: add jn and yn

diff --git a/src/math/gen/functions.h b/src/math/gen/functions.h
index bbc4864..f23fad1 100644 (file)
--- 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 (executable)
--- 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 (file)
--- 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 (file)
--- 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 (file)
index 0000000..96658dd
--- /dev/null
+++ b/src/math/jn.c
@@ -0,0 +1,41 @@
+#include <stdint.h>
+#include <stdio.h>
+#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 (file)
index 0000000..7b163d7
--- /dev/null
+++ b/src/math/jnf.c
@@ -0,0 +1,41 @@
+#include <stdint.h>
+#include <stdio.h>
+#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 (file)
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 (file)
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 (file)
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 (file)
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 (file)
index 0000000..eb7e5ce
--- /dev/null
+++ b/src/math/yn.c
@@ -0,0 +1,41 @@
+#include <stdint.h>
+#include <stdio.h>
+#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 (file)
index 0000000..fbd0359
--- /dev/null
+++ b/src/math/ynf.c
@@ -0,0 +1,41 @@
+#include <stdint.h>
+#include <stdio.h>
+#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;
+}