diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index 8d6c61b5..355dfcce 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -17,6 +17,11 @@ #include #include #include + +#include +#include +#include + #include #include #include diff --git a/include/kyosu/functions/airy.hpp b/include/kyosu/functions/airy.hpp new file mode 100644 index 00000000..fc777f30 --- /dev/null +++ b/include/kyosu/functions/airy.hpp @@ -0,0 +1,77 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_airy: eve::elementwise + { + using callable_tag_type = callable_airy; + + KYOSU_DEFERS_CALLABLE(airy_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept + { + return eve::airy(v); + } + + template + KYOSU_FORCEINLINE auto operator()(T const & target0) const noexcept + -> decltype(eve::tag_invoke(*this, target0)) + { + return eve::tag_invoke(*this, target0); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var airy +//! @brief Computes simultaneously the airy functions \f$Ai\f$ and \f$Bi\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto airy(T z) noexcept; +//! +//! template constexpr auto airy(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns a kumi pair containing \f$\{Ai(z)\}f$ and \f$\{Bi(z)\}f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/airy.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_airy airy = {}; +} diff --git a/include/kyosu/functions/airy_ai.hpp b/include/kyosu/functions/airy_ai.hpp new file mode 100644 index 00000000..d4d4a630 --- /dev/null +++ b/include/kyosu/functions/airy_ai.hpp @@ -0,0 +1,77 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_airy_ai: eve::elementwise + { + using callable_tag_type = callable_airy_ai; + + KYOSU_DEFERS_CALLABLE(airy_ai_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept + { + return eve::airy_ai(v); + } + + template + KYOSU_FORCEINLINE auto operator()(T const & target0) const noexcept + -> decltype(eve::tag_invoke(*this, target0)) + { + return eve::tag_invoke(*this, target0); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var airy_ai +//! @brief Computes the airy function \f$Ai\f$, +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto airy_ai(T z) noexcept; +//! +//! template constexpr auto airy_ai(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$\{Ai(z)\}f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/airy_ai.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_airy_ai airy_ai = {}; +} diff --git a/include/kyosu/functions/airy_bi.hpp b/include/kyosu/functions/airy_bi.hpp new file mode 100644 index 00000000..d3575bd1 --- /dev/null +++ b/include/kyosu/functions/airy_bi.hpp @@ -0,0 +1,77 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_airy_bi: eve::elementwise + { + using callable_tag_type = callable_airy_bi; + + KYOSU_DEFERS_CALLABLE(airy_bi_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept + { + return eve::airy_bi(v); + } + + template + KYOSU_FORCEINLINE auto operator()(T const & target0) const noexcept + -> decltype(eve::tag_invoke(*this, target0)) + { + return eve::tag_invoke(*this, target0); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var airy_bi +//! @brief Computes the airy function \f$Bi\f$, +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto airy_bi(T z) noexcept; +//! +//! template constexpr auto airy_bi(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$\{Bi(z)\}f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/airy_bi.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_airy_bi airy_bi = {}; +} diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp index 2d29e9b2..71866a6f 100644 --- a/include/kyosu/types/impl/bessel.hpp +++ b/include/kyosu/types/impl/bessel.hpp @@ -111,3 +111,5 @@ #include #include #include +// airy functions +#include diff --git a/include/kyosu/types/impl/besselr/airy.hpp b/include/kyosu/types/impl/besselr/airy.hpp new file mode 100644 index 00000000..6ade4c39 --- /dev/null +++ b/include/kyosu/types/impl/besselr/airy.hpp @@ -0,0 +1,220 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +namespace kyosu::_ +{ + ///////////////////////////////// + // contains implementations of + // airy_ai + // airy_bi + // airy + ///////////////////////////////// + // utilities + // zet + // ai + // bi + // aibi + ///////////////////////////////// + // implementation is done using j and y functions + + //===------------------------------------------------------------------------------------------- + // zet + //===------------------------------------------------------------------------------------------- + template KYOSU_FORCEINLINE + auto zet(Z z) noexcept + { + using u_t = eve::underlying_type_t; + auto sqz = sqrt(z); +// auto zeta = 2*sqz*z/3; + auto zeta = (pow(z, u_t(1.5))*2)/3; + return kumi::tuple{sqz/eve::sqrt_3(as(z)), zeta}; + } + + //===------------------------------------------------------------------------------------------- + // ai + //===------------------------------------------------------------------------------------------- + template KYOSU_FORCEINLINE + auto ai(Z z) noexcept + { + using u_t = eve::underlying_type_t; + 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 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); + }; + + auto notdone = eve::true_(as()); + Z r = eve::nan(as()); + auto re = real(z); + if( eve::any(notdone) ) + { + notdone = next_interval(br_im_eq_0, notdone, is_real(z), r, z); + if( eve::any(notdone) ) + { + notdone = next_interval(br_re_gt_0, notdone, eve::is_gtz(re), r, z); + if( eve::any(notdone) ) + { + notdone = last_interval(br_re_lt_0, notdone, r, z); + } + } + } + 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)); + } + + //===------------------------------------------------------------------------------------------- + // bi + //===------------------------------------------------------------------------------------------- + template KYOSU_FORCEINLINE + auto bi(Z z) noexcept + { + using u_t = eve::underlying_type_t; + 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); + + }; + + 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)); + }; + + auto notdone = eve::true_(as()); + Z r = eve::nan(as()); + auto re = real(z); + if( eve::any(notdone) ) + { + notdone = next_interval(br_im_eq_0, notdone, is_real(z), r, z); + if( eve::any(notdone) ) + { + notdone = next_interval(br_re_gt_0, notdone, eve::is_gtz(re), r, z); + if( eve::any(notdone) ) + { + notdone = last_interval(br_re_lt_0, notdone, r, z); + } + } + } + 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)); + } + + //===------------------------------------------------------------------------------------------- + // aibi + //===------------------------------------------------------------------------------------------- + template KYOSU_FORCEINLINE + auto aibi(Z z) noexcept + { + using u_t = eve::underlying_type_t; + auto [sqzo3, zeta] = zet(z); + auto ip = cyl_bessel_i(eve::third(as()), zeta); + auto im = cyl_bessel_i(-eve::third(as()), zeta); + + return kumi::tuple{ invpi(as())*sqzo3*(im-ip), sqzo3*(ip+im)}; + } + + //===------------------------------------------------------------------------------------------- + // airy_ai + //===------------------------------------------------------------------------------------------- + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(kyosu::concepts::complex ) + { + return ai(z); + } + else + { + return cayley_extend(airy_ai, z); + } + } + + //===------------------------------------------------------------------------------------------- + // airy_bi + //===------------------------------------------------------------------------------------------- + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return bi(z); + } + else + { + return cayley_extend(airy_bi, z); + } + } + + //===------------------------------------------------------------------------------------------- + // airy + //===------------------------------------------------------------------------------------------- + template KYOSU_FORCEINLINE + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return aibi(z); + } + else + { + return cayley_extend2(airy, z); + } + } +} diff --git a/test/unit/function/airy_ai.cpp b/test/unit/function/airy_ai.cpp new file mode 100644 index 00000000..520a1a86 --- /dev/null +++ b/test/unit/function/airy_ai.cpp @@ -0,0 +1,85 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::airy_ai over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + using e_t = eve::element_type_t; + if constexpr(sizeof(e_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 reres{ + 1.6763717755962929e+07, 1.0671669307692060e+06, 5.1770135828865488e+04, 2.7209127311234188e+03 + , 1.7464957947925365e+02, 1.5043848907202149e+01, 1.9864751505771967e+00, 5.2597973488369920e-01 + , 2.7880648195500490e-01, 6.0458308371838042e-02, -1.0036228447036250e-01, 2.4722152686112525e-02 + , 8.1141306614663667e-03, -1.1116711994387154e-02, 6.2362124668067535e-03, -1.1845517847442273e-03 + }; + std::array imres{ + 5.7419608040318936e+07, -1.2409242416987466e+06, 2.6427659633621028e+04, -3.4216030566322274e+02 + , -3.0804431341868103e+01, 6.6862293791164378e+00, -1.3118365573923054e+00, 3.2173225694177510e-01 + , 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_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 new file mode 100644 index 00000000..92739d80 --- /dev/null +++ b/test/unit/function/airy_bi.cpp @@ -0,0 +1,62 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j0 over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + using e_t = eve::element_type_t; + if constexpr(sizeof(e_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 reres{-5.7419608040318385e+07, -1.2409242416987633e+06, -2.6427659633595144e+04, -3.4216029692014837e+02, 3.0804691728218572e+01, 6.6903996777286583e+00, 1.3523669469741970e+00, 5.5870307892357129e-01, 7.5248558508731545e-01, 7.1665807338276855e-01, -7.3369695701857451e-01, -1.3293492209488911e-01, 1.1081359835209212e+00, -2.0805776844472024e+00, 3.1132412017028535e+00, -3.5952187284708805e+00}; + + std::array imres{1.6763717755961375e+07, -1.0671669307691315e+06, 5.1770135827841987e+04, -2.7209127089468357e+03, 1.7464925781992241e+02, -1.5040759898113270e+01, 1.9670200948712224e+00, -4.4241768244532709e-01, 0.0000000000000000e+00, 6.1988929040084473e-01, -4.7545485913890873e-01, -1.2208933006131721e+00, -1.3901622231770903e+00, -1.5528897868070297e+00, -2.4569086817625085e+00, -5.1871698800009600e+00}; + + for(int i=0; i < 16; ++i) + { + 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_CASE_WITH ( "Check kyosu::airy_bi 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_bi; +// auto j0c = airy_bi(c); +// TTS_EQUAL(j0c, kyosu::conj(airy_bi(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_bi(z), kyosu::complex(airy_bi(T(0)))); +// } +// }; diff --git a/test/unit/function/cyl_bessel_k.cpp b/test/unit/function/cyl_bessel_k.cpp index cb9a6ef2..fe438422 100644 --- a/test/unit/function/cyl_bessel_k.cpp +++ b/test/unit/function/cyl_bessel_k.cpp @@ -9,6 +9,76 @@ #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" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) @@ -26,38 +96,23 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real" , -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; + constexpr int N = 0; + T v1 = N+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{-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} + vc_t{-7.6339268874816966e+01, -9.2192136865580512e+00, 1.6508508118018792e+01, 1.2375096096000567e+01 + , 2.3503760861380147e+00, -2.5566294609268230e+00, -2.4913144969412047e+00, -6.0742679340229055e-01 + , 1.5091129245821360e+00, 7.5263258886450732e-02, -1.0402938701852732e-01, -5.4882419994366216e-02 + , -6.5488350696814700e-03, 7.4931851019796377e-03, 4.9792792996237569e-03, 8.9451308743070568e-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} + 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); @@ -65,13 +120,117 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k over real" { auto c = kyosu::complex(re[k], im[k]); kyosu::cyl_bessel_k(v1, c, is_03); - auto vv = T(-0.3); + 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, kyosu::cyl_bessel_k(vv, c), 1.e-6); - vv = eve::dec(vv); + +// 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_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 = N+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{1.1900413043351881e-03, 1.5495899619039162e-03, 9.2313922815213754e-03, 1.4640077997088552e-02, -1.0136111628580130e-02, -1.0257686840586897e-01, -2.1965377035090466e-01, 1.3644035629075407e-02, 7.5455646229106821e-01, -1.5515848980450035e+00, -4.2927833713114696e+00, -4.0556204704598322e+00, 4.7439243146286634e+00, 2.2123510870256215e+01, 2.8376292912481247e+01, -1.8674126596346056e+01}, + }; + a_t imresi_03{ + 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); + 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(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_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_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 = N+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{-7.6339268874816966e+01, -9.2192136865580512e+00, 1.6508508118018792e+01, 1.2375096096000567e+01, 2.3503760861380147e+00, -2.5566294609268230e+00, -2.4913144969412047e+00, -6.0742679340229055e-01, 1.5091129245821360e+00, 7.5263258886450732e-02, -1.0402938701852732e-01, -5.4882419994366216e-02, -6.5488350696814700e-03, 7.4931851019796377e-03, 4.9792792996237569e-03, 8.9451308743070568e-04}, + }; + a_t imresi_03{ + 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); + 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(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_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); } } }