From c721ab9281f33184dea7ec1df7e4a92bc680e86d Mon Sep 17 00:00:00 2001 From: jtlap Date: Sat, 16 Dec 2023 13:34:01 +0100 Subject: [PATCH] airy_ai airy_bi + bug bessel_j for v <0 --- include/kyosu/functions/cyl_bessel_y.hpp | 14 ++ include/kyosu/types/impl/besselr/airy.hpp | 49 ++---- include/kyosu/types/impl/besselr/cb_jyr.hpp | 43 +++-- test/doc/cyl_bessel_j.cpp | 13 +- test/unit/function/airy_ai.cpp | 34 +--- test/unit/function/airy_bi.cpp | 4 +- test/unit/function/cyl_bessel_i.cpp | 1 - test/unit/function/cyl_bessel_j.cpp | 119 ++++++++++++++ test/unit/function/cyl_bessel_k.cpp | 165 ++++++++++---------- test/unit/function/cyl_bessel_y.cpp | 102 ++++++++++++ 10 files changed, 363 insertions(+), 181 deletions(-) create mode 100644 test/unit/function/cyl_bessel_j.cpp create mode 100644 test/unit/function/cyl_bessel_y.cpp diff --git a/include/kyosu/functions/cyl_bessel_y.hpp b/include/kyosu/functions/cyl_bessel_y.hpp index 2b4c0e0e..33df0eea 100644 --- a/include/kyosu/functions/cyl_bessel_y.hpp +++ b/include/kyosu/functions/cyl_bessel_y.hpp @@ -25,6 +25,13 @@ namespace kyosu::tags return fnu(nu, complex(v)); } + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N nu, T const& v, R& ys) noexcept + { + auto fnu = callable_cyl_bessel_y{}; + return fnu(nu, complex(v), ys); + } + template KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1) const noexcept -> decltype(eve::tag_invoke(*this, target0, target1)) @@ -32,6 +39,13 @@ namespace kyosu::tags return eve::tag_invoke(*this, target0, target1); } + template + KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1, R & out0) const noexcept + -> decltype(eve::tag_invoke(*this, target0, target1, out0)) + { + return eve::tag_invoke(*this, target0, target1, out0); + } + template eve::unsupported_call operator()(T&&... x) const requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; diff --git a/include/kyosu/types/impl/besselr/airy.hpp b/include/kyosu/types/impl/besselr/airy.hpp index 6ade4c39..60001def 100644 --- a/include/kyosu/types/impl/besselr/airy.hpp +++ b/include/kyosu/types/impl/besselr/airy.hpp @@ -46,27 +46,24 @@ namespace kyosu::_ auto third = eve::third(as()); auto invpi = eve::inv_pi(as()); auto br_im_eq_0 = [](auto z){ -// std::cout << " === 3 "<< std::endl; -// return Z(3.550280538878172e-01); return Z(eve::airy_ai(real(z))); }; + auto br_re_lt_0 = [third, invpi](auto z){ auto izgt0 = eve::is_gez(imag(z)); z = if_else(izgt0, z, conj(z)); auto zet = 2*pow(z, u_t(1.5) )/3; -// std::cout << " === 2 "<< z << "arg z " << arg(z) << " arg zet " << arg(zet) << std::endl; auto argzetgtpio2 = arg(zet) > eve::pio_2(eve::as< u_t>()); - auto sqtmz = sqrt(-z/3); //pow(z, u_t(0.5)); + auto sqtmz = sqrt(-z/3); auto r = mulmi(invpi*sqtmz*cyl_bessel_k(third, zet)); -// return if_else(imag(z) > 0, r, -r); r = if_else(izgt0, r, conj(r)); return if_else(argzetgtpio2, -r, r); -// return (cyl_bessel_i(third, zet)+cyl_bessel_i(-third, zet))*sqrt(-z)*third; }; auto br_re_gt_0 = [third, invpi](auto z){ -// std::cout << " === 1 "<< std::endl; auto zet = 2*pow(z, u_t(1.5) )/3; - return invpi*sqrt(z/3)*cyl_bessel_k(third, zet); + std::array ks; + auto r = invpi*sqrt(z/3)*cyl_bessel_k(third, zet, ks); + return r; }; auto notdone = eve::true_(as()); @@ -84,14 +81,7 @@ namespace kyosu::_ } } } - return r; -// using u_t = eve::underlying_type_t; -// auto [sqzo3, zeta] = zet(z); -// auto r = eve::inv_pi(as(z))*sqzo3*cyl_bessel_k(eve::third(as()), zeta); -// return r; -// return if_else(eve::is_ltz(arg(z)), -r, r); -// return sqzo3*(cyl_bessel_i(eve::third(as()), zeta) -// +cyl_bessel_i(-eve::third(as()), zeta)); + return r; } //===------------------------------------------------------------------------------------------- @@ -104,28 +94,15 @@ namespace kyosu::_ auto third = eve::third(as()); auto invpi = eve::inv_pi(as()); auto br_im_eq_0 = [](auto z){ -// std::cout << " === 3 "<< std::endl; -// return Z(3.550280538878172e-01); return Z(eve::airy_bi(real(z))); }; + auto br_re_lt_0 = [third, invpi](auto z){ - std::cout << " === 1 "<< std::endl; auto zet = 2*pow(-z, u_t(1.5) )/3; - std::cout << " === 2 "<< z << "arg z " << arg(z) << " arg zet " << arg(zet) << std::endl; - return (cyl_bessel_j(third, zet)+cyl_bessel_j(-third, zet))*sqrt(-z/3) ; -// auto izgt0 = eve::is_gez(imag(z)); -// z = if_else(izgt0, z, conj(z)); -// auto zet = 2*pow(z, u_t(1.5) )/3; -// auto argzetgtpio2 = arg(zet) > eve::pio_2(eve::as< u_t>()); -// auto sqtmz = sqrt(-z/3); //pow(z, u_t(0.5)); -// auto r = mulmi(invpi*sqtmz*cyl_bessel_k(third, zet)); -// r = if_else(izgt0, r, conj(r)); -// return if_else(argzetgtpio2, -r, r); - + return (kyosu::cyl_bessel_j(-third, zet)-kyosu::cyl_bessel_j(third, zet))*sqrt(-z/3) ; }; auto br_re_gt_0 = [third, invpi](auto z){ - std::cout << " === 1 "<< std::endl; auto zet = 2*pow(z, u_t(1.5) )/3; return sqrt(z/3)*(cyl_bessel_i(third, zet)+cyl_bessel_i(-third, zet)); }; @@ -145,15 +122,7 @@ namespace kyosu::_ } } } - return r; - // using u_t = eve::underlying_type_t; -// auto [sqzo3, zeta] = zet(z); -// // std::cout << "z" << z << " sqzo3 "<< sqzo3 << " zeta "<< zeta << std::endl; -// // std::cout << "third " << eve::third(as()) << std::endl; -// // std::cout << "cyl_bessel_i(eve::third(as()), zeta) " << cyl_bessel_i(eve::third(as()), zeta) << std::endl; -// // std::cout << "cyl_bessel_i(-eve::third(as()), zeta) " << cyl_bessel_i(-eve::third(as()), zeta) << std::endl; -// return sqzo3*(cyl_bessel_i(-eve::third(as()), zeta) -// -cyl_bessel_i(eve::third(as()), zeta)); + return r; } //===------------------------------------------------------------------------------------------- diff --git a/include/kyosu/types/impl/besselr/cb_jyr.hpp b/include/kyosu/types/impl/besselr/cb_jyr.hpp index 0494dec2..e38f3cac 100644 --- a/include/kyosu/types/impl/besselr/cb_jyr.hpp +++ b/include/kyosu/types/impl/besselr/cb_jyr.hpp @@ -40,12 +40,13 @@ namespace kyosu::_ cjv[1] = cjv1; cyv[1] = cyv1; } - if(n == 0) { - return kumi::tuple{cjv0, cyv0}; + if(n == 0) + { + return kumi::tuple{cjv[0], cyv[0]}; } else if (n == 1) { - return kumi::tuple{cjv1, cyv1}; + return kumi::tuple{cjv[1], cyv[1]}; } else if (n >= 2) { @@ -126,16 +127,24 @@ namespace kyosu::_ else { // v is not flint and is less than 0; + // we use: + // J_{-v] = J_v(z) cos(v pi) - Y_v *sin(v pi) + // y_{-v] = J_v(z) sin(v pi) + Y_v *cos(v pi) cb_jy(-v, z, cjv, cyv); - for(int i=0; i < an ; ++i) + auto v0 = frac(-v); + auto [s, c] = sinpicospi(v0); + for(int i=0; i <= an ; ++i) { - auto [s, c] = sinpicospi(v); - auto a = cjv[i]*s; - auto b = cyv[i]*c; - cjv[i] = a-b; - cyv[i] = a+b; + auto aa = cjv[i]*c; + auto bb = cyv[i]*s; + auto cc = cjv[i]*s; + auto dd = cyv[i]*c; + cjv[i] = aa-bb; //mulmi(a-b); + cyv[i] = cc+dd; + s = -s; + c = -c; } - return kumi::tuple{cjv[n], cyv[n]}; + return kumi::tuple{cjv[an], cyv[an]}; } } @@ -154,10 +163,10 @@ namespace kyosu::_ template auto dispatch(eve::tag_of, N v, Z z, R & js) noexcept { - auto n = int(abs(v))+1; - auto doit = [v, z, &js](auto ys){ + auto n = int(abs(v)); + auto doit = [n, v, z, &js](auto ys){ auto [jvn, _] = cb_jy(v, z, js, ys); - return jvn; + return js[n]; }; return with_alloca(n+1, doit); } @@ -168,12 +177,12 @@ namespace kyosu::_ template auto dispatch(eve::tag_of, N v, Z z) noexcept { - auto n = int(abs(v))+1; + auto n = int(abs(v)); if constexpr(concepts::complex ) { - auto doit = [v, z](auto js, auto ys){ + auto doit = [n, v, z](auto js, auto ys){ auto [jv, _] = cb_jy(v, z, js, ys); - return jv; + return js[n]; }; return with_alloca(n+1, doit); } @@ -203,7 +212,7 @@ namespace kyosu::_ template auto dispatch(eve::tag_of, N v, Z z) noexcept { - auto n = int(abs(v))+1; + auto n = int(abs(v)); if constexpr(concepts::complex ) { auto doit = [v, z](auto js, auto ys){ diff --git a/test/doc/cyl_bessel_j.cpp b/test/doc/cyl_bessel_j.cpp index 4200b279..0d86ba39 100644 --- a/test/doc/cyl_bessel_j.cpp +++ b/test/doc/cyl_bessel_j.cpp @@ -7,13 +7,18 @@ int main() std::cout.precision(16); using w_t = eve::wide>; auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); - auto v = 3.3; + auto v = -(2+1/3.0); + int nb = int(eve::abs(v)+1); std::cout << "z " << z << std::endl; - std::vector js(4); + std::vector js(nb); kyosu::cyl_bessel_j(v, z, js); - for(int n=0; n <= 3; ++n) + auto inc = eve::sign(v); + auto v0 = eve::frac(v); + for(int n=0; n < nb; ++n) { - std::cout << "js[" << n << "] = " << js[n] << std::endl; + std::cout << "js[" << n << "] = " << js[n] << std::endl; + std::cout << "j(" << v0 << ", z) " << kyosu::cyl_bessel_j(v0, z) << std::endl; + v0 += inc; } return 0; } diff --git a/test/unit/function/airy_ai.cpp b/test/unit/function/airy_ai.cpp index 520a1a86..98d10292 100644 --- a/test/unit/function/airy_ai.cpp +++ b/test/unit/function/airy_ai.cpp @@ -41,45 +41,15 @@ TTS_CASE_WITH ( "Check kyosu::airy_ai over real" , 0.0000000000000000e+00, -1.5188956587718155e-01, 2.1867029634299019e-02, 6.0030132776007987e-02 , 3.9632838829414818e-02, 2.2033585182583651e-02, 1.3313384345125566e-02, 8.4631437627202440e-03}; - auto z = kyosu::complex(0.0); - std::cout << kyosu::airy_ai(z) << std::endl; for(int i=0; i < 16; ++i) { auto c = kyosu::complex(re[i], im[i]); auto res = kyosu::complex(reres[i], imres[i]); -// std::cout << "i " << i << " c " << c << " arg(c) " << kyosu::arg(c) << " sqrt(c) "<< kyosu::sqrt(c) << std::endl; TTS_RELATIVE_EQUAL(kyosu::airy_ai(c), res, 1.0e-3) << i << " <- " << c << "arg " << kyosu::arg(c) << '\n'; -// auto z = 2*kyosu::pow(c, T(1.5) )/3; -// auto kz = -kyosu::cyl_bessel_k(1/T(3), z)*kyosu::sqrt(c/3)/eve::pi(eve::as(c)); -// std::cout << "kz " << kz << std::endl; -// TTS_RELATIVE_EQUAL(kz, res, 1.0e-4) << i << " <- " << c << '\n'; TTS_RELATIVE_EQUAL(kyosu::airy_ai(re[i]), kyosu::real(kyosu::airy_ai(kyosu::complex(re[i], e_t(0.0)))), 1.0e-4)<< i << " -> " << re[i] << '\n'; TTS_RELATIVE_EQUAL(kyosu::airy_ai(im[i]), kyosu::real(kyosu::airy_ai(kyosu::complex(im[i], e_t(0.0)))), 1.0e-4)<< i << " -> " << im[i] << '\n'; + TTS_RELATIVE_EQUAL(kyosu::airy_ai(kyosu::conj(c)), kyosu::conj(res), 1.0e-3) << i << " <- " << c << "arg " << kyosu::arg(c) << '\n'; + } } }; - -// TTS_CASE_WITH ( "Check kyosu::airy_ai over real" -// , kyosu::real_types -// , tts::generate(tts::randoms(-10,10), -// tts::randoms(-10,10) -// ) -// ) -// (T a0, T a1) -// { -// using u_t = eve::underlying_type_t; -// if constexpr(sizeof(u_t) == 8) -// { -// auto c = kyosu::complex(a0, a1); -// auto cb= kyosu::conj(c); -// // auto cr= kyosu::complex(a0); -// // auto ci= kyosu::complex(T(0), a1); -// using kyosu::airy_ai; -// auto j0c = airy_ai(c); -// TTS_EQUAL(j0c, kyosu::conj(airy_ai(cb))) << "c = " << c << '\n'; -// // TTS_EXPECT(eve::all(kyosu::is_real(cr))); -// // TTS_EXPECT(eve::all(kyosu::is_pure(ci))); -// // auto z = kyosu::complex(T(0), T(0)); -// // TTS_IEEE_EQUAL(airy_ai(z), kyosu::complex(airy_ai(T(0)))); -// } -// }; diff --git a/test/unit/function/airy_bi.cpp b/test/unit/function/airy_bi.cpp index 92739d80..f346e4a5 100644 --- a/test/unit/function/airy_bi.cpp +++ b/test/unit/function/airy_bi.cpp @@ -30,8 +30,8 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j0 over real" auto c = kyosu::complex(re[i], im[i]); auto res = kyosu::complex(reres[i], imres[i]); TTS_RELATIVE_EQUAL(kyosu::airy_bi(c), res, 1.0e-3) << i << " <- " << c << '\n'; -// TTS_RELATIVE_EQUAL(kyosu::airy_bi(re[i]), kyosu::real(kyosu::airy_bi(kyosu::complex(re[i], e_t(0.0)))), 1.0e-4)<< re[i] << '\n'; -// TTS_RELATIVE_EQUAL(kyosu::airy_bi(im[i]), kyosu::real(kyosu::airy_bi(kyosu::complex(im[i], e_t(0.0)))), 1.0e-4)<< im[i] << '\n'; + TTS_RELATIVE_EQUAL(kyosu::airy_bi(re[i]), kyosu::real(kyosu::airy_bi(kyosu::complex(re[i], e_t(0.0)))), 1.0e-4)<< re[i] << '\n'; + TTS_RELATIVE_EQUAL(kyosu::airy_bi(im[i]), kyosu::real(kyosu::airy_bi(kyosu::complex(im[i], e_t(0.0)))), 1.0e-4)<< im[i] << '\n'; } } }; diff --git a/test/unit/function/cyl_bessel_i.cpp b/test/unit/function/cyl_bessel_i.cpp index 33b35b11..f552918b 100644 --- a/test/unit/function/cyl_bessel_i.cpp +++ b/test/unit/function/cyl_bessel_i.cpp @@ -65,7 +65,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i over real" auto vv = T(0.3); for (n = 0; n <= N; ++n) { - std::cout << "n = "<< n << " k = "<< k << std::endl; auto c = kyosu::complex(re[k], im[k]); kyosu::cyl_bessel_i(v1, c, is_03); auto refi_03= kyosu::complex_t(reresi_03[n][k], imresi_03[n][k]); diff --git a/test/unit/function/cyl_bessel_j.cpp b/test/unit/function/cyl_bessel_j.cpp new file mode 100644 index 00000000..01ba2710 --- /dev/null +++ b/test/unit/function/cyl_bessel_j.cpp @@ -0,0 +1,119 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k 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}; + + constexpr int N = 0; + T v1 = 1.0/3; + std::cout << "v1 " << v1 << std::endl; + auto n = int(v1); + using vc_t = std::array; + using a_t = std::array; + a_t reresi_03{ + vc_t{2.8090373106452773e+02, 1.7536493250871054e+01, -3.3400462158788272e+01, -2.3524142528283910e+01 + , -9.3274540332656208e+00, -1.9656734427524503e+00, 4.3413771695853465e-01, 7.4917804310983560e-01 + , 5.8501480583737409e-01, 1.0213243250718731e+00, 1.3227852389571668e+00, 1.9818890844578396e-01 + , -6.0472153767481238e+00, -2.3374276615174153e+01, -5.1769656721627001e+01, -5.2626473557055462e+01 + } + + }; + a_t imresi_03{ + vc_t{-2.6513665958745401e+02, 1.5122859663021515e+02, -5.0648785179729913e+01, 7.0192828558358924e+00 + , 4.0222943595848140e+00, -3.8995769510156886e+00, 1.9902620949926999e+00, -7.5756275269309670e-01 + , 0.0000000000000000e+00, -1.1158707376004939e-01, 1.3330326497253218e+00, -4.1355677304829221e+00 + , 7.6943469735731851e+00, -4.6838788665148172e+00, -2.8742293351241671e+01, 1.3952363114976131e+02 + } + }; + + std::vector> is_03(N+2); + for (int k = 0; k < 16; ++k) + { + auto c = kyosu::complex(re[k], im[k]); + kyosu::cyl_bessel_j(v1, c, is_03); + auto vv = T(N+1.0/3); + for (n = 0; n <= N; ++n) + { + auto refi_03= kyosu::complex_t(reresi_03[n][k], imresi_03[n][k]); + TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-2) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; + TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_j(vv, c), 1.e-2) << "n " << n << " k " << k << " c "<< c <<'\n'; + } + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k 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}; + + constexpr int N = 0; + T v1 = -1.0/3; + std::cout << "v1 " << v1 << std::endl; + auto n = int(v1); + using vc_t = std::array; + using a_t = std::array; + a_t reresi_03{ + vc_t{-8.9163143837610050e+01, -1.2219942089268812e+02, -6.0563369116851803e+01, -1.7842181975937766e+01 + , -1.1860080210885884e+00, 2.3789187887624883e+00, 1.9232059348460351e+00, 1.1368224659946016e+00 + , 1.3432948699326057e+00, 5.5139655214064687e-01, -4.8990051277299940e-01, -3.4926635951494012e+00 + , -9.6925116789155155e+00, -1.5745152038061331e+01, -9.9352493046059109e-01, 9.4517830645215696e+01 + } + }; + a_t imresi_03{ + vc_t{-3.7583808427575212e+02, 9.0801181249849066e+01, 3.6018913915299962e+00, -1.6864299541716473e+01 + , 1.0089846321224439e+01, -3.6425542565379301e+00, 5.6128150602785309e-01, 4.7851915025080155e-01 + , 0.0000000000000000e+00, -7.8534132361994757e-01, 1.7555311861860581e+00, -2.2254640329062294e+00 + , -1.3912074018251386e+00, 1.7899970399493917e+01, -5.9204417397002132e+01, 1.1533747474469016e+02 + } + }; + + std::vector> is_03(N+2); + for (int k = 0; k < 16; ++k) + { + auto c = kyosu::complex(re[k], im[k]); + kyosu::cyl_bessel_j(v1, c, is_03); + auto vv = -T(N+1.0/3); + for (n = 0; n <= N; ++n) + { + auto refi_03= kyosu::complex_t(reresi_03[n][k], imresi_03[n][k]); + TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-2) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; + TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_j(vv, c), 1.e-2); + } + } + } +}; diff --git a/test/unit/function/cyl_bessel_k.cpp b/test/unit/function/cyl_bessel_k.cpp index fe438422..99cfc8de 100644 --- a/test/unit/function/cyl_bessel_k.cpp +++ b/test/unit/function/cyl_bessel_k.cpp @@ -9,73 +9,74 @@ #include -// TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k 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}; - -// constexpr int N = 2; -// T v1 = N+0.3; -// auto n = int(v1); -// using vc_t = std::array; -// using a_t = std::array; -// a_t reresi_03{ -// vc_t{-7.6356085800614579e+01, -9.1777541957610591e+00, 1.6542613518775354e+01, 1.2381215528291099e+01 -// , 2.3387662997005862e+00, -2.5690360792159450e+00, -2.4930621764042549e+00, -5.9179800603848065e-01 -// , 1.4823411623387830e+00, 7.6222886510162097e-02, -1.0353447894367163e-01, -5.4817788296649773e-02 -// , -6.5755373416543517e-03, 7.4754253132452942e-03, 4.9759020245597188e-03, 8.9580881413898973e-04}, -// vc_t{-7.4858670276675767e+01, -1.2123306468557997e+01, 1.3991108787767311e+01, 1.1826801567594778e+01 -// , 3.1163376758999819e+00, -1.6392538522424551e+00, -2.2158615052476724e+00, -1.6399139926251056e+00 -// , 5.0252049759614961e+00, -3.0710975110808281e-02, -1.4448722663392172e-01, -5.9358498236853956e-02 -// , -4.2856644566998800e-03, 8.8712819627955915e-03, 5.2227131524972374e-03, 7.8878939645721496e-04}, -// vc_t{-7.0028359986020078e+01, -1.7325314484522657e+01, 8.6039335006386377e+00, 1.0034374609472982e+01 -// , 4.1203503622633431e+00, 1.8737546605673755e-01, -7.8661759114041829e-01, -1.1787359554492010e+00 -// , 4.5034117620671751e+01, -6.6359059112780450e-01, -2.6202186702526037e-01, -6.5475381627126655e-02 -// , 3.0805234748470314e-03, 1.2366476238369981e-02, 5.6902500521605311e-03, 4.7783649284957982e-04} -// }; -// a_t imresi_03{ -// vc_t{-2.8233953538916463e+01, 4.2311180773416339e+01, -1.6375448345023578e+01, -2.8400823313968524e+00 -// , 6.6879086265570535e+00, -3.1842556027568345e+00, -3.5061121052017286e-01, 1.6999326061984978e+00 -// , 0.0000000000000000e+00, -3.6592454260706442e-01, 9.2115150427228051e-02, 1.6505655277180269e-02 -// , -2.3966513298980401e-02, 8.2061105972084727e-03, 8.6732273722515746e-04, -2.1491225231985342e-03}, -// vc_t{-2.1578471449355344e+01, 3.9375141495795027e+01, -1.6966023116719221e+01, -1.2074088510183536e+00 -// , 5.7592731043626317e+00, -3.3418934981593567e+00, 4.9100238158007220e-01, 9.1783022063328734e-01 -// , 0.0000000000000000e+00, -5.3837631537993569e-01, 8.7181646569271742e-02, 2.7319352748734618e-02 -// , -2.7099439714577039e-02, 8.0965217457305613e-03, 1.3412607263892384e-03, -2.3322553819202355e-03}, -// vc_t{-8.0970040511656194e+00, 3.2454127323908530e+01, -1.7278120141267198e+01, 1.8078762957088461e+00 -// , 3.5487990981203326e+00, -2.9510953475366790e+00, 1.5914642244066768e+00, -2.7986189543910651e+00 -// , 0.0000000000000000e+00, -1.0258894849569298e+00, 3.8996035672480382e-02, 5.9423612418072189e-02 -// , -3.4235642328623889e-02, 7.3102426796495309e-03, 2.5947374159400196e-03, -2.7525951662690611e-03} -// }; - -// std::vector> is_03(N+2); -// for (int k = 0; k < 16; ++k) -// { -// auto c = kyosu::complex(re[k], im[k]); -// kyosu::cyl_bessel_k(v1, c, is_03); -// auto vv = T(-0.3); -// for (n = 0; n <= N; ++n) -// { -// auto refi_03= kyosu::complex_t(reresi_03[n][k], imresi_03[n][k]); -// TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-6) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; -// TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_k(vv, c), 1.e-6); -// vv = eve::dec(vv); -// } -// } -// } -// }; +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real v = 2.3" + , 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}; + + constexpr int N = 2; + T v1 = N+0.3; + auto n = int(v1); + using vc_t = std::array; + using a_t = std::array; + a_t reresi_03{ + vc_t{-7.6356085800614579e+01, -9.1777541957610591e+00, 1.6542613518775354e+01, 1.2381215528291099e+01 + , 2.3387662997005862e+00, -2.5690360792159450e+00, -2.4930621764042549e+00, -5.9179800603848065e-01 + , 1.4823411623387830e+00, 7.6222886510162097e-02, -1.0353447894367163e-01, -5.4817788296649773e-02 + , -6.5755373416543517e-03, 7.4754253132452942e-03, 4.9759020245597188e-03, 8.9580881413898973e-04}, + vc_t{-7.4858670276675767e+01, -1.2123306468557997e+01, 1.3991108787767311e+01, 1.1826801567594778e+01 + , 3.1163376758999819e+00, -1.6392538522424551e+00, -2.2158615052476724e+00, -1.6399139926251056e+00 + , 5.0252049759614961e+00, -3.0710975110808281e-02, -1.4448722663392172e-01, -5.9358498236853956e-02 + , -4.2856644566998800e-03, 8.8712819627955915e-03, 5.2227131524972374e-03, 7.8878939645721496e-04}, + vc_t{-7.0028359986020078e+01, -1.7325314484522657e+01, 8.6039335006386377e+00, 1.0034374609472982e+01 + , 4.1203503622633431e+00, 1.8737546605673755e-01, -7.8661759114041829e-01, -1.1787359554492010e+00 + , 4.5034117620671751e+01, -6.6359059112780450e-01, -2.6202186702526037e-01, -6.5475381627126655e-02 + , 3.0805234748470314e-03, 1.2366476238369981e-02, 5.6902500521605311e-03, 4.7783649284957982e-04} + }; + a_t imresi_03{ + vc_t{-2.8233953538916463e+01, 4.2311180773416339e+01, -1.6375448345023578e+01, -2.8400823313968524e+00 + , 6.6879086265570535e+00, -3.1842556027568345e+00, -3.5061121052017286e-01, 1.6999326061984978e+00 + , 0.0000000000000000e+00, -3.6592454260706442e-01, 9.2115150427228051e-02, 1.6505655277180269e-02 + , -2.3966513298980401e-02, 8.2061105972084727e-03, 8.6732273722515746e-04, -2.1491225231985342e-03}, + vc_t{-2.1578471449355344e+01, 3.9375141495795027e+01, -1.6966023116719221e+01, -1.2074088510183536e+00 + , 5.7592731043626317e+00, -3.3418934981593567e+00, 4.9100238158007220e-01, 9.1783022063328734e-01 + , 0.0000000000000000e+00, -5.3837631537993569e-01, 8.7181646569271742e-02, 2.7319352748734618e-02 + , -2.7099439714577039e-02, 8.0965217457305613e-03, 1.3412607263892384e-03, -2.3322553819202355e-03}, + vc_t{-8.0970040511656194e+00, 3.2454127323908530e+01, -1.7278120141267198e+01, 1.8078762957088461e+00 + , 3.5487990981203326e+00, -2.9510953475366790e+00, 1.5914642244066768e+00, -2.7986189543910651e+00 + , 0.0000000000000000e+00, -1.0258894849569298e+00, 3.8996035672480382e-02, 5.9423612418072189e-02 + , -3.4235642328623889e-02, 7.3102426796495309e-03, 2.5947374159400196e-03, -2.7525951662690611e-03} + }; + + std::vector> ks_03(N+2); + for (int k = 0; k < 16; ++k) + { + auto c = kyosu::complex(re[k], im[k]); + kyosu::cyl_bessel_k(v1, c, ks_03); + auto vv = T(-0.3); + for (n = 0; n <= N; ++n) + { + auto refi_03= kyosu::complex_t(reresi_03[n][k], imresi_03[n][k]); + TTS_RELATIVE_EQUAL(refi_03, ks_03[n], 1.e-6) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; + TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_k(vv, c), 1.e-6); + TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_k(-vv, c), 1.e-6); + vv = eve::dec(vv); + } + } + } +}; @@ -115,16 +116,16 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real" , -2.4007330768135104e-02, 8.2054910688589972e-03, 8.7325515356663085e-04, -2.1515168333318974e-03}, }; - std::vector> is_03(N+2); + std::vector> ks_03(N+2); for (int k = 0; k < 16; ++k) { auto c = kyosu::complex(re[k], im[k]); - kyosu::cyl_bessel_k(v1, c, is_03); + kyosu::cyl_bessel_k(v1, c, ks_03); auto vv = T(N+1.0/3); for (n = 0; n <= N; ++n) { auto refi_03= kyosu::complex_t(reresi_03[n][k], imresi_03[n][k]); - TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-6) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; + TTS_RELATIVE_EQUAL(refi_03, ks_03[n], 1.e-6) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_k(vv, c), 1.e-6); // auto z = 2*kyosu::pow(c, T(1.5) )/3; @@ -166,21 +167,18 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real" vc_t{1.6213718228760029e-03, -4.0459476346553763e-03, 1.9336808997634971e-03, 1.4818266297842454e-02, -4.5874488719344970e-02, 3.8478715861459151e-02, 1.5990522343845420e-01, -7.8491087169627471e-01, -3.2078969068337546e+00, 2.4846307644453542e+00, -1.2426224367177374e-01, -5.8516695697616443e+00, 1.1525345800362158e+01, -4.0835724859787863e+00, -3.0165300008239711e+01, 7.4997855427520179e+01} }; - std::vector> is_03(N+2); + std::vector> ks_03(N+2); for (int k = 0; k < 16; ++k) { auto c = -kyosu::complex(re[k], im[k]); - kyosu::cyl_bessel_k(v1, c, is_03); + kyosu::cyl_bessel_k(v1, c, ks_03); auto vv = T(N+1.0/3); for (n = 0; n <= N; ++n) { auto refi_03= kyosu::complex_t(reresi_03[n][k], imresi_03[n][k]); - TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-2) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; + TTS_RELATIVE_EQUAL(refi_03, ks_03[n], 1.e-2) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_k(vv, c), 1.e-2); - -// auto z = 2*kyosu::pow(c, T(1.5) )/3; -// std::cout << "c " << c << " z " << z << " -> " << kyosu::cyl_bessel_k(vv, z) << std::endl; -// vv = eve::dec(vv); + TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_k(-vv,c), 1.e-2); } } } @@ -204,7 +202,7 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real" , -4.0000000000000000e+00, 5.0000000000000000e+00, -6.0000000000000000e+00, 7.0000000000000000e+00}; constexpr int N = 0; - T v1 = N+1.0/3; + T v1 = N-1.0/3; std::cout << "v1 " << v1 << std::endl; auto n = int(v1); using vc_t = std::array; @@ -216,21 +214,18 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real" vc_t{2.8143596488672866e+01, -4.2273285530950140e+01, 1.6385214857913763e+01, 2.8174548075889092e+00, -6.6761389661094128e+00, 3.1878291715173703e+00, 3.3868350891838195e-01, -1.6945714309783351e+00, 0.0000000000000000e+00, 3.6796936276059572e-01, -9.2097126162947351e-02, -1.6635565810036414e-02, 2.4007330768135104e-02, -8.2054910688589972e-03, -8.7325515356663085e-04, 2.1515168333318974e-03} }; - std::vector> is_03(N+2); + std::vector> ks_03(N+2); for (int k = 0; k < 16; ++k) { auto c = kyosu::complex(re[k], -im[k]); - kyosu::cyl_bessel_k(v1, c, is_03); + kyosu::cyl_bessel_k(v1, c, ks_03); auto vv = T(N+1.0/3); for (n = 0; n <= N; ++n) { auto refi_03= kyosu::complex_t(reresi_03[n][k], imresi_03[n][k]); - TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-2) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; + TTS_RELATIVE_EQUAL(refi_03, ks_03[n], 1.e-2) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_k(vv, c), 1.e-2); - -// auto z = 2*kyosu::pow(c, T(1.5) )/3; -// std::cout << "c " << c << " z " << z << " -> " << kyosu::cyl_bessel_k(vv, z) << std::endl; -// vv = eve::dec(vv); + TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_k(-vv, c), 1.e-2); } } } diff --git a/test/unit/function/cyl_bessel_y.cpp b/test/unit/function/cyl_bessel_y.cpp new file mode 100644 index 00000000..8ccbbc20 --- /dev/null +++ b/test/unit/function/cyl_bessel_y.cpp @@ -0,0 +1,102 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k 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}; + + constexpr int N = 0; + T v1 = 1.0/3; + std::cout << "v1 " << v1 << std::endl; + auto n = int(v1); + using vc_t = std::array; + using a_t = std::array; + a_t reresi_03{ + vc_t{2.6513657493935034e+02, 1.5122843619345218e+02, 5.0648789106855794e+01, 7.0207071123160771e+00, -4.0157239964866642e+00, -3.8818208974565880e+00, -1.9700774006295092e+00, -8.8015136866517407e-01, -1.2133448538831417e+00, -4.7036021607109069e-02, 1.3293987996432377e+00, 4.1474049533381976e+00, 7.7005870282777540e+00, 4.6857906393290154e+00, -2.8742001471989823e+01, -1.3952369860713620e+02} + }; + a_t imresi_03{ + vc_t{2.8090371647178273e+02, -1.7536301901048642e+01, -3.3401195686627879e+01, 2.3525800606543957e+01, -9.3284782480155641e+00, 1.9546375587054139e+00, 5.0096629910926505e-01, -9.8992537961477600e-01, 0.0000000000000000e+00, 8.4240922212198033e-01, -1.2574860466731328e+00, 1.8207337449423835e-01, 6.0487612323155489e+00, -2.3373344181702222e+01, 5.1769001839281714e+01, -5.2626238180368290e+01} + }; + + std::vector> is_03(N+2); + for (int k = 0; k < 16; ++k) + { + auto c = kyosu::complex(re[k], im[k]); + kyosu::cyl_bessel_y(v1, c, is_03); + auto vv = T(N+1.0/3); + for (n = 0; n <= N; ++n) + { + auto refi_03= kyosu::complex_t(reresi_03[n][k], imresi_03[n][k]); + TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-2) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; + TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_y(vv, c), 1.e-2) << "n " << n << " k " << k << " c "<< c <<'\n'; + } + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k 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}; + + constexpr int N = 0; + T v1 = -1.0/3; + std::cout << "v1 " << v1 << std::endl; + auto n = int(v1); + using vc_t = std::array; + using a_t = std::array; + a_t reresi_03{ + vc_t{3.7583805458938821e+02, 9.0801266745274788e+01, -3.6012541742235697e+00, -1.6862151475581715e+01, -1.0085674143682983e+01, -3.6432335856963327e+00, -6.0906440868768552e-01, 2.0873153295804381e-01, -1.0003474349638408e-01, 8.6097480021168360e-01, 1.8102650205095943e+00, 2.2453391061314565e+00, -1.3867486242808820e+00, -1.7899822024160851e+01, -5.9204838602123715e+01, -1.1533771231556815e+02} + }; + a_t imresi_03{ + vc_t{-8.9163224441390696e+01, 1.2219965550991174e+02, -6.0563732479780832e+01, 1.7841777572774447e+01, -1.1808300271084744e+00, -2.3998139242391443e+00, 1.9741006840075483e+00, -1.1510312786004779e+00, 0.0000000000000000e+00, 3.2456737045081957e-01, 5.2569711539964548e-01, -3.4904700264222477e+00, 9.6878805608040661e+00, -1.5743030177502007e+01, 9.9294471444101262e-01, 9.4517889913758978e+01} + }; + + std::vector> is_03(N+2); + for (int k = 0; k < 16; ++k) + { + auto c = kyosu::complex(re[k], im[k]); + kyosu::cyl_bessel_y(v1, c, is_03); + auto vv = -T(N+1.0/3); + for (n = 0; n <= N; ++n) + { + auto refi_03= kyosu::complex_t(reresi_03[n][k], imresi_03[n][k]); + TTS_RELATIVE_EQUAL(refi_03, is_03[n], 1.e-2) << "n " << n << " k " << k << " c "<< c << " arg(c) "<< kyosu::arg(kyosu::conj(c)) << " < 1.57" << '\n'; + TTS_RELATIVE_EQUAL(refi_03, kyosu::cyl_bessel_y(vv, c), 1.e-2); + } + } + } +};