From fe69e00b2408c0539ff5e79c2133c0cd2be2008e Mon Sep 17 00:00:00 2001 From: jtlap Date: Fri, 3 Nov 2023 15:45:57 +0100 Subject: [PATCH] More Bessel functions rework --- include/kyosu/functions.hpp | 7 + include/kyosu/functions/cyl_bessel_h1n.hpp | 78 +++++++ include/kyosu/functions/cyl_bessel_h2n.hpp | 78 +++++++ include/kyosu/functions/cyl_bessel_i0.hpp | 9 +- include/kyosu/functions/cyl_bessel_i1.hpp | 9 +- include/kyosu/functions/cyl_bessel_in.hpp | 16 +- include/kyosu/functions/cyl_bessel_j0.hpp | 6 +- include/kyosu/functions/cyl_bessel_j1.hpp | 5 +- include/kyosu/functions/cyl_bessel_jn.hpp | 14 +- include/kyosu/functions/cyl_bessel_k0.hpp | 74 +++++++ include/kyosu/functions/cyl_bessel_k1.hpp | 73 +++++++ include/kyosu/functions/cyl_bessel_kn.hpp | 75 +++++++ include/kyosu/functions/cyl_bessel_y0.hpp | 13 +- include/kyosu/functions/cyl_bessel_y1.hpp | 77 +++++++ include/kyosu/functions/cyl_bessel_yn.hpp | 75 +++++++ include/kyosu/functions/if_else.hpp | 5 + include/kyosu/types/impl/bessel.hpp | 6 + .../kyosu/types/impl/detail/bessel_utils.hpp | 202 ++++++++++++++++++ include/kyosu/types/impl/detail/besselhn.hpp | 46 ++++ include/kyosu/types/impl/detail/besselin.hpp | 4 +- include/kyosu/types/impl/detail/besseljn.hpp | 15 +- include/kyosu/types/impl/detail/besselkn.hpp | 53 +++++ include/kyosu/types/impl/detail/bessely0.hpp | 38 ++++ include/kyosu/types/impl/detail/bessely1.hpp | 28 +++ include/kyosu/types/impl/detail/besselyn.hpp | 40 ++++ .../kyosu/types/impl/detail/cyl_bessely1.hpp | 28 +++ include/kyosu/types/impl/detail/li0.hpp | 29 +++ test/doc/besselin.hpp | 33 --- test/doc/cyl_bessel_h1n.cpp | 17 ++ test/doc/cyl_bessel_h2n.cpp | 17 ++ test/doc/cyl_bessel_j0.cpp | 8 + test/doc/cyl_bessel_kn.cpp | 17 ++ test/doc/cyl_bessel_y0.cpp | 33 +++ test/doc/cyl_bessel_y1.cpp | 33 +++ test/doc/cyl_bessel_yn.cpp | 34 +++ test/unit/function/cyl_bessel_i0.cpp | 6 +- test/unit/function/cyl_bessel_i1.cpp | 6 +- test/unit/function/cyl_bessel_in.cpp | 23 +- test/unit/function/cyl_bessel_j0.cpp | 6 +- test/unit/function/cyl_bessel_j1.cpp | 6 +- test/unit/function/cyl_bessel_jn.cpp | 27 ++- test/unit/function/cyl_bessel_k0.cpp | 56 +++++ test/unit/function/cyl_bessel_k1.cpp | 57 +++++ test/unit/function/cyl_bessel_kn.cpp | 64 ++++++ test/unit/function/cyl_bessel_y0.cpp | 59 +++++ test/unit/function/cyl_bessel_y1.cpp | 58 +++++ test/unit/function/cyl_bessel_yn.cpp | 64 ++++++ 47 files changed, 1613 insertions(+), 114 deletions(-) create mode 100644 include/kyosu/functions/cyl_bessel_h1n.hpp create mode 100644 include/kyosu/functions/cyl_bessel_h2n.hpp create mode 100644 include/kyosu/functions/cyl_bessel_k0.hpp create mode 100644 include/kyosu/functions/cyl_bessel_k1.hpp create mode 100644 include/kyosu/functions/cyl_bessel_kn.hpp create mode 100644 include/kyosu/functions/cyl_bessel_y1.hpp create mode 100644 include/kyosu/functions/cyl_bessel_yn.hpp create mode 100644 include/kyosu/types/impl/detail/bessel_utils.hpp create mode 100644 include/kyosu/types/impl/detail/besselhn.hpp create mode 100644 include/kyosu/types/impl/detail/besselkn.hpp create mode 100644 include/kyosu/types/impl/detail/bessely0.hpp create mode 100644 include/kyosu/types/impl/detail/bessely1.hpp create mode 100644 include/kyosu/types/impl/detail/besselyn.hpp create mode 100644 include/kyosu/types/impl/detail/cyl_bessely1.hpp create mode 100644 include/kyosu/types/impl/detail/li0.hpp delete mode 100644 test/doc/besselin.hpp create mode 100644 test/doc/cyl_bessel_h1n.cpp create mode 100644 test/doc/cyl_bessel_h2n.cpp create mode 100644 test/doc/cyl_bessel_kn.cpp create mode 100644 test/doc/cyl_bessel_y0.cpp create mode 100644 test/doc/cyl_bessel_y1.cpp create mode 100644 test/doc/cyl_bessel_yn.cpp create mode 100644 test/unit/function/cyl_bessel_k0.cpp create mode 100644 test/unit/function/cyl_bessel_k1.cpp create mode 100644 test/unit/function/cyl_bessel_kn.cpp create mode 100644 test/unit/function/cyl_bessel_y0.cpp create mode 100644 test/unit/function/cyl_bessel_y1.cpp create mode 100644 test/unit/function/cyl_bessel_yn.cpp diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index 70bdc83c..92af6a8f 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -54,6 +54,13 @@ #include #include #include +#include +#include +#include +#include +#include +#include +#include #include #include #include diff --git a/include/kyosu/functions/cyl_bessel_h1n.hpp b/include/kyosu/functions/cyl_bessel_h1n.hpp new file mode 100644 index 00000000..70dea117 --- /dev/null +++ b/include/kyosu/functions/cyl_bessel_h1n.hpp @@ -0,0 +1,78 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_bessel_h1n: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_h1n; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_h1n_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& z) noexcept + { + using e_t = eve::element_type_t; + return complex(cyl_bessel_jn(n, z), cyl_bessel_yn(n, z)); + } + + template + KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target0, target1)) + { + return eve::tag_invoke(*this, target0, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_bessel_h1n +//! @brief Computes the Bessel/Hankel functions of the third kind, +//! \f$ H_n^{(1)}(z) = J_n(z)+iY_n(z)\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_h1n(int n, T z) noexcept; +//! template constexpr auto cyl_bessel_h1n(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$H_n^{(1)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_h1n.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_h1n cyl_bessel_h1n = {}; +} diff --git a/include/kyosu/functions/cyl_bessel_h2n.hpp b/include/kyosu/functions/cyl_bessel_h2n.hpp new file mode 100644 index 00000000..d0f04f0d --- /dev/null +++ b/include/kyosu/functions/cyl_bessel_h2n.hpp @@ -0,0 +1,78 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_bessel_h2n: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_h1n; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_h2n_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& z) noexcept + { + using e_t = eve::element_type_t; + return complex(cyl_bessel_jn(n, z), -cyl_bessel_yn(n, z)); + } + + template + KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target0, target1)) + { + return eve::tag_invoke(*this, target0, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_bessel_h2n +//! @brief Computes the Bessel/Hankel functions of the third kind , +//! \f$ H_n^{(2)} = J_n(z)-iY_n(z)\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_h2n(int n, T z) noexcept; +//! template constexpr auto cyl_bessel_h2n(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return \f$H_n^{(2)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_h2n.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_h2n cyl_bessel_h2n = {}; +} diff --git a/include/kyosu/functions/cyl_bessel_i0.hpp b/include/kyosu/functions/cyl_bessel_i0.hpp index cab13917..a5f1a21c 100644 --- a/include/kyosu/functions/cyl_bessel_i0.hpp +++ b/include/kyosu/functions/cyl_bessel_i0.hpp @@ -39,11 +39,10 @@ namespace kyosu //! @addtogroup functions //! @{ //! @var cyl_bessel_i0 -//! @brief Computes the Bessel function of the first kind, -//! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) -//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. +//! @brief Computes the modified Bessel function of the first kind \f$I_{0}(x)=J_{0}(ix)\f$ +//! extended to the complex plane and cayley_dickson algebras. //! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! It is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. //! //! @code //! #include @@ -65,7 +64,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_0(z)\f$. +//! * returns \f$I_0(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_i1.hpp b/include/kyosu/functions/cyl_bessel_i1.hpp index dbc8870e..ecbfe37c 100644 --- a/include/kyosu/functions/cyl_bessel_i1.hpp +++ b/include/kyosu/functions/cyl_bessel_i1.hpp @@ -39,11 +39,10 @@ namespace kyosu //! @addtogroup functions //! @{ //! @var cyl_bessel_i1 -//! @brief Computes the Bessel function of the first kind, -//! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) -//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. +//! @brief Computes the modified Bessel function of the first kind, +//! \f$ I_1(x)= _iJ_1(ix) \f$ extended to the complex plane and cayley_dickson algebras. //! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! It is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 0\f$. //! //! @code //! #include @@ -65,7 +64,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_0(z)\f$. +//! * returns \f$I_1(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_in.hpp b/include/kyosu/functions/cyl_bessel_in.hpp index 0002b675..7bef2fcd 100644 --- a/include/kyosu/functions/cyl_bessel_in.hpp +++ b/include/kyosu/functions/cyl_bessel_in.hpp @@ -39,11 +39,11 @@ namespace kyosu //====================================================================================================================== //! @addtogroup functions //! @{ -//! @brief Computes the Bessel functions of the first kind, -//! \f$ J_{n}(x)=\sum_{p=0}^{\infty}{\frac{(-1)^p}{p!\,\Gamma (p+n +1)}} -//! {\left({x \over 2}\right)}^{2p+n }\f$. +//! @var cyl_bessel_in +//! @brief Computes the modified Bessel functions of the first kind \f$I_{n}(x)=i^{-n}J_{n }(ix)\f$, +//! extended to the complex plane and cayley_dickson algebras. //! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which +//! It is the solution of \f$ x^{2}y''+xy'+(x^2+n^2)y=0\f$ for which //! \f$ y(0) = 0\f$ if \f$n \ne 0\f$ else \f$1\f$. //! //! @code @@ -55,8 +55,8 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto cyl_bessel_in(int n, T z) noexcept; -//! template constexpr T cyl_bessel_in(N n, T z) noexcept; +//! template constexpr auto cyl_bessel_in(int n, T z) noexcept; +//! template constexpr T cyl_bessel_in(int n, T z) noexcept; //! } //! @endcode //! @@ -66,9 +66,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_n(z)\f$. -//! -//! @warning Up to now the cayley_dickson versions have only been imlemented dor scalar int values of n. +//! * returns \f$J_n(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_j0.hpp b/include/kyosu/functions/cyl_bessel_j0.hpp index abb35ae6..c15741cd 100644 --- a/include/kyosu/functions/cyl_bessel_j0.hpp +++ b/include/kyosu/functions/cyl_bessel_j0.hpp @@ -41,9 +41,9 @@ namespace kyosu //! @var cyl_bessel_j0 //! @brief Computes the Bessel function of the first kind, //! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) -//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. +//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson algebras. //! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! It is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. //! //! @code //! #include @@ -65,7 +65,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_0(z)\f$. +//! * returns \f$J_0(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_j1.hpp b/include/kyosu/functions/cyl_bessel_j1.hpp index b6442abb..36613d2f 100644 --- a/include/kyosu/functions/cyl_bessel_j1.hpp +++ b/include/kyosu/functions/cyl_bessel_j1.hpp @@ -42,8 +42,7 @@ namespace kyosu //! @brief Computes the Bessel function of the first kind, //! \f$ J_1(x)=\frac1{\pi }\int _{0}^{\pi}\cos(\tau-x\sin \tau )\,\mathrm {d} \tau \f$ //! extended to the complex plane and cayley_dickson values. -//! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! It is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. //! //! @code //! #include @@ -65,7 +64,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_0(z)\f$. +//! * returns \f$J_0(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_jn.hpp b/include/kyosu/functions/cyl_bessel_jn.hpp index f225676d..09e3ffc0 100644 --- a/include/kyosu/functions/cyl_bessel_jn.hpp +++ b/include/kyosu/functions/cyl_bessel_jn.hpp @@ -39,11 +39,13 @@ namespace kyosu //====================================================================================================================== //! @addtogroup functions //! @{ +//! @var cyl_bessel_jn //! @brief Computes the Bessel functions of the first kind, //! \f$ J_{n}(x)=\sum_{p=0}^{\infty}{\frac{(-1)^p}{p!\,\Gamma (p+n +1)}} -//! {\left({x \over 2}\right)}^{2p+n }\f$. +//! {\left({x \over 2}\right)}^{2p+n }\f$ +//! extended to the complex plane and cayley_dickson values. //! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which +//! It is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which //! \f$ y(0) = 0\f$ if \f$n \ne 0\f$ else \f$1\f$. //! //! @code @@ -55,8 +57,8 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto cyl_bessel_jn(int n, T z) noexcept; -//! template constexpr T cyl_bessel_jn(N n, T z) noexcept; +//! template constexpr auto cyl_bessel_jn(int n, T z) noexcept; +//! template constexpr T cyl_bessel_jn(int n, T z) noexcept; //! } //! @endcode //! @@ -66,9 +68,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_n(z)\f$. -//! -//! @warning Up to now the cayley_dickson versions have only been imlemented dor scalar int values of n. +//! * returns \f$J_n(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_k0.hpp b/include/kyosu/functions/cyl_bessel_k0.hpp new file mode 100644 index 00000000..a7f38da3 --- /dev/null +++ b/include/kyosu/functions/cyl_bessel_k0.hpp @@ -0,0 +1,74 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_bessel_k0: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_k0; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_k0_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::cyl_bessel_k0(v); } + + template + KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) + { + return eve::tag_invoke(*this, target); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_bessel_k0 +//! @brief Computes the modified Bessel function of the second kind, +//! \f$ K_0(x)=\lim_{\alpha\to 0}{\frac {\pi }{2}}{\frac {I_{-\alpha }(x)-I_{\alpha }(x)}{\sin \alpha \pi }}\f$. +//! extended to the complex plane and cayley_dickson values. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_k0(T z) noexcept; +//! template constexpr auto cyl_bessel_k0(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$K_0(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_k0.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_k0 cyl_bessel_k0 = {}; +} diff --git a/include/kyosu/functions/cyl_bessel_k1.hpp b/include/kyosu/functions/cyl_bessel_k1.hpp new file mode 100644 index 00000000..d5c57d8e --- /dev/null +++ b/include/kyosu/functions/cyl_bessel_k1.hpp @@ -0,0 +1,73 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_bessel_k1: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_k1; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_k1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::cyl_bessel_k1(v); } + + template + KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) + { + return eve::tag_invoke(*this, target); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_bessel_k1 +//! @brief Computes the Bessel function of the second kind, +//! \f$ K_1(x)\lim_{\alpha\to 1}{\frac {\pi }{2}}{\frac {I_{-\alpha }(x)-I_{\alpha }(x)}{\sin \alpha \pi }}\f$ extended to the complex plane and cayley_dickson values. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_k1(T z) noexcept; +//! template constexpr auto cyl_bessel_k1(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$K_1(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_k1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_k1 cyl_bessel_k1 = {}; +} diff --git a/include/kyosu/functions/cyl_bessel_kn.hpp b/include/kyosu/functions/cyl_bessel_kn.hpp new file mode 100644 index 00000000..3283b61c --- /dev/null +++ b/include/kyosu/functions/cyl_bessel_kn.hpp @@ -0,0 +1,75 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_bessel_kn: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_kn; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_kn_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& v) noexcept { return eve::cyl_bessel_kn(n, v); } + + template + KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target0, target1)) + { + return eve::tag_invoke(*this, target0, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_bessel_kn +//! @brief Computes the modified Bessel functions of the second kind, +//! \f$ K_{n}(x)=\lim_{\alpha\to n}{\frac {\pi }{2}}{\frac {I_{-\alpha }(x)-I_{\alpha }(x)}{\sin \alpha \pi }}\f$. +//! extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_kn(int n, T z) noexcept; +//! template constexpr T cyl_bessel_kn(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$K_n(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_kn.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_kn cyl_bessel_kn = {}; +} diff --git a/include/kyosu/functions/cyl_bessel_y0.hpp b/include/kyosu/functions/cyl_bessel_y0.hpp index b597c144..1511d819 100644 --- a/include/kyosu/functions/cyl_bessel_y0.hpp +++ b/include/kyosu/functions/cyl_bessel_y0.hpp @@ -19,7 +19,10 @@ namespace kyosu::tags KYOSU_DEFERS_CALLABLE(cyl_bessel_y0_); template - static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::cyl_bessel_y0(v); } + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept + { + return eve::cyl_bessel_y0(v); + } template KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) @@ -40,10 +43,8 @@ namespace kyosu //! @{ //! @var cyl_bessel_y0 //! @brief Computes the Bessel function of the second kind, -//! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) -//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. -//! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! \f$ Y_0(x)=\lim_{\alpha\to 0}{{\frac {J_{\alpha }(x)\cos(\alpha\pi)-J_{-\alpha }(x)}{\sin(\alpha\pi)}}}\f$, +//! extended to the complex plane and cayley_dickson algebras. //! //! @code //! #include @@ -65,7 +66,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_0(z)\f$. +//! * returns \f$Y_0(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_y1.hpp b/include/kyosu/functions/cyl_bessel_y1.hpp new file mode 100644 index 00000000..dc179389 --- /dev/null +++ b/include/kyosu/functions/cyl_bessel_y1.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_cyl_bessel_y1: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_y1; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_y1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept + { + return eve::cyl_bessel_y1(v); + } + + template + KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) + { + return eve::tag_invoke(*this, target); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_bessel_y1 +//! @brief Computes the Bessel function of the second kind, +//! \f$ Y_1(x)=\lim_{\alpha\to 1}{{\frac {J_{\alpha }(x)\cos(\alpha\pi)-J_{-\alpha }(x)}{\sin(\alpha\pi)}}}\f$, +//! extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_y1(T z) noexcept; +//! template constexpr T cyl_bessel_y1(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$Y_1(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_y1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_y1 cyl_bessel_y1 = {}; +} diff --git a/include/kyosu/functions/cyl_bessel_yn.hpp b/include/kyosu/functions/cyl_bessel_yn.hpp new file mode 100644 index 00000000..075be242 --- /dev/null +++ b/include/kyosu/functions/cyl_bessel_yn.hpp @@ -0,0 +1,75 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_bessel_yn: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_yn; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_yn_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& v) noexcept { return eve::cyl_bessel_yn(n, v); } + + template + KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target0, target1)) + { + return eve::tag_invoke(*this, target0, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_bessel_yn +//! @brief Computes the modified Bessel functions of the second kind, +//! \f$ Y_n(x)=\lim_{\alpha\to n}{{\frac {J_{\alpha }(x)\cos(\alpha\pi)-J_{-\alpha }(x)}{\sin(\alpha\pi)}}}\f$, +//! extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_yn(int n, T z) noexcept; +//! template constexpr T cyl_bessel_yn(N n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$Y_n(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_yn.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_yn cyl_bessel_yn = {}; +} diff --git a/include/kyosu/functions/if_else.hpp b/include/kyosu/functions/if_else.hpp index 3361cdda..7708eb05 100644 --- a/include/kyosu/functions/if_else.hpp +++ b/include/kyosu/functions/if_else.hpp @@ -15,6 +15,11 @@ namespace kyosu::tags { using callable_tag_type = callable_if_else; + KYOSU_DEFERS_CALLABLE(if_else_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, auto c, T const& t, U u) noexcept { return eve::if_else(c, t, u); } + KYOSU_FORCEINLINE auto operator()(auto const& m, auto const& t, auto const& f) const noexcept -> decltype(eve::tag_invoke(*this,m,t,f)) { diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp index 71257611..01cb5975 100644 --- a/include/kyosu/types/impl/bessel.hpp +++ b/include/kyosu/types/impl/bessel.hpp @@ -12,9 +12,15 @@ #include #include +#include #include #include #include #include #include #include +#include +#include +#include +#include +#include diff --git a/include/kyosu/types/impl/detail/bessel_utils.hpp b/include/kyosu/types/impl/detail/bessel_utils.hpp new file mode 100644 index 00000000..6798bdef --- /dev/null +++ b/include/kyosu/types/impl/detail/bessel_utils.hpp @@ -0,0 +1,202 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // bound + //===------------------------------------------------------------------------------------------- + template + auto bound(Z const & z) noexcept + { + using u_t = eve::underlying_type_t; + return int(eve::ceil(eve::maximum(abs(z)*u_t(1.1)+2))); + } + + //===------------------------------------------------------------------------------------------- + // R use lentz + //===------------------------------------------------------------------------------------------- + template + struct fraction_traits_simple + { + using value_type = typename Gen::result_type; + using result_type = typename Gen::result_type; + using scalar_type = typename Gen::scalar_type; + + static result_type a(const value_type&) noexcept + { + return result_type(1); + } + static result_type b(const value_type& v) noexcept + { + return v; + } + + }; + + template + struct fraction_traits_pair + { + using value_type = typename Gen::value_type; + using result_type = typename kumi::element<0, value_type>; + using scalar_type = typename Gen::scalar_type; + + static result_type a(const value_type& v) noexcept + { + return kumi::get<1>(v); + } + static result_type b(const value_type& v) noexcept + { + return kumi::get<0>(v); + } + }; + + template + struct fraction_traits + : std::conditional_t< + std::same_as, + fraction_traits_simple, + fraction_traits_pair> + { + }; + + template + inline auto lentz_b(Gen& g, const U& eps, size_t const& max_terms = 1000) noexcept + { + using traits = fraction_traits; + using value_type = typename traits::value_type; + using result_type = typename traits::result_type; + using scalar_type = typename traits::scalar_type; + using u_t = scalar_type; + + u_t tiny = 16*eve::smallestposval(eve::as()) ; + u_t terminator(eps); + constexpr auto b_only = std::same_as; + + auto v = g(); + + result_type f, C, delta; + f = traits::b(v); + C = f; + result_type D{}; + + std::uintmax_t counter(max_terms); + do{ + v = g(); + if constexpr(b_only) + { + D = v + D; + C = v + kyosu::rec(C); + } + else + { + D = kumi::get<0>(v) + kumi::get<1>(v) * D; + C = kumi::get<0>(v) + kumi::get<1>(v) / C; + } + D = kyosu::if_else(kyosu::is_eqz(D), tiny, D); + C = kyosu::if_else(kyosu::is_eqz(C), tiny, C); + D = kyosu::rec(D); + delta = C*D; + f = f * delta; + } while(eve::any((kyosu::abs(kyosu::dec(delta)) > terminator) && --counter)); + g.reset(); + return f; + } + + template + struct R_estimate + { + using result_type = T; + using value_type = T; + using scalar_type = eve::underlying_type_t; + + int n; + T rz; + int i; + int sgn; + + R_estimate(size_t n_, T z_) + : n(n_), rz(kyosu::rec(z_)) + { + reset(); + } + + void reset() + { + i = 0; sgn = -1; + } + + value_type operator()() + { + sgn = -sgn; + ++i; + auto r = sgn*2*(n+i-1)*rz; + return r; + } + }; + + template + inline auto R(size_t n, Z z) noexcept + // compute the ratio Jn(n, z)/Jn(n-1, z)J + { + using u_t = eve::underlying_type_t; + R_estimate re(n, z); + auto r = lentz_b(re, eve::eps(eve::as())); + return r; + } + + template + inline auto Rs(size_t n, Z z) noexcept + // compute the ratio Jn(i, z)/Jn(i+1, z)J for i = n-1:-1:0 + { + std::vector rs(n); + using u_t = eve::underlying_type_t; + if (is_eqz(n)) return rs; + Z r{}, s{}; + auto rz = rec(z); + rs[n-1] = R(n, z); + for(int i=n-1; i >= 1; --i) + { + rs[i-1] = 2*i*rz-rec(rs[i]); + } + return rs; + } + + + template + inline auto Js(size_t n, Z z) noexcept + // compute Jn(i, z) for i = n:-1:0 + { + auto rs = Rs(n+1, z); + std::vector js(n+1); + js[n] = cyl_bessel_jn(n, z); + for(int i=n-1; i >= 0; --i) + { + js[i] = js[i+1]*rs[i]; + } + return js; + } + + + + template + inline auto arg_adjust(Z z) noexcept + // modify ipart of z modulo 2*pi fit in ]-pi +pi] + { + if constexpr(kyosu::concepts::complex) + { + ipart(z) = eve::to_nearest(eve::rem)(ipart(z), 2*eve::pi(eve::as(ipart(z)))); + return z; + } + else return z; + } + +} diff --git a/include/kyosu/types/impl/detail/besselhn.hpp b/include/kyosu/types/impl/detail/besselhn.hpp new file mode 100644 index 00000000..95ddd8be --- /dev/null +++ b/include/kyosu/types/impl/detail/besselhn.hpp @@ -0,0 +1,46 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_h1n + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, int n, Z z) noexcept + { + if constexpr(concepts::complex) + { + auto muli = [](auto z){ auto [r, i] = z; return complex(-i, r); }; + return cyl_bessel_jn(n, z)+muli(cyl_bessel_yn(n, z)); + } + else + { + return cayley_extend_rev(cyl_bessel_h1n, n, z); + } + } + + template + auto dispatch(eve::tag_of, int n, Z z) noexcept + { + if constexpr(concepts::complex) + { + auto muli = [](auto z){ auto [r, i] = z; return complex(-i, r); }; + return cyl_bessel_jn(n, z)-muli(cyl_bessel_yn(n, z)); + } + else + { + return cayley_extend_rev(cyl_bessel_h2n, n, z); + } + } +} diff --git a/include/kyosu/types/impl/detail/besselin.hpp b/include/kyosu/types/impl/detail/besselin.hpp index a1659e8b..3e33b171 100644 --- a/include/kyosu/types/impl/detail/besselin.hpp +++ b/include/kyosu/types/impl/detail/besselin.hpp @@ -26,8 +26,8 @@ namespace kyosu::_ else if (n%4 == 2) return complex(eve::mone(eve::as())); else return complex(eve::zero(eve::as()), eve::one(eve::as())); }; - - return miton(n)*cyl_bessel_jn(n,complex(-ipart(z), real(z))); + auto an = eve::abs(n); + return miton(an)*cyl_bessel_jn(an,complex(-ipart(z), real(z))); } else { diff --git a/include/kyosu/types/impl/detail/besseljn.hpp b/include/kyosu/types/impl/detail/besseljn.hpp index 21a96158..f976dcbf 100644 --- a/include/kyosu/types/impl/detail/besseljn.hpp +++ b/include/kyosu/types/impl/detail/besseljn.hpp @@ -37,8 +37,7 @@ namespace kyosu::_ { using e_t = as_real_t; using u_t = eve::underlying_type_t; - auto n = u_t(nn); - auto an = eve::abs(n); + auto n = u_t(eve::abs(nn)); auto az = kyosu::abs(z); auto forward = [n](auto z){ @@ -66,7 +65,7 @@ namespace kyosu::_ // Starting point for backward recurrence // for when |Jn(x)|~10e-mg // using the secant method. - auto n0 = inc(eve::nearest( u_t(1.1)*az)); + auto n0 = inc(eve::ceil( u_t(1.1)*az)); auto f0 = minus_log10_cyl_j_at_infinity(n0, az) - mg; auto n1 = n0 + 5; auto f1 = minus_log10_cyl_j_at_infinity(n1, az) - mg; @@ -94,7 +93,7 @@ namespace kyosu::_ auto ejn = minus_log10_cyl_j_at_infinity(n, az); auto t = ejn <= hmp; auto obj = eve::if_else(t, sd, hmp+ejn); - auto n0 = eve::if_else(t, eve::nearest(e_t(1.1)*az), n); + auto n0 = eve::if_else(t, eve::ceil(e_t(1.1)*az), n); auto f0 = minus_log10_cyl_j_at_infinity(n0, az) - obj; auto n1 = n0 + 5; auto f1 = minus_log10_cyl_j_at_infinity(n1, az) - obj; @@ -115,8 +114,9 @@ namespace kyosu::_ }; auto backward = [az, n, ini_for_br_1, ini_for_br_2](auto z){ - auto m = ini_for_br_1(az, e_t(200)); - m = eve::if_else ( m >= n, ini_for_br_2(n, az, e_t(15)), m); + auto m1 = ini_for_br_1(az, e_t(200)); + auto m2 = ini_for_br_2(n, az, e_t(15)); + auto m = eve::if_else( m1 >= n && eve::is_not_nan(m2), m2, m1); auto cf2 = Z(0); auto cf1 = complex(eve::sqrtsmallestposval(eve::as< e_t>())); Z cf(cf2); @@ -153,7 +153,8 @@ namespace kyosu::_ } } auto sgnaltern = [n](auto x){return eve::if_else(eve::is_ltz(x), eve::one, eve::sign_alternate(n));}; - return sgnaltern(srz)*sgnaltern(n)*r; + r = sgnaltern(srz)*sgnaltern(n)*r; + return nn < 0 ? r*eve::sign_alternate(u_t(nn)) : r; } } else diff --git a/include/kyosu/types/impl/detail/besselkn.hpp b/include/kyosu/types/impl/detail/besselkn.hpp new file mode 100644 index 00000000..f9bf7b4c --- /dev/null +++ b/include/kyosu/types/impl/detail/besselkn.hpp @@ -0,0 +1,53 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_kn + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, int n, Z z) noexcept + { + if constexpr(concepts::complex) + { + using u_t = eve::underlying_type_t; + auto argz = arg(z); + n = eve::abs(n); + auto piotwo = eve::pio_2(eve::as()); + auto i = complex(u_t(0), u_t(1)); + auto cpi = piotwo*i*exp_ipi(u_t(n)/2); + auto cmi = -piotwo*i*exp_ipi(-u_t(n)/2); + auto epio2 = exp_ipi(eve::half(eve::as())); + auto empio2 = exp_ipi(eve::mhalf(eve::as())); + auto r = if_else(eve::is_ltz(argz) + , cpi*cyl_bessel_h1n(n, z*epio2) + , cmi*cyl_bessel_h2n(n, z*empio2)); + return if_else(is_eqz(z), complex(eve::inf(eve::as())), r); + } + else + { + return cayley_extend_rev(cyl_bessel_kn, n, z); + } + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return cyl_bessel_kn(0, z); + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return cyl_bessel_kn(1, z); + } +} diff --git a/include/kyosu/types/impl/detail/bessely0.hpp b/include/kyosu/types/impl/detail/bessely0.hpp new file mode 100644 index 00000000..17ef8ad4 --- /dev/null +++ b/include/kyosu/types/impl/detail/bessely0.hpp @@ -0,0 +1,38 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_y0 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + using u_t = eve::underlying_type_t; + auto twoopi = eve::two_o_pi(eve::as()); + auto egamma = eve::egamma(eve::as()); + auto eps = eve::eps(eve::as()); + auto j0z = cyl_bessel_j0(z); + auto bd = bound(z); + auto js = Js(2*bd, z); + Z s{}, sk{}; + auto sgn = -rec(j0z); + int k = 1; + do { + sk = sgn*js[2*k]/k; + ++k; + sgn = -sgn; + s+= sk; + } while (k < bd && eve::any(kyosu::abs(sk) > abs(s)*eps)); + return if_else(is_eqz(z), complex(eve::minf(eve::as())), twoopi*((log(z/2)+egamma)-2*s)*js[0]); //cyl_bessel_j0(z)); + } +} diff --git a/include/kyosu/types/impl/detail/bessely1.hpp b/include/kyosu/types/impl/detail/bessely1.hpp new file mode 100644 index 00000000..d2705a28 --- /dev/null +++ b/include/kyosu/types/impl/detail/bessely1.hpp @@ -0,0 +1,28 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_y1 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + using u_t = eve::underlying_type_t; + auto twoopi = eve::two_o_pi(eve::as()); + auto r1 = R(1, z); + auto y0 = cyl_bessel_y0(z); + auto recs1 = rec(r1)-twoopi/(z*cyl_bessel_j0(z)*y0); + return if_else(is_eqz(z), complex(eve::minf(eve::as())), y0*recs1); + } +} diff --git a/include/kyosu/types/impl/detail/besselyn.hpp b/include/kyosu/types/impl/detail/besselyn.hpp new file mode 100644 index 00000000..7b484b4b --- /dev/null +++ b/include/kyosu/types/impl/detail/besselyn.hpp @@ -0,0 +1,40 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_yn + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, int nn, Z z) noexcept + { + using u_t = eve::underlying_type_t; + auto y = cyl_bessel_y0(z); + if (nn != 0) + { + auto n = eve::abs(nn); + auto js = Js(n, z); + using u_t = eve::underlying_type_t; + auto twoopi = eve::two_o_pi(eve::as()); + for(int i=1; i <= n ; ++i) + { + y = (js[i]*y-twoopi*rec(z))/js[i-1]; + } + auto r = if_else(eve::is_gtz(real(z)) && is_real(z) && is_nan(y), complex(real(y)), y); + r = if_else(is_eqz(z), complex(eve::minf(eve::as())), r); + return nn < 0 ? r*eve::sign_alternate(u_t(n)) : r; + } + else return y; + } +} diff --git a/include/kyosu/types/impl/detail/cyl_bessely1.hpp b/include/kyosu/types/impl/detail/cyl_bessely1.hpp new file mode 100644 index 00000000..6b5b023a --- /dev/null +++ b/include/kyosu/types/impl/detail/cyl_bessely1.hpp @@ -0,0 +1,28 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_y1 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + using u_t = eve::underlying_type_t; + auto twoopi = eve::two_o_pi(eve::as()); + auto r1 = R(1, z); + auto y0 = cyl_bessel_y0(z); + auto recs1 = rec(r1)-twoopi/(z*cyl_bessel_j0(z)*y0) + + return y0*recs1; + } +} diff --git a/include/kyosu/types/impl/detail/li0.hpp b/include/kyosu/types/impl/detail/li0.hpp new file mode 100644 index 00000000..fea51dad --- /dev/null +++ b/include/kyosu/types/impl/detail/li0.hpp @@ -0,0 +1,29 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_lbessel_i0 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return cyl_lbessel_in(0, z); + } + else + { + return cayley_extend(cyl_lbessel_i0, z); + } + } +} diff --git a/test/doc/besselin.hpp b/test/doc/besselin.hpp deleted file mode 100644 index 3e05cf88..00000000 --- a/test/doc/besselin.hpp +++ /dev/null @@ -1,33 +0,0 @@ -#include -#include -#include - -int main() -{ - using kyosu::cyl_bessel_j0; - std::cout.precision(16); - std::cout << std::scientific << std::endl; - std::cout << "1.0 " << " -> " << cyl_bessel_j0(1.0) << "\n"; - std::cout << "1+0i " << " -> " << cyl_bessel_j0(kyosu::complex(1.0, 0.0)) << "\n"; - std::cout << "15.0 " << " -> " << cyl_bessel_j0(15.0) << "\n"; - std::cout << "15+0i " << " -> " << cyl_bessel_j0(kyosu::complex(15.0, 0.0)) << "\n"; - std::cout << "40.0 " << " -> " << cyl_bessel_j0(40.0) << "\n"; - std::cout << "40+0i " << " -> " << cyl_bessel_j0(kyosu::complex(40.0, 0.0)) << "\n"; - std::cout << "60.0 " << " -> " << cyl_bessel_j0(60.0) << "\n"; - std::cout << "60+0i " << " -> " << cyl_bessel_j0(kyosu::complex(60.0, 0.0)) << "\n"; - - eve::wide> z(1.0, 15.0, 40.0, 60.0); - auto zz = kyosu::complex(z); - std::cout << z << " -> \n" << cyl_bessel_j0(z) << "\n"; - std::cout << zz << " -> \n" << cyl_bessel_j0(zz) << "\n"; - - std::cout << "1.0 " << " -> " << cyl_bessel_j0(1.0) << "\n"; - std::cout << "1.0+36.0i " << " -> " << cyl_bessel_j0(kyosu::complex(1.0, 36.0)) << "\n"; - std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_bessel_j0(kyosu::quaternion(1.0, 36.0, 2.0, 1.5)) << "\n"; - - - eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); - auto z0 = kyosu::complex(z1, z2); - std::cout << z0 << "\n -> " << cyl_bessel_j0(z0) << "\n"; - return 0; -} diff --git a/test/doc/cyl_bessel_h1n.cpp b/test/doc/cyl_bessel_h1n.cpp new file mode 100644 index 00000000..03a21143 --- /dev/null +++ b/test/doc/cyl_bessel_h1n.cpp @@ -0,0 +1,17 @@ +#include +#include +#include + +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)); + for(int n=2; n <= 2; ++n) + { + auto h1n = kyosu::cyl_bessel_h1n(n, z); + std::cout << "n " << n << " z " << z << std::endl; + std::cout << "yl_bessel_h1n(n, z) " << h1n << std::endl; + } + return 0; +} diff --git a/test/doc/cyl_bessel_h2n.cpp b/test/doc/cyl_bessel_h2n.cpp new file mode 100644 index 00000000..b6bf7701 --- /dev/null +++ b/test/doc/cyl_bessel_h2n.cpp @@ -0,0 +1,17 @@ +#include +#include +#include + +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)); + for(int n=2; n <= 2; ++n) + { + auto h1n = kyosu::cyl_bessel_h2n(n, z); + std::cout << "n " << n << " z " << z << std::endl; + std::cout << "yl_bessel_h2n(n, z) " << h1n << std::endl; + } + return 0; +} diff --git a/test/doc/cyl_bessel_j0.cpp b/test/doc/cyl_bessel_j0.cpp index c6fc42f9..ef0deab7 100644 --- a/test/doc/cyl_bessel_j0.cpp +++ b/test/doc/cyl_bessel_j0.cpp @@ -29,5 +29,13 @@ int main() eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); auto z0 = kyosu::complex(z1, z2); std::cout << z0 << "\n -> " << cyl_bessel_j0(z0) << "\n"; + + eve::wide> l(1.0, 15.0, 40.0, 60.0, eve::inf(eve::as()),eve::minf(eve::as()), eve::nan(eve::as()), 0.0); + auto lims = kyosu::complex(l); + std::cout << lims << " -> \n" << cyl_bessel_j0(lims) << "\n"; + std::cout << l << " -> \n" << cyl_bessel_j0(l) << "\n"; + auto nanan = kyosu::complex(eve::nan(eve::as()), eve::nan(eve::as())); + std::cout << nanan << " -> \n" << cyl_bessel_j0(nanan) << std::endl; + return 0; } diff --git a/test/doc/cyl_bessel_kn.cpp b/test/doc/cyl_bessel_kn.cpp new file mode 100644 index 00000000..f255f069 --- /dev/null +++ b/test/doc/cyl_bessel_kn.cpp @@ -0,0 +1,17 @@ +#include +#include +#include + +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)); + for(int n=2; n <= 10; ++n) + { + auto kn = kyosu::cyl_bessel_kn(n, z); + std::cout << "n " << n << ", z " << z << std::endl; + std::cout << "cyl_bessel_kn(n, z) " << kn << std::endl; + } + return 0; +} diff --git a/test/doc/cyl_bessel_y0.cpp b/test/doc/cyl_bessel_y0.cpp new file mode 100644 index 00000000..b4350f33 --- /dev/null +++ b/test/doc/cyl_bessel_y0.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::cyl_bessel_y0; + std::cout.precision(16); + std::cout << std::scientific; + std::cout << "1.0 " << " -> " << cyl_bessel_y0(1.0) << "\n"; + std::cout << "1+0i " << " -> " << cyl_bessel_y0(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << cyl_bessel_y0(15.0) << "\n"; + std::cout << "15+0i " << " -> " << cyl_bessel_y0(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << cyl_bessel_y0(40.0) << "\n"; + std::cout << "40+0i " << " -> " << cyl_bessel_y0(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << cyl_bessel_y0(60.0) << "\n"; + std::cout << "60+0i " << " -> " << cyl_bessel_y0(kyosu::complex(60.0, 0.0)) << "\n"; + + eve::wide> z(1.0, 15.0, 40.0, 60.0); + auto zz = kyosu::complex(z); + std::cout << z << " -> \n" << cyl_bessel_y0(z) << "\n"; + std::cout << zz << " -> \n" << cyl_bessel_y0(zz) << "\n"; + + std::cout << "1.0 " << " -> " << cyl_bessel_y0(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << cyl_bessel_y0(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_bessel_y0(kyosu::quaternion(1.0, 36.0, 2.0, 1.5)) << "\n"; + + + eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); + auto z0 = kyosu::complex(z1, z2); + std::cout << z0 << "\n -> " << cyl_bessel_y0(z0) << "\n"; + return 0; +} diff --git a/test/doc/cyl_bessel_y1.cpp b/test/doc/cyl_bessel_y1.cpp new file mode 100644 index 00000000..0ca817f2 --- /dev/null +++ b/test/doc/cyl_bessel_y1.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::cyl_bessel_y1; + std::cout.precision(16); + std::cout << std::scientific; + std::cout << "1.0 " << " -> " << cyl_bessel_y1(1.0) << "\n"; + std::cout << "1+0i " << " -> " << cyl_bessel_y1(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << cyl_bessel_y1(15.0) << "\n"; + std::cout << "15+0i " << " -> " << cyl_bessel_y1(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << cyl_bessel_y1(40.0) << "\n"; + std::cout << "40+0i " << " -> " << cyl_bessel_y1(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << cyl_bessel_y1(60.0) << "\n"; + std::cout << "60+0i " << " -> " << cyl_bessel_y1(kyosu::complex(60.0, 0.0)) << "\n"; + + eve::wide> z(1.0, 15.0, 40.0, 60.0); + auto zz = kyosu::complex(z); + std::cout << z << " -> \n" << cyl_bessel_y1(z) << "\n"; + std::cout << zz << " -> \n" << cyl_bessel_y1(zz) << "\n"; + + std::cout << "1.0 " << " -> " << cyl_bessel_y1(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << cyl_bessel_y1(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_bessel_y1(kyosu::quaternion(1.0, 36.0, 2.0, 1.5)) << "\n"; + + + eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); + auto z0 = kyosu::complex(z1, z2); + std::cout << z0 << "\n -> " << cyl_bessel_y1(z0) << "\n"; + return 0; +} diff --git a/test/doc/cyl_bessel_yn.cpp b/test/doc/cyl_bessel_yn.cpp new file mode 100644 index 00000000..bc8bf246 --- /dev/null +++ b/test/doc/cyl_bessel_yn.cpp @@ -0,0 +1,34 @@ +#include +#include +#include + +int main() +{ + size_t n = 12; + using kyosu::cyl_bessel_yn; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << cyl_bessel_yn(n,-1.0) << "\n"; + std::cout << "1+0i " << " -> " << cyl_bessel_yn(n,kyosu::complex(-1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << cyl_bessel_yn(n,15.0) << "\n"; + std::cout << "15+0i " << " -> " << cyl_bessel_yn(n,kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << cyl_bessel_yn(n,40.0) << "\n"; + std::cout << "40+0i " << " -> " << cyl_bessel_yn(n,kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << cyl_bessel_yn(n,60.0) << "\n"; + std::cout << "60+0i " << " -> " << cyl_bessel_yn(n,kyosu::complex(60.0, 0.0)) << "\n"; + + eve::wide> z(1.0, 15.0, 40.0, 60.0); + auto zz = kyosu::complex(z); + std::cout << z << "\n -> " << cyl_bessel_yn(n,z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(cyl_bessel_yn(n,zz)) << "\n"; + + std::cout << "1.0 " << " -> " << cyl_bessel_yn(n,1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << cyl_bessel_yn(n,kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_bessel_yn(n,kyosu::quaternion(1.0, 36.0, 2.0, 1.5)) << "\n"; + + + eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); + auto z0 = kyosu::complex(z1, z2); + std::cout << z0 << " \n-> " << cyl_bessel_yn(n,z0) << "\n"; + return 0; +} diff --git a/test/unit/function/cyl_bessel_i0.cpp b/test/unit/function/cyl_bessel_i0.cpp index 78ce839e..e9814629 100644 --- a/test/unit/function/cyl_bessel_i0.cpp +++ b/test/unit/function/cyl_bessel_i0.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i0 over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{1.3556565953650257e+00, 3.3014253671480853e+00, -2.5475601159592665e+00, 1.1681025625526764e+00, 2.0457564384951734e+00, -1.2532453407268573e+00, -4.7295944112409973e+00, -2.2840688135673504e+00}; std::array im{4.0914290638525896e+00, -5.6356545999116547e-01, -3.6064274257830387e+00, 1.6951171853616265e+00, -1.3703722677989172e+00, -1.2563658033049907e+00, -6.3364633643107848e-01, 3.7665729797577785e+00}; @@ -34,7 +33,6 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_i0(c), res, 1.0e-7) << i << " <- " << c << '\n'; } } - std::cout.precision(17); std::array re{1.1554018562666468e-01, 5.3044189227319072e-01, 3.5245310048347744e-01, 3.2527350287655354e-01, 2.4688936247411308e+00, -5.1146428618899231e-01, 1.6635701653597774e+00, -1.7452778188732947e+00}; @@ -56,7 +54,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i0 over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) diff --git a/test/unit/function/cyl_bessel_i1.cpp b/test/unit/function/cyl_bessel_i1.cpp index b5db1d5e..84f2bf87 100644 --- a/test/unit/function/cyl_bessel_i1.cpp +++ b/test/unit/function/cyl_bessel_i1.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i1 over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{-4.1836770228521267e+01, -4.1682207527601065e+01, -4.0619579092445747e+01, -4.4234864207453732e+01, -4.3110416350751848e+01, -4.3626638817587562e+01, -4.4580540281737981e+01, -4.1708122005195847e+01}; std::array im{-4.4724569983549614e+01, -4.6490400566267638e+01, -4.7317594550520738e+01, -4.1337120165516517e+01, -4.1276020002765804e+01, -4.6670524375443286e+01, -4.9467141247871425e+01, -4.1317766188944873e+01}; @@ -33,7 +32,6 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_i1(c), res, 1.0e-7) << i << " <- " << c << '\n'; } } - std::cout.precision(17); std::array re{-2.0045932610658785e+00, 1.5610767253921298e-02, 1.2016964006673581e+00, -1.4746256231752448e+00, 3.8509464414269778e-01, 4.3119169358173060e-01, 1.8117761375602393e-01, 2.2535543048838065e+00}; std::array im{-1.3262936608790188e+00, 7.3358653580584687e-01, -7.7829729505626277e-01, 1.8039445371160530e+00, -3.9048140529175879e-01, 2.3733805921203750e+00, 4.1558677759532259e-01, -2.2418145868867645e-01}; @@ -53,7 +51,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" } }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i1 over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) diff --git a/test/unit/function/cyl_bessel_in.cpp b/test/unit/function/cyl_bessel_in.cpp index a92ce65c..96aea3b3 100644 --- a/test/unit/function/cyl_bessel_in.cpp +++ b/test/unit/function/cyl_bessel_in.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_in over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -16,7 +16,6 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" { if constexpr(sizeof(T) == 8) { - std::cout.precision(17); std::array re{-4.2214088790474964e+00, -7.8055061036477191e-01, 3.4677730339081023e+00, -3.5357366395596843e+00, 2.6849636588863888e+00, 4.1415226198517496e+00, -4.3623200824031505e+00, 1.0341209529062589e-01}; std::array im{1.0956977560286330e+00, -7.5270747683236494e-01, -2.3663854475056922e+00, 2.1367696130628167e-01, 3.8676287626125005e+00, 1.2360189302485036e+00, -1.6272444761141791e+00, 2.3032771326096082e+00}; @@ -52,7 +51,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_in over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) @@ -60,17 +59,25 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" ) (T a0, T a1) { + using u_t = eve::underlying_type_t; auto c = kyosu::complex(a0, a1); auto cb= kyosu::conj(c); auto cm= -c; auto cr= kyosu::complex(a0); auto ci= kyosu::complex(T(0), a1); + auto zer = kyosu::complex(T(0), T(0)); + auto one = kyosu::complex(T(1), T(0)); + using kyosu::cyl_bessel_in; - auto inc = cyl_bessel_in(3, c); - TTS_IEEE_EQUAL(inc, -cyl_bessel_in(3, cm)); - TTS_IEEE_EQUAL(inc, kyosu::conj(cyl_bessel_in(3, cb))); + + for(int i=0; i < 10; ++i) + { + auto inc = cyl_bessel_in(i, c); + TTS_IEEE_EQUAL(inc, eve::sign_alternate(u_t(i))*cyl_bessel_in(i, cm)); + TTS_IEEE_EQUAL(inc, kyosu::conj(cyl_bessel_in(i, cb))); 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(cyl_bessel_in(3, z), z); + TTS_IEEE_EQUAL(cyl_bessel_in(i, zer), i ? zer : one); + TTS_IEEE_EQUAL(inc, cyl_bessel_in(-i, c)) << i << '\n'; + } }; diff --git a/test/unit/function/cyl_bessel_j0.cpp b/test/unit/function/cyl_bessel_j0.cpp index ce52eef9..581e33a0 100644 --- a/test/unit/function/cyl_bessel_j0.cpp +++ b/test/unit/function/cyl_bessel_j0.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j0 over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{ e_t(-9.8403264736233908e+01), e_t(8.6287324687322851e+01), e_t(2.8677183208357349e+01), e_t(7.5642931882481321e+01), e_t(9.6351793077594323e+01) , e_t(-7.8716695947540288e+01), e_t(-5.0019457963533974e+01), e_t(5.9030075032843811e+01), e_t(-7.5596124248257235e+00), e_t(8.5810875418338028e+01)}; @@ -37,7 +36,6 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_j0(im[i]), kyosu::real(kyosu::cyl_bessel_j0(kyosu::complex(im[i], e_t(0.0)))), 1.0e-4)<< re[i] << '\n'; } } - std::cout.precision(17); std::array re{2.175102683907268e+00, 3.007738475083643e-01, 3.571791197760752e+00, 1.916834892035037e+00, 2.309329011574615e+00, 1.116302757926106e+00, 3.454188251414261e+00, 2.504516854443353e+00}; std::array im{ 3.008362577802520e+00, 3.733804275914350e+00, 2.988312632706299e+00, 4.411110052057489e+00, 3.055249958745497e+00, 2.508994160634028e+00, 3.100976745146938e-01, 3.854880352108083e+00}; @@ -54,7 +52,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" } }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j0 over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) diff --git a/test/unit/function/cyl_bessel_j1.cpp b/test/unit/function/cyl_bessel_j1.cpp index d3730aec..2de8f892 100644 --- a/test/unit/function/cyl_bessel_j1.cpp +++ b/test/unit/function/cyl_bessel_j1.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j1 over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{ -9.8403264736233908e+01, 8.6287324687322851e+01, 2.8677183208357349e+01, 7.5642931882481321e+01, 9.6351793077594323e+01 , -7.8716695947540288e+01, -5.0019457963533974e+01, 5.9030075032843811e+01, -7.5596124248257235e+00, 8.5810875418338028e+01}; std::array im{-1.9271037994246587e+01, 7.9043790847200569e+01, -3.0599555605285602e+01, -7.1720588923445590e+01, -6.1110206702938363e+01 @@ -34,7 +33,6 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_j1(c), res, 1.0e-7) << i << " <- " << c << '\n'; } } - std::cout.precision(17); std::array re{2.175102683907268e+00, 3.007738475083643e-01, 3.571791197760752e+00, 1.916834892035037e+00, 2.309329011574615e+00, 1.116302757926106e+00, 3.454188251414261e+00, 2.504516854443353e+00}; std::array im{ 3.008362577802520e+00, 3.733804275914350e+00, 2.988312632706299e+00, 4.411110052057489e+00, 3.055249958745497e+00, 2.508994160634028e+00, 3.100976745146938e-01, 3.854880352108083e+00}; @@ -52,7 +50,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j1 over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) diff --git a/test/unit/function/cyl_bessel_jn.cpp b/test/unit/function/cyl_bessel_jn.cpp index 02cdd7e4..a8c7c70c 100644 --- a/test/unit/function/cyl_bessel_jn.cpp +++ b/test/unit/function/cyl_bessel_jn.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_jn over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -16,7 +16,6 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" { if constexpr(sizeof(T) == 8) { - std::cout.precision(17); std::array re{ -9.8403264736233908e+01, 8.6287324687322851e+01, 2.8677183208357349e+01, 7.5642931882481321e+01, 9.6351793077594323e+01 , -7.8716695947540288e+01, -5.0019457963533974e+01, 5.9030075032843811e+01, -7.5596124248257235e+00, 8.5810875418338028e+01}; std::array im{-1.9271037994246587e+01, 7.9043790847200569e+01, -3.0599555605285602e+01, -7.1720588923445590e+01, -6.1110206702938363e+01 @@ -53,7 +52,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_jn over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) @@ -61,17 +60,25 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" ) (T a0, T a1) { + using u_t = eve::underlying_type_t; auto c = kyosu::complex(a0, a1); auto cb= kyosu::conj(c); auto cm= -c; auto cr= kyosu::complex(a0); auto ci= kyosu::complex(T(0), a1); + auto zer = kyosu::complex(T(0), T(0)); + auto one = kyosu::complex(T(1), T(0)); + using kyosu::cyl_bessel_jn; - auto jnc = cyl_bessel_jn(3, c); - TTS_IEEE_EQUAL(jnc, -cyl_bessel_jn(3, cm)) << "c " << c << "\n"; - TTS_IEEE_EQUAL(jnc, kyosu::conj(cyl_bessel_jn(3, cb))); - 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(cyl_bessel_jn(3, z), z); + + for(int i=0; i < 10 ; ++i) + { + auto jnc = cyl_bessel_jn(i, c); + TTS_IEEE_EQUAL(jnc, eve::sign_alternate(u_t(i))*cyl_bessel_jn(i, cm)) << "c " << c << "\n"; + TTS_IEEE_EQUAL(jnc, kyosu::conj(cyl_bessel_jn(i, cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + TTS_IEEE_EQUAL(cyl_bessel_jn(i, zer), i ? zer : one); + TTS_IEEE_EQUAL(jnc, eve::sign_alternate(u_t(i))*cyl_bessel_jn(-i, c)) << i << '\n'; + } }; diff --git a/test/unit/function/cyl_bessel_k0.cpp b/test/unit/function/cyl_bessel_k0.cpp new file mode 100644 index 00000000..b4ac046b --- /dev/null +++ b/test/unit/function/cyl_bessel_k0.cpp @@ -0,0 +1,56 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k0 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{2.3584877349134130e+00, 9.0164385767932564e-01, -1.5393599740361523e+00, 1.5668095462351217e+00, 1.6179702389500128e+00, -5.5681137226619104e-01, -2.2990386896893691e+00, 1.9300383451287466e+00}; + + std::array im{-8.0071372624748560e-01, 1.9823595245805736e+00, -1.3978761922193099e+00, -1.4654991890321578e+00, -5.3416883666219606e-01, -1.9265301646062505e+00, -1.6132600067568537e+00, -2.3575615216282779e+00}; + + std::array reres{4.1782770293350736e-02, -2.7083150384809346e-01, -3.4769747409921696e+00, -4.0882933553504530e-02, 1.4002108258106824e-01, -1.4548439736185870e+00, -7.3832781780344741e+00, -9.4738063303879627e-02}; + + std::array imres{5.8709121187673685e-02, -1.9749705340577356e-01, 2.3664976923056558e+00, 1.6665625358420463e-01, 1.1276125697599239e-01, 5.3154381678401141e-01, 2.2763972412185578e+00, 3.6519862966414102e-02}; + 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::cyl_bessel_k0(c), res, 1.0e-4) << i << " <- " << c << '\n'; + } + } +}; + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k0 over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + auto c = kyosu::complex(a0, a1); + auto cb= kyosu::conj(c); + auto cr= kyosu::complex(a0); + using kyosu::cyl_bessel_k0; + auto k0c = cyl_bessel_k0(c); + auto inf = kyosu::complex(eve::inf(eve::as())); + auto z = kyosu::complex(T(0), T(0)); + + TTS_IEEE_EQUAL(k0c, kyosu::conj(cyl_bessel_k0(cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_IEEE_EQUAL(cyl_bessel_k0(z), inf); + TTS_IEEE_EQUAL(k0c, cyl_bessel_k0(c)); +}; diff --git a/test/unit/function/cyl_bessel_k1.cpp b/test/unit/function/cyl_bessel_k1.cpp new file mode 100644 index 00000000..64afaefe --- /dev/null +++ b/test/unit/function/cyl_bessel_k1.cpp @@ -0,0 +1,57 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k1 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{-7.3954299670486945e-01, -3.7928780215327196e-02, -1.4509312401953633e+00, 2.2964197111982774e+00, -7.9699549838500194e-01, 5.4099472045802843e-01, -5.8957905551605738e-01, -5.1146527014834298e-01}; + + std::array im{-2.1345822747239773e+00, -3.3484907524396035e-01, 2.2665940322875286e+00, -1.0045032564810230e+00, -2.4392625060316648e+00, 2.6547418070570761e-01, -2.3352892270792931e+00, 1.3879256025736524e+00}; + + std::array reres{-1.7027762686598518e+00, -5.8170412824877404e-01, -2.9141064941382457e+00, 2.7872334701787467e-02, -1.5779496640047279e+00, 1.1145462290717538e+00, -1.3982594001564612e+00, -1.4628928686130445e+00}; + + std::array imres{-1.9261773906655510e-01, 3.2589115393837513e+00, 1.0075935641554303e+00, 8.5373399021791460e-02, -6.4553389686343909e-01, -7.8086954815494913e-01, -3.6962699681190719e-01, -8.7620662002477645e-01}; + 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::cyl_bessel_k1(c), res, 1.1e-6) << i << " <- " << c << " -- "<< kyosu::arg(c) << '\n'; + } + } +}; + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k1 over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + auto c = kyosu::complex(a0, a1); + auto cb= kyosu::conj(c); + auto cr= kyosu::complex(a0); + using kyosu::cyl_bessel_k1; + auto k1c = cyl_bessel_k1(c); + auto inf = kyosu::complex(eve::inf(eve::as())); + + TTS_IEEE_EQUAL(k1c, kyosu::conj(cyl_bessel_k1(cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + auto z = kyosu::complex(T(0), T(0)); + TTS_IEEE_EQUAL(cyl_bessel_k1(z), inf); + TTS_IEEE_EQUAL(k1c, cyl_bessel_k1(c)); + +}; diff --git a/test/unit/function/cyl_bessel_kn.cpp b/test/unit/function/cyl_bessel_kn.cpp new file mode 100644 index 00000000..401e0417 --- /dev/null +++ b/test/unit/function/cyl_bessel_kn.cpp @@ -0,0 +1,64 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_in over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + if constexpr(sizeof(T) == 8) + { + std::array re{1.4376454063874640e+00, -1.0673204704071071e-01, -1.0828666774506268e+00, -7.9765980000519210e-01, 4.2409189866361996e-01,7.6608496996954323e-01, -2.2849830522845975e+00, 4.3813679108628711e-01}; + + std::array im{-2.5175040147125849e-01, 6.2364116715495577e-01, 1.1653810042557073e+00, 1.3466103438716175e+00, -1.8438471344336298e+00, -2.0952595273684760e+00, -7.7012218588605907e-01, 1.3705948670934771e+00}; + + std::array reres{1.6867419926930389e+00, 1.5635267649025545e+01, 1.7030997128713221e+00, 2.2880149616855170e+00, -6.6737532774450159e-01, -5.4294479427715847e-01, -1.2185670369281292e+00, -2.2366901615818047e+00}; + + std::array imres{1.1484897053630303e+00, 2.9235356194719746e+01, -4.8470591465484553e-01, 8.1442309859938877e-01, -1.6108318991995840e+00, -9.0484737739480448e-01, 7.8990385458658041e-01, 2.4137425589017183e+00}; + + for(size_t i=0; i < re.size(); ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_kn(3, c), res, 1.0e-6) << i << " <- " << c << '\n'; + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_kn 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; + auto c = kyosu::complex(a0, a1); + auto cb= kyosu::conj(c); + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + auto zer = kyosu::complex(T(0), T(0)); + auto inf = kyosu::complex(eve::inf(eve::as())); + + using kyosu::cyl_bessel_kn; + + for(int i=0; i < 10; ++i) + { + auto knc = cyl_bessel_kn(i, c); + TTS_IEEE_EQUAL(knc, kyosu::conj(cyl_bessel_kn(i, cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + TTS_IEEE_EQUAL(cyl_bessel_kn(i, zer), inf) << i << '\n'; + TTS_IEEE_EQUAL(knc, cyl_bessel_kn(-i, c)); + } +}; diff --git a/test/unit/function/cyl_bessel_y0.cpp b/test/unit/function/cyl_bessel_y0.cpp new file mode 100644 index 00000000..7b8f373f --- /dev/null +++ b/test/unit/function/cyl_bessel_y0.cpp @@ -0,0 +1,59 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_y0 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{2.4955917718518910e+00, -1.6982259765731111e+00, 1.2315182267355251e+00, 2.1258824905764118e+00, -2.1711533910447676e+00, 9.1742698952666013e-01, -2.5479731910545422e-02, 1.4924891723416751e+00}; + + std::array im{-5.0148281002365624e-01, 2.0253551112861641e+00, -5.0175898965530097e-01, -4.0455814438323523e-01, -2.1923016231568171e-01, -1.3185119199506317e+00, 1.0779142409536568e+00, 1.2528910061509635e-01}; + + std::array reres{5.5469839057084858e-01, -1.8165994943959307e+00, 3.3680718591187470e-01, 5.6362426783532593e-01, 2.8645154466364819e-01, 6.0960351781866107e-01, -2.5569311956619334e-01, 3.8450178879111341e-01}; + + std::array imres{7.8480890054575703e-02, 5.8352584787301165e-01, -2.9298142302554275e-01, -1.2803370889276175e-02, -2.4399346361135993e-01, -9.9087866505382638e-01, 1.3205924969979352e+00, 5.2208242352489359e-02}; + + + 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::cyl_bessel_yn(0, c), res, 1.0e-4) << i << " <- " << c << '\n'; + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_y0 over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + 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::cyl_bessel_y0; + auto y0c = cyl_bessel_y0(c); + TTS_IEEE_EQUAL(y0c, kyosu::conj(cyl_bessel_y0(cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + auto z = kyosu::complex(T(0), T(0)); + auto minf= kyosu::complex(eve::minf(eve::as())); + TTS_IEEE_EQUAL(cyl_bessel_y0(z), minf); +}; diff --git a/test/unit/function/cyl_bessel_y1.cpp b/test/unit/function/cyl_bessel_y1.cpp new file mode 100644 index 00000000..04593e67 --- /dev/null +++ b/test/unit/function/cyl_bessel_y1.cpp @@ -0,0 +1,58 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_y1 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{6.5758573026995004e-01, -2.4649387808542378e+00, 1.9205766369481703e+00, -6.9739731404600436e-01, 4.4532863650021348e-01, 2.4757649834958393e+00, -5.3435309039155832e-01, -2.4349028761454212e+00}; + + std::array im{1.2035170114779326e+00, -3.5932055457567835e-01, 8.3662582734559288e-01, 2.2469103234867700e+00, -9.1023149008937376e-01, 2.4818838321237546e-01, 1.5847819595769823e+00, -1.3967310652706204e+00}; + + std::array reres{-8.1045682169396482e-01, 2.3820060953407891e-02, -9.3267933254955498e-02, -1.5790753753124993e+00, -7.5013324227681488e-01, 1.4423021239534972e-01, -8.5123355830924607e-01, 4.3341743806150423e-01}; + + std::array imres{6.3354487728865994e-01, 8.9036262959787016e-01, 5.0186176942328398e-01, -1.1543670775735402e+00, -5.8526040753707131e-01, 1.1142422613715078e-01, -4.3418909281427642e-01, 1.0069485325105600e+00}; + + 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::cyl_bessel_y1(c), res, 1.0e-4) << i << " <- " << c << '\n'; + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_y1 over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + 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::cyl_bessel_y1; + auto y1c = cyl_bessel_y1(c); + TTS_IEEE_EQUAL(y1c, kyosu::conj(cyl_bessel_y1(cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + auto z = kyosu::complex(T(0), T(0)); + auto minf= kyosu::complex(eve::minf(eve::as())); + TTS_IEEE_EQUAL(cyl_bessel_y1(z), minf); +}; diff --git a/test/unit/function/cyl_bessel_yn.cpp b/test/unit/function/cyl_bessel_yn.cpp new file mode 100644 index 00000000..21e43ba2 --- /dev/null +++ b/test/unit/function/cyl_bessel_yn.cpp @@ -0,0 +1,64 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_yn over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + if constexpr(sizeof(T) == 8) + { + std::array re{2.0309144503737393e+00, -1.5992515729134549e+00, 1.6982551054108397e+00, -1.3327549270623136e+00, 2.6585368919111341e-01, 7.3163244235853131e-01, -1.5491151507087859e+00, 2.4846599242220186e+00}; + + std::array im{-1.4137922349350807e+00, -1.7357125614937789e+00, 4.0762417134967421e-02, -1.4645033242641059e+00, 1.5914290999536518e+00, 2.3057243858262124e+00, -2.2075167685839387e+00, -5.7866716195750623e-01}; + + std::array reres{-4.4062036092080337e-01, -3.3967954461649441e-01, -1.5667619074526007e+00, -4.0477676250027339e-01, 5.9271179286244735e-01, 3.7250752149022126e-01, -4.0842410865380668e-01, -6.8424644495325704e-01}; + + std::array imres{-3.6401836760736062e-01, -4.6948726732250373e-01, 8.0811726440696674e-02, -6.2890286218778024e-01, -8.1624047120865173e-01, -4.2529834100001956e-01,-5.2072912299547713e-01, -2.6673624676594038e-01}; + + 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::cyl_bessel_yn(3, c), res, 1.0e-6) << i << " <- " << c << '\n'; + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_yn 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; + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(eve::abs(a0)); + auto zer = kyosu::complex(T(0), T(0)); + auto minf= kyosu::complex(eve::minf(eve::as())); + { + auto conjz = kyosu::conj(z); + + for(int i=0; i <10 ; ++i) + { + auto ynz = kyosu::cyl_bessel_yn(i, z); + TTS_IEEE_EQUAL(ynz, kyosu::conj(kyosu::cyl_bessel_yn(i, conjz))); + auto yre = kyosu::cyl_bessel_yn(i, kyosu::abs(re)); + TTS_EXPECT(eve::all(kyosu::is_real(yre)))<< i << " -> " << yre << '\n'; + auto yzer = kyosu::cyl_bessel_yn(i, zer); + TTS_IEEE_EQUAL(yzer, minf); + TTS_IEEE_EQUAL(ynz, eve::sign_alternate(u_t(i))*kyosu::cyl_bessel_yn(-i, z)) << i << '\n'; + } + } +};