From f33814144a7fc8f79659f7beb7873bc24ea62b2a Mon Sep 17 00:00:00 2001 From: jtlap Date: Mon, 4 Dec 2023 20:16:27 +0100 Subject: [PATCH] sph h i and k --- include/kyosu/functions.hpp | 32 +- include/kyosu/types/impl/bessel.hpp | 2 + include/kyosu/types/impl/bessel/sb_hn.hpp | 162 +++++++++ include/kyosu/types/impl/bessel/sb_ikn.hpp | 331 ++++++++++++++++++ ...h_bessel_h0.inactive => sph_bessel_h0.cpp} | 0 ...h_bessel_h1.inactive => sph_bessel_h1.cpp} | 0 ...h_bessel_hn.inactive => sph_bessel_hn.cpp} | 0 ...h_bessel_i0.inactive => sph_bessel_i0.cpp} | 0 ...h_bessel_i1.inactive => sph_bessel_i1.cpp} | 0 test/unit/function/sph_bessel_in.cpp | 95 +++++ test/unit/function/sph_bessel_in.inactive | 76 ---- ...h_bessel_k0.inactive => sph_bessel_k0.cpp} | 0 ...h_bessel_k1.inactive => sph_bessel_k1.cpp} | 0 ...h_bessel_kn.inactive => sph_bessel_kn.cpp} | 0 14 files changed, 606 insertions(+), 92 deletions(-) create mode 100644 include/kyosu/types/impl/bessel/sb_hn.hpp create mode 100644 include/kyosu/types/impl/bessel/sb_ikn.hpp rename test/unit/function/{sph_bessel_h0.inactive => sph_bessel_h0.cpp} (100%) rename test/unit/function/{sph_bessel_h1.inactive => sph_bessel_h1.cpp} (100%) rename test/unit/function/{sph_bessel_hn.inactive => sph_bessel_hn.cpp} (100%) rename test/unit/function/{sph_bessel_i0.inactive => sph_bessel_i0.cpp} (100%) rename test/unit/function/{sph_bessel_i1.inactive => sph_bessel_i1.cpp} (100%) create mode 100644 test/unit/function/sph_bessel_in.cpp delete mode 100644 test/unit/function/sph_bessel_in.inactive rename test/unit/function/{sph_bessel_k0.inactive => sph_bessel_k0.cpp} (100%) rename test/unit/function/{sph_bessel_k1.inactive => sph_bessel_k1.cpp} (100%) rename test/unit/function/{sph_bessel_kn.inactive => sph_bessel_kn.cpp} (100%) diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index ca98e931..719ba164 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -184,13 +184,6 @@ #include #include -// #include -// #include -// #include - -// #include -// #include -// #include #include #include @@ -200,17 +193,24 @@ #include #include -// #include -// #include -// #include +#include +#include +#include + +#include +#include +#include -// #include -// #include -// #include +#include +#include +#include +#include +#include +#include +#include +#include +#include -// #include -// #include -// #include #include #include diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp index 46f8a4a3..5feafaf5 100644 --- a/include/kyosu/types/impl/bessel.hpp +++ b/include/kyosu/types/impl/bessel.hpp @@ -42,7 +42,9 @@ #include #include #include +#include #include +#include // #include diff --git a/include/kyosu/types/impl/bessel/sb_hn.hpp b/include/kyosu/types/impl/bessel/sb_hn.hpp new file mode 100644 index 00000000..1f29dd39 --- /dev/null +++ b/include/kyosu/types/impl/bessel/sb_hn.hpp @@ -0,0 +1,162 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +namespace kyosu::_ +{ + ///////////////////////////////// + // contains implementations of + // sph_bessel_h1_0 + // sph_bessel_h1_1 + // sph_bessel_h2_0 + // sph_bessel_h2_1 + // sph_bessel_h1n + // sph_bessel_h2n + // sph_bessel_h12n + ///////////////////////////////// + // utilities + // sb_h1_0 + // sb_h1_1 + // sb_h2_0 + // sb_h2_1 + // sb_h1n + // sb_h2n + // sb_h12n + //////////////////////////////// + + + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + // spherical bessel of the third kind + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + + template KYOSU_FORCEINLINE + auto sb_h2n(N n, Z z) noexcept + { + auto imzlt0 = eve::is_ltz(imag(z)); + z = if_else(imzlt0, z, -z); + using u_t = eve::underlying_type_t; + auto i = complex(u_t(0), u_t(1)); + auto rz = rec(z); + auto miz = mulmi(z); + auto h0 = if_else(imzlt0, i*exp(miz), -i*exp(-miz))*rz; + if(n == 0) return h0; + auto h1 = if_else(imzlt0, (rz+i),(rz-i))*h0 ; + if(n == 1) return if_else(imzlt0,h1, -h1); + + auto h2 = h1; + for(int i=1; i < n; ++i) + { + auto f = (2*i+1)*rz; + h2 = f*h1-h0; + h0 = h1; + h1 = h2; + } + return if_else(imzlt0, h2, eve::sign_alternate(u_t(n))*h2); + } + + + template KYOSU_FORCEINLINE + auto sb_h1n(N n, Z z) noexcept + { + using u_t = eve::underlying_type_t; + return eve::sign_alternate(u_t(n))*sph_bessel_h2n(n, -z); + } + + template KYOSU_FORCEINLINE + auto sb_h2_0(Z z) noexcept + { + return sb_h2n(0, z); + } + + template KYOSU_FORCEINLINE + auto sb_h2_1(Z z) noexcept + { + return sb_h2n(1, z); + } + + template KYOSU_FORCEINLINE + auto sb_h1n(int n, Z z) noexcept + { + using u_t = eve::underlying_type_t; + return eve::sign_alternate(u_t(n))*sph_bessel_h2n(n, -z); + } + + template KYOSU_FORCEINLINE + auto sb_h1_0(Z z) noexcept + { + return sb_h2_0(-z); + } + + template KYOSU_FORCEINLINE + auto sb_h1_1(Z z) noexcept + { + return -sb_h2_1(-z); + } + + //===------------------------------------------------------------------------------------------- + // sph_bessel_h1n + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, int n, Z z) noexcept + { + if constexpr(concepts::complex) + { + EVE_ASSERT(n >= 0, "spherical functions need positive order"); + return sb_h1n(n, z); + } + else + { + return cayley_extend_rev(sph_bessel_h1n, n, z); + } + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return sb_h1_0(z); + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return sb_h1_1(z); + } + + //===------------------------------------------------------------------------------------------- + // sph_bessel_h2n + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, N n, Z z) noexcept + { + EVE_ASSERT(n >= 0, "spherical functions need positive order"); + if constexpr(concepts::complex) + { + return sb_h2n(n, z); + } + else + { + return cayley_extend_rev(sph_bessel_h2n, n, z); + } + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return sb_h2n(0, z); + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return sb_h2n(1, z); + } + + +} diff --git a/include/kyosu/types/impl/bessel/sb_ikn.hpp b/include/kyosu/types/impl/bessel/sb_ikn.hpp new file mode 100644 index 00000000..f09ee5e4 --- /dev/null +++ b/include/kyosu/types/impl/bessel/sb_ikn.hpp @@ -0,0 +1,331 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +namespace kyosu::_ +{ + + ///////////////////////////////// + // contains implementations of + // sph_bessel_i0 + // sph_bessel_i1 + // sph_bessel_k0 + // sph_bessel_k1 + // sph_bessel_ikn + ///////////////////////////////// + // utilities + // sb_i1_0 + // sb_i1_1 + // sb_i2_0 + // sb_i2_1 + // sb_k0 + // sb_k1 + // sb_ikn + //////////////////////////////// + + template KYOSU_FORCEINLINE + auto sb_i1_0(Z z) noexcept + { + return sinhc(z); + } + + template KYOSU_FORCEINLINE + auto sb_i1_1(Z z) noexcept + { + auto [sh, ch] = sinhcosh(z); + auto rz = rec(z); + return fnma(sh, rz, ch)*rz; + } + + template KYOSU_FORCEINLINE + auto sb_i2n(N n, Z z) noexcept + { + using e_t = as_real_type_t; + auto iton = [](N n){ + if (n%4 == 0) return complex(eve::one(eve::as())); + else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::one(eve::as())); + else if (n%4 == 2) return complex(eve::mone(eve::as())); + else return complex(eve::zero(eve::as()), eve::mone(eve::as())); + }; + return rec(iton(n+1))*sph_bessel_yn(n,muli(z)); + } + + template KYOSU_FORCEINLINE + auto sb_i2n(N n, Z z, R & sis) noexcept + { + using e_t = as_real_type_t; + auto iton = [](N n){ + if (n%4 == 0) return complex(eve::one(eve::as())); + else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::one(eve::as())); + else if (n%4 == 2) return complex(eve::mone(eve::as())); + else return complex(eve::zero(eve::as()), eve::mone(eve::as())); + }; + return rec(iton(n +1))*sph_bessel_yn(n,muli(z), sis); + } + + + template KYOSU_FORCEINLINE + auto sb_i2_0(Z z) noexcept + { + return cosh(z)/z; + } + + template KYOSU_FORCEINLINE + auto sb_i2_1(Z z) noexcept + { + auto [sh, ch] = sinhcosh(z); + auto rz = rec(z); + return eve::fnma(ch, rz, sh)*rz; + } + + template KYOSU_FORCEINLINE + auto sb_i1n(N n, Z z) noexcept + { + if constexpr(concepts::complex ) + { + using e_t = as_real_type_t; + auto iton = [](N n){ + if (n%4 == 0) return complex(eve::one(eve::as())); + else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::one(eve::as())); + else if (n%4 == 2) return complex(eve::mone(eve::as())); + else return complex(eve::zero(eve::as()), eve::mone(eve::as())); + }; + return rec(iton(n))*sph_bessel_jn(n,muli(z)); + } + } + + template KYOSU_FORCEINLINE + auto sb_i1n(N n, Z z, R & sys) noexcept + { + using e_t = as_real_type_t; + auto iton = [](N n){ + if (n%4 == 0) return complex(eve::one(eve::as())); + else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::one(eve::as())); + else if (n%4 == 2) return complex(eve::mone(eve::as())); + else return complex(eve::zero(eve::as()), eve::mone(eve::as())); + }; + return rec(iton(n))*sph_bessel_yn(n,muli(z), sys); + } + + + + //===------------------------------------------------------------------------------------------- + // sb_kn + //===------------------------------------------------------------------------------------------- + template KYOSU_FORCEINLINE + auto sb_kn(int n, Z z) noexcept + { + using u_t = eve::underlying_type_t; + auto iton = [](int n){ + if (n%4 == 0) return complex(eve::one(eve::as())); + else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::one(eve::as())); + else if (n%4 == 2) return complex(eve::mone(eve::as())); + else return complex(eve::zero(eve::as()), eve::mone(eve::as())); + }; + return -eve::pio_2(eve::as())*iton(n)*sph_bessel_h1n(n, complex(-imag(z), real(z))); + } + + template KYOSU_FORCEINLINE + auto sb_kn(int n, Z z, R& sks) noexcept + { + using u_t = eve::underlying_type_t; + auto iton = [](int n){ + if (n%4 == 0) return complex(eve::one(eve::as())); + else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::one(eve::as())); + else if (n%4 == 2) return complex(eve::mone(eve::as())); + else return complex(eve::zero(eve::as()), eve::mone(eve::as())); + }; + return -eve::pio_2(eve::as())*iton(n)*sph_bessel_h1n(n, complex(-imag(z), real(z)), sks); + } + + template KYOSU_FORCEINLINE + auto sb_ikn(int n, Z z, R1& sis, R2& sks) noexcept + { + using u_t = eve::underlying_type_t; + auto iton = [](int n){ + if (n%4 == 0) return complex(eve::one(eve::as())); + else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::one(eve::as())); + else if (n%4 == 2) return complex(eve::mone(eve::as())); + else return complex(eve::zero(eve::as()), eve::mone(eve::as())); + }; + sph_bessel_h12n(n, complex(-imag(z), real(z)), sis, sks); + for(int i=0; i <= n; ++i) + { + auto fac = -eve::pio_2(eve::as())*iton(n); + sis[i]*= fac; + sks[i]*= fac; + } + return kumi::tuple{sis[n], sks[n]}; + } + + //===------------------------------------------------------------------------------------------- + // sb_k0 + //===------------------------------------------------------------------------------------------- + template KYOSU_FORCEINLINE + auto sb_k0(Z z) noexcept + { + using u_t = eve::underlying_type_t; + return eve::pio_2(eve::as())*exp(-z)/z; + } + + //===------------------------------------------------------------------------------------------- + // sb_k1 + //===------------------------------------------------------------------------------------------- + template + auto sb_k1(Z z) noexcept + { + using u_t = eve::underlying_type_t; + auto rz = rec(z); + return eve::pio_2(eve::as())*exp(-z)*fma(rz, rz, rz); + } + + //////////////////////////////////////////////////////////// + // end utilities + //////////////////////////////////////////////////////////// + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return sb_i2_0(z); + } + else + { + return cakley_extend(sph_bessel_i2_0, z); + } + } + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return sb_i2_1(z); + } + else + { + return cakley_extend(sph_bessel_i2_1, z); + } + } + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return sb_i1_0(z); + } + else + { + return cakley_extend(sph_bessel_i1_0, z); + } + } + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return sb_i1_1(z); + } + else + { + return cakley_extend(sph_bessel_i1_1, z); + } + } + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, N n, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return sb_i1n(n, z); + } + else + { + return cakley_extend_rev(sph_bessel_i1n, n, z); + } + } + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, N n, Z z, R& sis) noexcept + { + return sb_i1n(n, z, sis); + } + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, N n, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return sb_i2n(n, z); + } + else + { + return cakley_extend_rev(sph_bessel_i2n, n, z); + } + } + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, N n, Z z, R& sis) noexcept + { + return sb_i2n(n, z, sis); + } + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return sb_k0(z); + } + else + { + return cayley_extend(sph_bessel_k0, z); + } + } + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return sb_k1(z); + } + else + { + return cayley_extend(sph_bessel_k1, z); + } + } + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, N n, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return sb_kn(n, z); + } + else + { + return cayley_extend_rev(sph_bessel_kn, n, z); + } + } + + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, N n, Z z, R& sks) noexcept + { + return sb_kn(n, z, sks); + } + +// template KYOSU_FORCEINLINE +// auto dispatch(eve::tag_of, N n, Z z, R1& sis, R2& sks) noexcept +// { +// return sb_ikn(n, z, sis, sks); +// } +} diff --git a/test/unit/function/sph_bessel_h0.inactive b/test/unit/function/sph_bessel_h0.cpp similarity index 100% rename from test/unit/function/sph_bessel_h0.inactive rename to test/unit/function/sph_bessel_h0.cpp diff --git a/test/unit/function/sph_bessel_h1.inactive b/test/unit/function/sph_bessel_h1.cpp similarity index 100% rename from test/unit/function/sph_bessel_h1.inactive rename to test/unit/function/sph_bessel_h1.cpp diff --git a/test/unit/function/sph_bessel_hn.inactive b/test/unit/function/sph_bessel_hn.cpp similarity index 100% rename from test/unit/function/sph_bessel_hn.inactive rename to test/unit/function/sph_bessel_hn.cpp diff --git a/test/unit/function/sph_bessel_i0.inactive b/test/unit/function/sph_bessel_i0.cpp similarity index 100% rename from test/unit/function/sph_bessel_i0.inactive rename to test/unit/function/sph_bessel_i0.cpp diff --git a/test/unit/function/sph_bessel_i1.inactive b/test/unit/function/sph_bessel_i1.cpp similarity index 100% rename from test/unit/function/sph_bessel_i1.inactive rename to test/unit/function/sph_bessel_i1.cpp diff --git a/test/unit/function/sph_bessel_in.cpp b/test/unit/function/sph_bessel_in.cpp new file mode 100644 index 00000000..ebd25572 --- /dev/null +++ b/test/unit/function/sph_bessel_in.cpp @@ -0,0 +1,95 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_i1n over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + if constexpr(sizeof(T) == 8) + { + std::array re{-5.2999999999999998e+00, -4.5999999999999996e+00, -3.8999999999999999e+00, -3.2000000000000002e+00, -2.5000000000000000e+00, -1.7999999999999998e+00, -1.1000000000000005e+00, -4.0000000000000036e-01, 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00, 3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00}; + std::array im{8.0000000000000000e+00, -7.0000000000000000e+00, 6.0000000000000000e+00, -5.0000000000000000e+00, 4.0000000000000000e+00, -3.0000000000000000e+00, 2.0000000000000000e+00, -1.0000000000000000e+00, -0.0000000000000000e+00, 1.0000000000000000e+00, -2.0000000000000000e+00, 3.0000000000000000e+00, -4.0000000000000000e+00, 5.0000000000000000e+00, -6.0000000000000000e+00, 7.0000000000000000e+00}; + { + std::array reres0{7.7706189084700474e+00, 5.7193348954956411e+00, 9.9651419983548073e-01, -1.3562427634458576e+00, -1.2786706440683431e+00, -3.2086069648420901e-01, 4.6505781935223922e-01, 8.6074414214508577e-01, 1.0150676448238098e+00, 9.6671074810035795e-01, 4.7487788451853485e-01, -7.2053339349585488e-01, -2.1909908780305312e+00, -2.1079169900772778e+00, 2.1152290966576488e+00, 1.0154509786953152e+01}; + std::array imres0{-6.9698029235714092e+00, -1.5982850495762253e+00, 3.3035431065212375e+00, -1.5627342120843337e+00, -1.8950023844531727e-01, 7.7839373841130077e-01, -5.3369453536091727e-01, 1.2236722808607153e-01, 0.0000000000000000e+00, 3.3174683331562105e-01, -9.5413188180798625e-01, 1.2274152539491650e+00, -1.1199496047458246e-01, -2.8693859068716674e+00, 5.6153371971630630e+00, -2.2178922189196260e+00}; + std::array reres1{-6.7174528679796257e+00, -5.5036519799541068e+00, -1.3098238628283609e+00, 1.0047646742110703e+00, 1.1639277971364117e+00, 4.9424086716566795e-01, -1.6557649475684527e-02, -9.7080721799298941e-02, 1.0090289768350588e-01, 2.6208507473901876e-01, 1.3771516732295430e-02, -8.7005346095573977e-01, -1.9415108978882807e+00, -1.5377246143859755e+00, 2.5461644883125758e+00, 9.6642158761147634e+00}; + std::array imres1{7.2434985699405363e+00, 9.2405428066530826e-01, -2.9370454059224973e+00, 1.6113493268362018e+00, -4.5039559309533850e-02, -6.2196029018627774e-01, 5.8880997583441919e-01, -3.1547053555014831e-01, 0.0000000000000000e+00, 3.9506579770822658e-01, -8.3764199676086581e-01, 8.9755809983712842e-01, 2.3504424348042030e-01, -2.8588922477687770e+00, 4.9413752253773460e+00, -1.1320604331709858e+00}; + std::array reres6{-2.6240029171547561e+00, -5.1076195004248492e-01, 1.2968488456394686e-01, 1.3990664186260438e-01, 5.3442223694414565e-02, 1.0867465917187520e-02, 9.1904729456308994e-04, 7.1102071274073375e-06, 5.4108106474652957e-09, 3.9446343583993648e-06, 1.5302374589411690e-03, 1.9752158517365788e-02, 1.0456001605232421e-01, 3.0094545968908643e-01, 3.9527685917385558e-01, -5.5947069481577105e-01}; + std::array imres6{-1.3782512229893624e+00, 1.2635768993261309e+00, -5.6286102987787712e-01, 1.5843241471028197e-01, -2.7888540684483137e-02, 2.9738694621436607e-03, -2.5713535324299396e-04, 8.6945533993184342e-06, 0.0000000000000000e+00, -5.9084004894842560e-05, 1.7696008192200146e-03, -8.8585161213433159e-03, -3.0090818333100081e-03, 1.8254279393427686e-01, -8.6690629805255193e-01, 2.1924764971833897e+00}; + for(int i=0; i < 16; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res0 = kyosu::complex(reres0[i], imres0[i]); + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1n(0, c), res0, 1.0e-7) << i << " <- " << c << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1n(0, c), kyosu::sph_bessel_i1_0(c), 1.0e-7) << i << " <- " << c << '\n'; + auto res1 = kyosu::complex(reres1[i], imres1[i]); + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1n(1, c), res1, 1.0e-7) << i << " <- " << c << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1n(1, c), kyosu::sph_bessel_i1_1(c), 1.0e-7) << i << " <- " << c << '\n'; + auto res6 = kyosu::complex(reres6[i], imres6[i]); + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1n(6, c), res6, 1.0e-7) << i << " <- " << c << '\n'; + } + } + { + std::array reres0{-7.7701480965979615e+00, -5.7191728632529424e+00, -9.9865703486958179e-01, 1.3496468508128983e+00, 1.2735311974235841e+00, 3.5064348468508577e-01, -3.1961958268284379e-01, -4.9937803640628725e-01, 3.4844617137628719e+00, 9.1131386544700832e-01, 4.0790158455486170e-01, -7.3773871722556050e-01, -2.1892302215898161e+00, -2.1045860265728922e+00, 2.1164135150144610e+00, 1.0154460554701284e+01}; + std::array imres0{6.9695817953414272e+00, 1.5994741138540221e+00, -3.3053895440682828e+00, 1.5608253671722003e+00, 2.0612597652806902e-01, -8.1507228516936348e-01, 5.2296512942223727e-01, 3.8435468075214141e-01, 0.0000000000000000e+00, 7.7583840315858354e-02, -9.3521361215077858e-01, 1.2435876926532230e+00, -1.2072100015852272e-01, -2.8681235230153779e+00, 5.6162266381272072e+00, -2.2185229272143054e+00}; + std::array reres1{6.7179699859518553e+00, 5.5039432710089606e+00, 1.3077341730520082e+00, -1.0122303696188917e+00, -1.1726336064302898e+00, -4.6906809706222657e-01, 1.9682138944905253e-01, 1.0198850612740733e+00, -1.0599804734385765e+01, 4.7226189521892448e-01, 1.0276468886616588e-01, -8.5333760464885822e-01, -1.9448475868827635e+00, -1.5415365510292587e+00, 2.5449801901900613e+00, 9.6643265361121937e+00}; + std::array imres1{-7.2436915244778488e+00, -9.2280342234562251e-01, 2.9348072847855646e+00, -1.6124956501976551e+00, 6.2609412386795293e-02, 5.7258813650279483e-01, -5.4597429025558630e-01, 6.8540163214395333e-01, 0.0000000000000000e+00, 7.4861184588119567e-01, -8.4178645778090444e-01, 8.7525897928490914e-01, 2.4455154428335113e-01, -2.8598539785097512e+00, 4.9402882911780104e+00, -1.1313911260155840e+00}; + std::array reres6{2.6246088813895847e+00, 5.1558271591018912e-01, -1.1332387711202019e-01, -1.0055046114692631e-01, 2.6616265103540510e-03, -2.6736332326717582e-01, -6.2624240645191396e+00, 3.1625262777005564e+03, 4.7336904824867837e+07, 7.0532075943362565e+02, 1.2546826026901138e+01, 7.5732241106251996e-01, 1.5046516286344094e-01, 2.9063629792862511e-01, 3.8744481722959184e-01, -5.6229910703386254e-01}; + std::array imres6{1.3800676683211495e+00, -1.2660426058258256e+00, 5.5883548768306357e-01, -1.1196873227250402e-01, -2.4894865246898401e-01, 2.0323483453640900e+00, -3.5895780247972894e+01, 5.5916670355854685e+03, 0.0000000000000000e+00, 5.8750719298833349e+02, 1.1624042874403822e+00, -5.7381216048713535e-01, 1.3154820430720915e-01, 1.5222870262636137e-01, -8.6241979472821317e-01, 2.1932193362333159e+00}; + + for(int i=0; i < 16; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res0 = kyosu::complex(reres0[i], imres0[i]); + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2n(0, c), res0, 1.0e-7) << i << " <- " << c << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2n(0, c), kyosu::sph_bessel_i2_0(c), 1.0e-7) << i << " <- " << c << '\n'; + auto res1 = kyosu::complex(reres1[i], imres1[i]); + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2n(1, c), res1, 1.0e-7) << i << " <- " << c << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2n(1, c), kyosu::sph_bessel_i2_1(c), 1.0e-7) << i << " <- " << c << '\n'; + auto res6 = kyosu::complex(reres6[i], imres6[i]); + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2n(6, c), res6, 1.0e-7) << i << " <- " << c << '\n'; + } + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_jn over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + if constexpr(sizeof(eve::underlying_type_t) == 8) + { + using u_t = eve::underlying_type_t; + auto c = kyosu::complex(a0, a1); + auto cm = -c; + auto zer = kyosu::complex(T(0), T(0)); + // auto one = kyosu::complex(T(1), T(0)); + + using kyosu::sph_bessel_i1n; + using kyosu::sph_bessel_i2n; + + for(int i=1; i < 4 ; ++i) + { + auto i1nc = sph_bessel_i1n(i, c); + auto i2nc = sph_bessel_i2n(i, c); + TTS_IEEE_EQUAL(i1nc, eve::sign_alternate(u_t(i))*sph_bessel_i1n(i, cm)) << "c " << c << "\n"; + TTS_IEEE_EQUAL(i2nc, eve::sign_alternate(u_t(i+1))*sph_bessel_i2n(i, cm)) << "c " << c << "\n"; + TTS_EXPECT(eve::all(kyosu::is_eqz(sph_bessel_i1n(i, zer)))) << i << sph_bessel_i1n(i, zer) << '\n'; + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_i2n(i, zer)))) << i << sph_bessel_i2n(i, zer) << '\n'; + } + } +}; diff --git a/test/unit/function/sph_bessel_in.inactive b/test/unit/function/sph_bessel_in.inactive deleted file mode 100644 index 0b50d34d..00000000 --- a/test/unit/function/sph_bessel_in.inactive +++ /dev/null @@ -1,76 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#include -#include - -TTS_CASE_WITH ( "Check kyosu::sph_bessel_i1nn over real" - , kyosu::scalar_real_types - , tts::generate(tts::randoms(-10,10)) - ) -(T ) -{ - if constexpr(sizeof(T) == 8) - { - { - std::array re{-1.3102964777073500e+00, 2.6429797151271517e+00, -1.2991425964097203e+00, -2.0616201216776595e+00, -3.5692327918927735e-01, 4.5409169416450617e+00, 1.3915394762511490e+00, -1.5742664543118334e+00}; - std::array im{-4.9978344205200456e+00, 4.6966388867930240e+00, 4.5575687436111947e+00, -4.3856054937978231e+00, 1.7496198325295331e+00, -4.8976206253867760e+00, -2.2292456733273172e+00, -3.6547854593577123e+00}; - std::array reres{-8.9448633564799367e-04, -1.3550451109077518e-02, -1.4079674805043841e-03, -5.9344800455695375e-03, -1.6670926554323534e-08, 1.0879026036746578e-02, -2.4460612590389117e-05, -1.3445907340862448e-03}; - std::array imres{-7.8433567383746679e-03, 3.9123036962636004e-03, 4.2071250786199152e-03, -7.1820574884934665e-04, 2.7760904475739103e-06, -1.1171678180080313e-01, -5.5883425389321553e-05, -2.7135427306845860e-04}; - for(int i=0; i < 8; ++i) - { - auto c = kyosu::complex(re[i], im[i]); - auto res = kyosu::complex(reres[i], imres[i]); - TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1n(8, c), res, 1.0e-7) << i << " <- " << c << '\n'; - } - } - { - std::array re{-1.0353298912411590e+00, -4.0590492982139348e-01, -3.7533860841866939e+00, 2.5026965372926147e-01, 1.3027543743434511e-01,-3.9640689770100854e+00, -1.0995476278801330e+00, -3.8763377787749520e-01}; - std::array im{-4.2831047668722677e+00, 7.7061553129385785e-01, 1.3657001488056097e+00, 4.1210952543771051e-01, 4.4818969533279871e+00, -3.7525013370467541e+00, 4.8289663184328813e-01, -7.0094140409912353e-01}; - std::array reres{-5.7006149471406644e+00, 6.6457154244604167e+06, 4.9722413633652538e+00, -1.4179514732481084e+09, 1.2132650006681511e+00, -4.7090187937440880e-02, 3.0782381172903924e+05, 1.4851214398026228e+07}; - std::array imres{-1.3655417114292883e+00, 2.5782447371366373e+06, 1.6560730551132823e+00, -2.7650888434981310e+08, -5.5042068213386841e+00, 4.0322394841215375e-01, 2.1882447301295877e+05, -2.7656213800463108e+06}; - for(int i=0; i < 8; ++i) - { - auto c = kyosu::complex(re[i], im[i]); - auto res = kyosu::complex(reres[i], imres[i]); - TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2n(8, c), res, 1.0e-7) << i << " <- " << c << '\n'; - } - } - } -}; - - -TTS_CASE_WITH ( "Check kyosu::sph_bessel_jn over real" - , kyosu::real_types - , tts::generate(tts::randoms(-10,10), - tts::randoms(-10,10) - ) - ) -(T a0, T a1) -{ - if constexpr(sizeof(eve::underlying_type_t) == 8) - { - using u_t = eve::underlying_type_t; - auto c = kyosu::complex(a0, a1); - auto cm = -c; - auto zer = kyosu::complex(T(0), T(0)); - // auto one = kyosu::complex(T(1), T(0)); - - using kyosu::sph_bessel_i1n; - using kyosu::sph_bessel_i2n; - - for(int i=1; i < 4 ; ++i) - { - auto i1nc = sph_bessel_i1n(i, c); - auto i2nc = sph_bessel_i2n(i, c); - TTS_IEEE_EQUAL(i1nc, eve::sign_alternate(u_t(i))*sph_bessel_i1n(i, cm)) << "c " << c << "\n"; - TTS_IEEE_EQUAL(i2nc, eve::sign_alternate(u_t(i+1))*sph_bessel_i2n(i, cm)) << "c " << c << "\n"; - TTS_EXPECT(eve::all(kyosu::is_eqz(sph_bessel_i1n(i, zer)))) << i << sph_bessel_i1n(i, zer) << '\n'; - TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_i2n(i, zer)))) << i << sph_bessel_i2n(i, zer) << '\n'; - } - } -}; diff --git a/test/unit/function/sph_bessel_k0.inactive b/test/unit/function/sph_bessel_k0.cpp similarity index 100% rename from test/unit/function/sph_bessel_k0.inactive rename to test/unit/function/sph_bessel_k0.cpp diff --git a/test/unit/function/sph_bessel_k1.inactive b/test/unit/function/sph_bessel_k1.cpp similarity index 100% rename from test/unit/function/sph_bessel_k1.inactive rename to test/unit/function/sph_bessel_k1.cpp diff --git a/test/unit/function/sph_bessel_kn.inactive b/test/unit/function/sph_bessel_kn.cpp similarity index 100% rename from test/unit/function/sph_bessel_kn.inactive rename to test/unit/function/sph_bessel_kn.cpp