From 21f9f3adc34e0c1339cdfa501afc705bf74dc16e Mon Sep 17 00:00:00 2001 From: jtlap Date: Mon, 13 Nov 2023 20:31:30 +0100 Subject: [PATCH] Spherical Bessel functions --- doc/index.md | 15 +- include/kyosu/details/cayleyify.hpp | 1 + .../besseli0.hpp => details/with_alloca.hpp} | 13 +- include/kyosu/functions.hpp | 30 ++ include/kyosu/functions/cyl_bessel_i1.hpp | 2 +- include/kyosu/functions/from_angle_axis.hpp | 4 - include/kyosu/functions/sinhc.hpp | 74 +++ include/kyosu/functions/sph_bessel_h1_0.hpp | 77 +++ include/kyosu/functions/sph_bessel_h1_1.hpp | 77 +++ include/kyosu/functions/sph_bessel_h1n.hpp | 77 +++ include/kyosu/functions/sph_bessel_h2_0.hpp | 77 +++ include/kyosu/functions/sph_bessel_h2_1.hpp | 77 +++ include/kyosu/functions/sph_bessel_h2n.hpp | 77 +++ include/kyosu/functions/sph_bessel_i1_0.hpp | 76 +++ include/kyosu/functions/sph_bessel_i1_1.hpp | 78 +++ include/kyosu/functions/sph_bessel_i1n.hpp | 78 +++ include/kyosu/functions/sph_bessel_i2_0.hpp | 75 +++ include/kyosu/functions/sph_bessel_i2_1.hpp | 78 +++ include/kyosu/functions/sph_bessel_i2n.hpp | 78 +++ include/kyosu/functions/sph_bessel_j0.hpp | 73 +++ include/kyosu/functions/sph_bessel_j1.hpp | 74 +++ include/kyosu/functions/sph_bessel_jn.hpp | 74 +++ include/kyosu/functions/sph_bessel_k0.hpp | 76 +++ include/kyosu/functions/sph_bessel_k1.hpp | 77 +++ include/kyosu/functions/sph_bessel_kn.hpp | 79 +++ include/kyosu/functions/sph_bessel_y0.hpp | 73 +++ include/kyosu/functions/sph_bessel_y1.hpp | 74 +++ include/kyosu/functions/sph_bessel_yn.hpp | 74 +++ include/kyosu/types/impl/bessel.hpp | 65 ++- include/kyosu/types/impl/detail/bessel_h.hpp | 151 ++++++ include/kyosu/types/impl/detail/bessel_i.hpp | 182 +++++++ include/kyosu/types/impl/detail/bessel_j.hpp | 465 ++++++++++++++++++ .../detail/{besselkn.hpp => bessel_k.hpp} | 47 +- include/kyosu/types/impl/detail/bessel_y.hpp | 169 +++++++ include/kyosu/types/impl/detail/besselhn.hpp | 46 -- include/kyosu/types/impl/detail/besseli1.hpp | 30 -- include/kyosu/types/impl/detail/besselin.hpp | 37 -- include/kyosu/types/impl/detail/besselj0.hpp | 5 +- include/kyosu/types/impl/detail/besselj1.hpp | 5 +- include/kyosu/types/impl/detail/besseljn.hpp | 2 - 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 | 2 +- include/kyosu/types/impl/detail/cylseq.hpp | 51 ++ include/kyosu/types/impl/trigo.hpp | 17 + test/doc/cyl_bessel_yn.cpp | 6 + test/doc/sinhc.cpp | 25 + test/doc/sph_bessel_h1_0.cpp | 33 ++ test/doc/sph_bessel_h1_1.cpp | 33 ++ test/doc/sph_bessel_h1n.cpp | 41 ++ test/doc/sph_bessel_h2_0.cpp | 33 ++ test/doc/sph_bessel_h2_1.cpp | 33 ++ test/doc/sph_bessel_h2n.cpp | 40 ++ test/doc/sph_bessel_i1_0.cpp | 33 ++ test/doc/sph_bessel_i1_1.cpp | 33 ++ test/doc/sph_bessel_i1n.cpp | 40 ++ test/doc/sph_bessel_i2_0.cpp | 33 ++ test/doc/sph_bessel_i2_1.cpp | 33 ++ test/doc/sph_bessel_i2n.cpp | 40 ++ test/doc/sph_bessel_j0.cpp | 33 ++ test/doc/sph_bessel_j1.cpp | 33 ++ test/doc/sph_bessel_jn.cpp | 41 ++ test/doc/sph_bessel_k0.cpp | 33 ++ test/doc/sph_bessel_k1.cpp | 34 ++ test/doc/sph_bessel_kn.cpp | 41 ++ test/doc/sph_bessel_y0.cpp | 33 ++ test/doc/sph_bessel_y1.cpp | 33 ++ test/doc/sph_bessel_yn.cpp | 40 ++ test/unit/function/sinhc.cpp | 36 ++ test/unit/function/sph_bessel_h0.cpp | 46 ++ test/unit/function/sph_bessel_h1.cpp | 46 ++ test/unit/function/sph_bessel_hn.cpp | 77 +++ test/unit/function/sph_bessel_i0.cpp | 62 +++ test/unit/function/sph_bessel_i1.cpp | 70 +++ test/unit/function/sph_bessel_in.cpp | 76 +++ test/unit/function/sph_bessel_j0.cpp | 50 ++ test/unit/function/sph_bessel_j1.cpp | 50 ++ test/unit/function/sph_bessel_jn.cpp | 67 +++ test/unit/function/sph_bessel_k0.cpp | 48 ++ test/unit/function/sph_bessel_k1.cpp | 51 ++ test/unit/function/sph_bessel_kn.cpp | 59 +++ test/unit/function/sph_bessel_y0.cpp | 50 ++ test/unit/function/sph_bessel_y1.cpp | 50 ++ test/unit/function/sph_bessel_yn.cpp | 64 +++ 85 files changed, 4539 insertions(+), 258 deletions(-) rename include/kyosu/{types/impl/detail/besseli0.hpp => details/with_alloca.hpp} (50%) create mode 100644 include/kyosu/functions/sinhc.hpp create mode 100644 include/kyosu/functions/sph_bessel_h1_0.hpp create mode 100644 include/kyosu/functions/sph_bessel_h1_1.hpp create mode 100644 include/kyosu/functions/sph_bessel_h1n.hpp create mode 100644 include/kyosu/functions/sph_bessel_h2_0.hpp create mode 100644 include/kyosu/functions/sph_bessel_h2_1.hpp create mode 100644 include/kyosu/functions/sph_bessel_h2n.hpp create mode 100644 include/kyosu/functions/sph_bessel_i1_0.hpp create mode 100644 include/kyosu/functions/sph_bessel_i1_1.hpp create mode 100644 include/kyosu/functions/sph_bessel_i1n.hpp create mode 100644 include/kyosu/functions/sph_bessel_i2_0.hpp create mode 100644 include/kyosu/functions/sph_bessel_i2_1.hpp create mode 100644 include/kyosu/functions/sph_bessel_i2n.hpp create mode 100644 include/kyosu/functions/sph_bessel_j0.hpp create mode 100644 include/kyosu/functions/sph_bessel_j1.hpp create mode 100644 include/kyosu/functions/sph_bessel_jn.hpp create mode 100644 include/kyosu/functions/sph_bessel_k0.hpp create mode 100644 include/kyosu/functions/sph_bessel_k1.hpp create mode 100644 include/kyosu/functions/sph_bessel_kn.hpp create mode 100644 include/kyosu/functions/sph_bessel_y0.hpp create mode 100644 include/kyosu/functions/sph_bessel_y1.hpp create mode 100644 include/kyosu/functions/sph_bessel_yn.hpp create mode 100644 include/kyosu/types/impl/detail/bessel_h.hpp create mode 100644 include/kyosu/types/impl/detail/bessel_i.hpp create mode 100644 include/kyosu/types/impl/detail/bessel_j.hpp rename include/kyosu/types/impl/detail/{besselkn.hpp => bessel_k.hpp} (53%) create mode 100644 include/kyosu/types/impl/detail/bessel_y.hpp delete mode 100644 include/kyosu/types/impl/detail/besselhn.hpp delete mode 100644 include/kyosu/types/impl/detail/besseli1.hpp delete mode 100644 include/kyosu/types/impl/detail/besselin.hpp delete mode 100644 include/kyosu/types/impl/detail/bessely0.hpp delete mode 100644 include/kyosu/types/impl/detail/bessely1.hpp delete mode 100644 include/kyosu/types/impl/detail/besselyn.hpp create mode 100644 include/kyosu/types/impl/detail/cylseq.hpp create mode 100644 test/doc/sinhc.cpp create mode 100644 test/doc/sph_bessel_h1_0.cpp create mode 100644 test/doc/sph_bessel_h1_1.cpp create mode 100644 test/doc/sph_bessel_h1n.cpp create mode 100644 test/doc/sph_bessel_h2_0.cpp create mode 100644 test/doc/sph_bessel_h2_1.cpp create mode 100644 test/doc/sph_bessel_h2n.cpp create mode 100644 test/doc/sph_bessel_i1_0.cpp create mode 100644 test/doc/sph_bessel_i1_1.cpp create mode 100644 test/doc/sph_bessel_i1n.cpp create mode 100644 test/doc/sph_bessel_i2_0.cpp create mode 100644 test/doc/sph_bessel_i2_1.cpp create mode 100644 test/doc/sph_bessel_i2n.cpp create mode 100644 test/doc/sph_bessel_j0.cpp create mode 100644 test/doc/sph_bessel_j1.cpp create mode 100644 test/doc/sph_bessel_jn.cpp create mode 100644 test/doc/sph_bessel_k0.cpp create mode 100644 test/doc/sph_bessel_k1.cpp create mode 100644 test/doc/sph_bessel_kn.cpp create mode 100644 test/doc/sph_bessel_y0.cpp create mode 100644 test/doc/sph_bessel_y1.cpp create mode 100644 test/doc/sph_bessel_yn.cpp create mode 100644 test/unit/function/sinhc.cpp create mode 100644 test/unit/function/sph_bessel_h0.cpp create mode 100644 test/unit/function/sph_bessel_h1.cpp create mode 100644 test/unit/function/sph_bessel_hn.cpp create mode 100644 test/unit/function/sph_bessel_i0.cpp create mode 100644 test/unit/function/sph_bessel_i1.cpp create mode 100644 test/unit/function/sph_bessel_in.cpp create mode 100644 test/unit/function/sph_bessel_j0.cpp create mode 100644 test/unit/function/sph_bessel_j1.cpp create mode 100644 test/unit/function/sph_bessel_jn.cpp create mode 100644 test/unit/function/sph_bessel_k0.cpp create mode 100644 test/unit/function/sph_bessel_k1.cpp create mode 100644 test/unit/function/sph_bessel_kn.cpp create mode 100644 test/unit/function/sph_bessel_y0.cpp create mode 100644 test/unit/function/sph_bessel_y1.cpp create mode 100644 test/unit/function/sph_bessel_yn.cpp diff --git a/doc/index.md b/doc/index.md index 87031da1..66a05526 100644 --- a/doc/index.md +++ b/doc/index.md @@ -106,7 +106,7 @@ Operators --------- Operators (as said before) `+`, `-`, `*` and `/` can be used in infix form and can mix cayley-dickson values of - different dimensinalities. Of course the biggest dimensionlity is recovered in the output. + different dimensionalities. Of course the biggest dimensionlity is recovered in the output. Prefix forms are also provided as `add`, `sub`, `multiply` and `div`. Also plus and minus for unary versions. @@ -115,10 +115,9 @@ The left division sometimes necessary if the dimensionality is greater than 2 is Functions --------- -Most **KYOSU** callables are usable with all cayley_dickson types. The exception being mainly special -complex functions and rotation related quaternion usage. +Most **KYOSU** callables are usable with all cayley_dickson types. The exceptions mainly being rotation related quaternion usage. -@warning: **EVE** callables that correspond to mathematical functions that +@warning **EVE** callables that correspond to mathematical functions that are only defined on a proper part of the real axis as, for example, `acos` DOES NOT ever provide the same result if called in **EVE** or **KYOSU** context. @@ -163,8 +162,14 @@ complex functions and rotation related quaternion usage. | [sqr](@ref kyosu::sqr ) | [sqr_abs](@ref kyosu::sqr_abs ) | [sqrt](@ref kyosu::sqrt ) | [tan](@ref kyosu::tan ) | [tanpi](@ref kyosu::tanpi ) | | [tanh](@ref kyosu::tanh ) | [to_polar](@ref kyosu::to_polar ) | [tgamma](@ref kyosu::tgamma ) | [trunc](@ref kyosu::trunc ) | [zeta](@ref kyosu::zeta ) | + * Bessel functions - * callables usable with quaternion complex and real only + They deserve a separate list to avoid verbosity: + + * cylindrical Bessel functions: cyl_bessel_xxx are all defined for xxx being j0, j1, jn, y0, y1, yn, i0, i1, in, h1_0, h1_1, h1n, h2_0, h2_1, h2n, k1, k0, kn. + * spherical Bessel functions: sph_bessel_xxx are all defined for xxx being j0, j1, jn, y0, y1, yn, i1_0, i1_1, i1n, i2_0, i2_1, i2n, h1_0, h1_1, h1n, h2_0, h2_1 h2n, k1, k0, kn. + + * Callables usable with quaternion complex and real only These functions are related to \f$\mathbb{R}^3\f$ rotations diff --git a/include/kyosu/details/cayleyify.hpp b/include/kyosu/details/cayleyify.hpp index 8ae31006..8d6a38e3 100644 --- a/include/kyosu/details/cayleyify.hpp +++ b/include/kyosu/details/cayleyify.hpp @@ -23,6 +23,7 @@ auto cayley_extend_rev(auto f, auto z1, auto z2) return kyosu::real(c) + kyosu::ipart(c)*kyosu::sign(p); } + auto cayley_extend2(auto f, auto z, auto... rs) { auto p = kyosu::pure(z); diff --git a/include/kyosu/types/impl/detail/besseli0.hpp b/include/kyosu/details/with_alloca.hpp similarity index 50% rename from include/kyosu/types/impl/detail/besseli0.hpp rename to include/kyosu/details/with_alloca.hpp index f821855c..1c320591 100644 --- a/include/kyosu/types/impl/detail/besseli0.hpp +++ b/include/kyosu/details/with_alloca.hpp @@ -6,17 +6,16 @@ */ //====================================================================================================================== #pragma once -#include +#include namespace kyosu::_ { - //===------------------------------------------------------------------------------------------- - // cyl_bessel_i0 - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, Z z) noexcept + template F> + decltype(auto) with_alloca(auto size, F f) { - return cyl_bessel_j0(complex(-ipart(z), real(z))); + T* p = (T*)(__builtin_alloca_with_align(size*sizeof(T), 8*alignof(T))); + return f(p); } + } diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index 92af6a8f..e0b5f199 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -163,8 +163,38 @@ #include #include #include +#include #include #include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + #include #include #include diff --git a/include/kyosu/functions/cyl_bessel_i1.hpp b/include/kyosu/functions/cyl_bessel_i1.hpp index ecbfe37c..ba089334 100644 --- a/include/kyosu/functions/cyl_bessel_i1.hpp +++ b/include/kyosu/functions/cyl_bessel_i1.hpp @@ -40,7 +40,7 @@ namespace kyosu //! @{ //! @var cyl_bessel_i1 //! @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. +//! \f$ I_1(x)= iJ_1(ix) \f$ extended to the complex plane and cayley_dickson algebras. //! //! It is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 0\f$. //! diff --git a/include/kyosu/functions/from_angle_axis.hpp b/include/kyosu/functions/from_angle_axis.hpp index 4ea32a40..64658ef0 100644 --- a/include/kyosu/functions/from_angle_axis.hpp +++ b/include/kyosu/functions/from_angle_axis.hpp @@ -26,10 +26,6 @@ namespace kyosu::tags { auto fn = callable_from_angle_axis{}; return fn(angle, axis, normalize); -// using e_t = decltype(axis[0]+angle); -// auto q = quaternion(e_t(0), e_t(axis[0]), e_t(axis[1]), e_t(axis[2])); -// auto [s, c] = eve::sincos(angle*eve::half(eve::as(angle))); -// return c+s*q; } template diff --git a/include/kyosu/functions/sinhc.hpp b/include/kyosu/functions/sinhc.hpp new file mode 100644 index 00000000..4fbc24b7 --- /dev/null +++ b/include/kyosu/functions/sinhc.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_sinhc: eve::elementwise + { + using callable_tag_type = callable_sinhc; + + KYOSU_DEFERS_CALLABLE(sinhc_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::sinhc(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 sinhc +//! @brief Computes the hyperbolic sine cardinal of the argument. +//! +//! **Defined in Header** +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sinhc(T z) noexcept; +//! template constexpr auto sinhc(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! Returns the hyperbolic sine cardinal of the argument. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sinhc.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sinhc sinhc = {}; +} diff --git a/include/kyosu/functions/sph_bessel_h1_0.hpp b/include/kyosu/functions/sph_bessel_h1_0.hpp new file mode 100644 index 00000000..19b4bd19 --- /dev/null +++ b/include/kyosu/functions/sph_bessel_h1_0.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_sph_bessel_h1_0: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_h1_0; + + KYOSU_DEFERS_CALLABLE(sph_bessel_h1_0_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept + { + return complex(eve::sph_bessel_j0(z), eve::sph_bessel_y0(z)); + } + + template + KYOSU_FORCEINLINE auto operator()(T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target1)) + { + return eve::tag_invoke(*this, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var sph_bessel_h1_0 +//! @brief Computes the spherical Bessel/hankel functions of the third kind, +//! \f$ h_0^{(1)}(z) = j_1(z)+iy_1(z)\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_h1_0(int n, T z) noexcept; +//! template constexpr auto sph_bessel_h1_0(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$h_0^{(1)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_h1_0.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_h1_0 sph_bessel_h1_0 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_h1_1.hpp b/include/kyosu/functions/sph_bessel_h1_1.hpp new file mode 100644 index 00000000..5d30ee2d --- /dev/null +++ b/include/kyosu/functions/sph_bessel_h1_1.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_sph_bessel_h1_1: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_h1_1; + + KYOSU_DEFERS_CALLABLE(sph_bessel_h1_1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept + { + return complex(eve::sph_bessel_j1(z), eve::sph_bessel_y1(z)); + } + + template + KYOSU_FORCEINLINE auto operator()(T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target1)) + { + return eve::tag_invoke(*this, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var sph_bessel_h1_1 +//! @brief Computes the spherical Bessel/hankel functions of the third kind, +//! \f$ h_1^{(1)}(z) = j_1(z)+iy_1(z)\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_h1_1(int n, T z) noexcept; +//! template constexpr auto sph_bessel_h1_1(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$h_1^{(1)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_h1_1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_h1_1 sph_bessel_h1_1 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_h1n.hpp b/include/kyosu/functions/sph_bessel_h1n.hpp new file mode 100644 index 00000000..e4c4aeca --- /dev/null +++ b/include/kyosu/functions/sph_bessel_h1n.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_sph_bessel_h1n: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_h1n; + + KYOSU_DEFERS_CALLABLE(sph_bessel_h1n_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& z) noexcept + { + return complex(eve::sph_bessel_jn(n, z), eve::sph_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 sph_bessel_h1n +//! @brief Computes the spherical 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 sph_bessel_h1n(int n, T z) noexcept; +//! template constexpr auto sph_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/sph_bessel_h1n.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_h1n sph_bessel_h1n = {}; +} diff --git a/include/kyosu/functions/sph_bessel_h2_0.hpp b/include/kyosu/functions/sph_bessel_h2_0.hpp new file mode 100644 index 00000000..b8b90f37 --- /dev/null +++ b/include/kyosu/functions/sph_bessel_h2_0.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_sph_bessel_h2_0: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_h2_0; + + KYOSU_DEFERS_CALLABLE(sph_bessel_h2_0_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept + { + return complex(eve::sph_bessel_j1(z), -eve::sph_bessel_y1(z)); + } + + template + KYOSU_FORCEINLINE auto operator()(T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target1)) + { + return eve::tag_invoke(*this, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var sph_bessel_h2_0 +//! @brief Computes the spherical Bessel/hankel functions of the third kind, +//! \f$ h_0^{(2)}(z) = j_0(z)-iy_0(z)\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_h2_1(int n, T z) noexcept; +//! template constexpr auto sph_bessel_h2_1(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$h_1^{(2)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_h2_1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_h2_0 sph_bessel_h2_0 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_h2_1.hpp b/include/kyosu/functions/sph_bessel_h2_1.hpp new file mode 100644 index 00000000..1a8ae30c --- /dev/null +++ b/include/kyosu/functions/sph_bessel_h2_1.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_sph_bessel_h2_1: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_h2_1; + + KYOSU_DEFERS_CALLABLE(sph_bessel_h2_1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept + { + return complex(eve::sph_bessel_j1(z), -eve::sph_bessel_y1(z)); + } + + template + KYOSU_FORCEINLINE auto operator()(T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target1)) + { + return eve::tag_invoke(*this, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var sph_bessel_h2_1 +//! @brief Computes the spherical Bessel/hankel functions of the third kind, +//! \f$ h_1^{(2)}(z) = j_1(z)-iy_1(z)\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_h2_1(int n, T z) noexcept; +//! template constexpr auto sph_bessel_h2_1(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$h_1^{(2)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_h2_1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_h2_1 sph_bessel_h2_1 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_h2n.hpp b/include/kyosu/functions/sph_bessel_h2n.hpp new file mode 100644 index 00000000..19e480dd --- /dev/null +++ b/include/kyosu/functions/sph_bessel_h2n.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_sph_bessel_h2n: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_h2n; + + KYOSU_DEFERS_CALLABLE(sph_bessel_h2n_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& z) noexcept + { + return complex(eve::sph_bessel_jn(n, z), -eve::sph_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 sph_bessel_h2n +//! @brief Computes the spherical Bessel/hankel functions of the third kind, +//! \f$ h_n^{(2)}(z) = j_n(z)-iy_n(z)\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_h2n(int n, T z) noexcept; +//! template constexpr auto sph_bessel_h2n(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$h_n^{(2)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_h2n.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_h2n sph_bessel_h2n = {}; +} diff --git a/include/kyosu/functions/sph_bessel_i1_0.hpp b/include/kyosu/functions/sph_bessel_i1_0.hpp new file mode 100644 index 00000000..25bb32ed --- /dev/null +++ b/include/kyosu/functions/sph_bessel_i1_0.hpp @@ -0,0 +1,76 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_sph_bessel_i1_0: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_i1_0; + + KYOSU_DEFERS_CALLABLE(sph_bessel_i1_0_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept + { + return eve::sinhc(z)/z; + } + + 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 sph_bessel_i1_0 +//! @brief Computes the Bessel function, +//! \f$ i_0^{(1)}(z) = j_n(iz)\f$ extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_i1_0(T z) noexcept; +//! template constexpr T sph_bessel_i1_0(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$i_0^{(1)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_i1_0.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_i1_0 sph_bessel_i1_0 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_i1_1.hpp b/include/kyosu/functions/sph_bessel_i1_1.hpp new file mode 100644 index 00000000..964cff78 --- /dev/null +++ b/include/kyosu/functions/sph_bessel_i1_1.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_sph_bessel_i1_1: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_i1_1; + + KYOSU_DEFERS_CALLABLE(sph_bessel_i1_1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept + { + auto [sh, ch] = eve::sinhcosh(z); + auto rz = rec(z); + return eve::fnma(sh, rz, ch)*rz; + } + + 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 sph_bessel_i1_1 +//! @brief Computes the Bessel function, +//! \f$ i_1^{(1)}(z) = -i j_1(iz)\f$ extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_i1_1(T z) noexcept; +//! template constexpr T sph_bessel_i1_1(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$i_1^{(1)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_i1_1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_i1_1 sph_bessel_i1_1 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_i1n.hpp b/include/kyosu/functions/sph_bessel_i1n.hpp new file mode 100644 index 00000000..023fcaaa --- /dev/null +++ b/include/kyosu/functions/sph_bessel_i1n.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_sph_bessel_i1n: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_i1n; + + KYOSU_DEFERS_CALLABLE(sph_bessel_i1n_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& z) noexcept + { + auto fn = callable_sph_bessel_i1n{}; + return fn(n, complex(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 sph_bessel_i1n +//! @brief Computes the spherical Bessel functions +//! \f$ i_n^{(1)}(z) = i^{-n}j_n(iz)\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_i1n(int n, T z) noexcept; +//! template constexpr auto sph_bessel_i1n(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$i_n^{(1)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_i1n.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_i1n sph_bessel_i1n = {}; +} diff --git a/include/kyosu/functions/sph_bessel_i2_0.hpp b/include/kyosu/functions/sph_bessel_i2_0.hpp new file mode 100644 index 00000000..40894190 --- /dev/null +++ b/include/kyosu/functions/sph_bessel_i2_0.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_sph_bessel_i2_0: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_i2_0; + + KYOSU_DEFERS_CALLABLE(sph_bessel_i2_0_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept + { return eve::cosh(z)/z; + } + + 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 sph_bessel_i2_0 +//! @brief Computes the Bessel function, +//! \f$ i_0^{(2)}(z) = -i y_n(iz)\f$ extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_i2_0(T z) noexcept; +//! template constexpr T sph_bessel_i2_0(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$i_0^{(2)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_i2_0.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_i2_0 sph_bessel_i2_0 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_i2_1.hpp b/include/kyosu/functions/sph_bessel_i2_1.hpp new file mode 100644 index 00000000..bc6770d8 --- /dev/null +++ b/include/kyosu/functions/sph_bessel_i2_1.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_sph_bessel_i2_1: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_i2_1; + + KYOSU_DEFERS_CALLABLE(sph_bessel_i2_1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept + { + auto [sh, ch] = eve::sinhcosh(z); + auto rz = eve::rec(z); + return eve::fnma(ch, rz, sh)*rz; + } + + 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 sph_bessel_i2_1 +//! @brief Computes the Bessel function, +//! \f$ i_1^{(2)}(z) = -y_1(iz)\f$ extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_i2_1(T z) noexcept; +//! template constexpr T sph_bessel_i2_1(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$i_1^{(2)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_i2_1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_i2_1 sph_bessel_i2_1 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_i2n.hpp b/include/kyosu/functions/sph_bessel_i2n.hpp new file mode 100644 index 00000000..7fdf560d --- /dev/null +++ b/include/kyosu/functions/sph_bessel_i2n.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_sph_bessel_i2n: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_i2n; + + KYOSU_DEFERS_CALLABLE(sph_bessel_i2n_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& z) noexcept + { + auto fn = callable_sph_bessel_i2n{}; + return fn(n, complex(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 sph_bessel_i2n +//! @brief Computes the spherical Bessel functions +//! \f$ i_n^{(2)}(z) = i^{-n-1}y_n(iz)\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_i2n(int n, T z) noexcept; +//! template constexpr auto sph_bessel_i2n(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$i_n^{(2)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_i2n.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_i2n sph_bessel_i2n = {}; +} diff --git a/include/kyosu/functions/sph_bessel_j0.hpp b/include/kyosu/functions/sph_bessel_j0.hpp new file mode 100644 index 00000000..6f03d936 --- /dev/null +++ b/include/kyosu/functions/sph_bessel_j0.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_sph_bessel_j0: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_j0; + + KYOSU_DEFERS_CALLABLE(sph_bessel_j0_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::sph_bessel_j0(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 sph_bessel_j0 +//! @brief Computes the Bessel function of the first kind, +//! \f$ j_0(x)=\sin z/z \f$ extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_j0(T z) noexcept; +//! template constexpr T sph_bessel_j0(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$j_0(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_j0.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_j0 sph_bessel_j0 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_j1.hpp b/include/kyosu/functions/sph_bessel_j1.hpp new file mode 100644 index 00000000..0b4c3c2c --- /dev/null +++ b/include/kyosu/functions/sph_bessel_j1.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_sph_bessel_j1: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_j1; + + KYOSU_DEFERS_CALLABLE(sph_bessel_j1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::sph_bessel_j1(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 sph_bessel_j1 +//! @brief Computes the Spherical Bessel function of the first kind, +//! \f$ j_1(x)= \sin z/z^2 -\cos z/zf$ +//! extended to the complex plane and cayley_dickson values. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_j1(T z) noexcept; +//! template constexpr T sph_bessel_j1(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$j_1(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_j1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_j1 sph_bessel_j1 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_jn.hpp b/include/kyosu/functions/sph_bessel_jn.hpp new file mode 100644 index 00000000..ea147ccb --- /dev/null +++ b/include/kyosu/functions/sph_bessel_jn.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_sph_bessel_jn: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_jn; + + KYOSU_DEFERS_CALLABLE(sph_bessel_jn_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& v) noexcept { return eve::sph_bessel_jn(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 sph_bessel_jn +//! @brief Computes the spherical Bessel functions of the first kind \f$j_{n}(x)\f$, +//! extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_jn(int n, T z) noexcept; +//! template constexpr T sph_bessel_jn(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$j_n(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_jn.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_jn sph_bessel_jn = {}; +} diff --git a/include/kyosu/functions/sph_bessel_k0.hpp b/include/kyosu/functions/sph_bessel_k0.hpp new file mode 100644 index 00000000..78f59f96 --- /dev/null +++ b/include/kyosu/functions/sph_bessel_k0.hpp @@ -0,0 +1,76 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_sph_bessel_k0: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_k0; + + KYOSU_DEFERS_CALLABLE(sph_bessel_k0_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept + { + return eve::pio_2(eve::as(z))*exp(-z)/z; + } + + 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 sph_bessel_k0 +//! @brief Computes the spherical Bessel function of the first kind, +//! \f$ k_0(x)= \pi e^{-z}/(2z) \f$ extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_k0(T z) noexcept; +//! template constexpr T sph_bessel_k0(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$k_0(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_k0.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_k0 sph_bessel_k0 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_k1.hpp b/include/kyosu/functions/sph_bessel_k1.hpp new file mode 100644 index 00000000..43d3134b --- /dev/null +++ b/include/kyosu/functions/sph_bessel_k1.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_sph_bessel_k1: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_k1; + + KYOSU_DEFERS_CALLABLE(sph_bessel_k1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& z) noexcept + { + auto rz = rec(z); + return eve::pio_2(eve::as(z))*exp(-z)*fma(rz, rz, rz); + } + + 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 sph_bessel_k1 +//! @brief Computes the spherical Bessel function of the first kind, +//! \f$ k_1(x)= (\pi/2) e^{-z}(1/z+1/z^2)\f$ extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_k1(T z) noexcept; +//! template constexpr T sph_bessel_k1(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$k_1(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_k1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_k1 sph_bessel_k1 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_kn.hpp b/include/kyosu/functions/sph_bessel_kn.hpp new file mode 100644 index 00000000..f4aff271 --- /dev/null +++ b/include/kyosu/functions/sph_bessel_kn.hpp @@ -0,0 +1,79 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_sph_bessel_kn: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_kn; + + KYOSU_DEFERS_CALLABLE(sph_bessel_kn_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, int n, T const& v) noexcept + { + callable_sph_bessel_kn fn{}; + return real(fn(n, complex(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 sph_bessel_kn +//! @brief Computes the spherical Bessel functions \f$k_{n}(x)\f$, +//! extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_kn(int n, T z) noexcept; +//! template constexpr T sph_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/sph_bessel_kn.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_kn sph_bessel_kn = {}; +} diff --git a/include/kyosu/functions/sph_bessel_y0.hpp b/include/kyosu/functions/sph_bessel_y0.hpp new file mode 100644 index 00000000..b1ff9c46 --- /dev/null +++ b/include/kyosu/functions/sph_bessel_y0.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_sph_bessel_y0: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_y0; + + KYOSU_DEFERS_CALLABLE(sph_bessel_y0_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::sph_bessel_y0(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 sph_bessel_y0 +//! @brief Computes the spherical Bessel function of the first kind, +//! \f$ y_0(x)=\cos z/z \f$ extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_y0(T z) noexcept; +//! template constexpr T sph_bessel_y0(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$y_0(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_y0.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_y0 sph_bessel_y0 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_y1.hpp b/include/kyosu/functions/sph_bessel_y1.hpp new file mode 100644 index 00000000..02410ef2 --- /dev/null +++ b/include/kyosu/functions/sph_bessel_y1.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_sph_bessel_y1: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_y1; + + KYOSU_DEFERS_CALLABLE(sph_bessel_y1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::sph_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 sph_bessel_y1 +//! @brief Computes the spherical Bessel function of the first kind, +//! \f$ y_1(x)= -\coz z/z^2 -\sin z/zf$ +//! extended to the complex plane and cayley_dickson values. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_y1(T z) noexcept; +//! template constexpr T sph_bessel_y1(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$y_1(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_y1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_y1 sph_bessel_y1 = {}; +} diff --git a/include/kyosu/functions/sph_bessel_yn.hpp b/include/kyosu/functions/sph_bessel_yn.hpp new file mode 100644 index 00000000..cb43f581 --- /dev/null +++ b/include/kyosu/functions/sph_bessel_yn.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_sph_bessel_yn: eve::elementwise + { + using callable_tag_type = callable_sph_bessel_yn; + + KYOSU_DEFERS_CALLABLE(sph_bessel_yn_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, int n, T const& v) noexcept { return eve::sph_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 sph_bessel_yn +//! @brief Computes the spherical Bessel functions of the first kind \f$y_{n}(x)\f$, +//! extended to the complex plane and cayley_dickson algebras. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto sph_bessel_yn(int n, T z) noexcept; +//! template constexpr T sph_bessel_yn(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * returns \f$y_n(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/sph_bessel_yn.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_sph_bessel_yn sph_bessel_yn = {}; +} diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp index 01cb5975..2f061e57 100644 --- a/include/kyosu/types/impl/bessel.hpp +++ b/include/kyosu/types/impl/bessel.hpp @@ -11,16 +11,59 @@ #include #include #include +//====================================================================================================================== +//! Up to now bessel functions of integer orders are only implemented here for scalar integral orders +//! +//! cylindical functions +//! +//! * \f$J_n\f$ many foloow the design of Jhonas Olivati de Sarro, with adaptations to simd cases. +//! Its copyrigth notice is reproduced at the end of this file +//! * \f$I_n\f$ are obtained from \f$J_n\f$ by various connections formulas that can be found in DLMF. +//! * \f$Y_0\f$ is obtained using the formula \f$Y_0(z) = \frac2\pi\left[\log(z/2)+\gamma)\right] J_0(z) - \frac4\pi\sum_{k = 1}_\infty (-1)^k \frac{J_{2k}}k \f$. +//! * \f$Y_n\f$ is then obtained by backward recurrence from the \f$(J_i(z))_{i \le n}\f$ and \f$Y_0\f$ +//! * \f$H_n^{(1)}\f$ is \f$J_n+ i Y_n\f$ +//! * \f$H_n^{(2)}\f$ is \f$J_n- i Y_n\f$ +//! * \f$K_n\f$ is computed from \f$H_n^{(1)}\f$ and \f$H_n^{(2)}\f$ with DLMF 10.27.8 +//! +//! spherical functions +//! +//! As in DMLF we restrict to non-negative order +//! +//! The computation scheme of \f$j_n\f$, \f$y_n\f$, \f$h_n^{(1)}\f$ and \f$h_n^{(2)}\f$ follows the article of Liang-Wu Kai +//! (On the computation of spherical Bessel functions of complex arguments, +//! Comput. Phys. Commun., 182 (3) (2011), pp. 663-668)). +//! \f$k_n\f$ is computed from f$h_n^{(1)}\f$ using DLMF 10.47.13 +//! * \f$i_n^{(1)}\f$ and \f$i_n^{(2)}\f$ are given by 10.47.12 + #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include + +#include +#include +#include +#include +#include + +//================================================================================ +// MIT License + +// Copyright (c) 2021 Jhonas Olivati de Sarro + +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: + +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. + +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//================================================================================ diff --git a/include/kyosu/types/impl/detail/bessel_h.hpp b/include/kyosu/types/impl/detail/bessel_h.hpp new file mode 100644 index 00000000..0f3173aa --- /dev/null +++ b/include/kyosu/types/impl/detail/bessel_h.hpp @@ -0,0 +1,151 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include + +namespace kyosu::_ +{ + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + // cylindrical bessel of the third kind + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + + //===------------------------------------------------------------------------------------------- + // 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); + } + } + + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + // spherical bessel of the third kind + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + + + //===------------------------------------------------------------------------------------------- + // sph_bessel_h1n + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, int n, Z z) noexcept + { + if constexpr(concepts::complex) + { + if (n < 0) + { + return eve::nan(as(z)); + } + else + { + using u_t = eve::underlying_type_t; + return eve::sign_alternate(u_t(n))*sph_bessel_h2n(n, -z); + } + } + else + { + return cayley_extend_rev(sph_bessel_h1n, n, z); + } + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return sph_bessel_h1n(0, z); + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return sph_bessel_h1n(1, z); + } + + //===------------------------------------------------------------------------------------------- + // sph_bessel_h2n + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, N n, Z z) noexcept + { + if constexpr(concepts::complex) + { + if (n < 0) + { + return eve::nan(as(z)); + } + else + { + auto imzlt0 = eve::is_ltz(imag(z)); + z = if_else(imzlt0, z, -z); + using u_t = eve::underlying_type_t; + auto i = complex(u_t(0), u_t(1)); + auto rz = rec(z); + auto miz = complex(ipart(z), -real(z)); + auto h0 = if_else(imzlt0, i*exp(miz), -i*exp(-miz))*rz; + if(n == 0) return h0; + auto h1 = if_else(imzlt0, (rz+i),(rz-i))*h0 ; + if(n == 1) return if_else(imzlt0,h1, -h1); + + auto h2 = h1; + for(int i=1; i < n; ++i) + { + auto f = (2*i+1)*rz; + h2 = f*h1-h0; + h0 = h1; + h1 = h2; + } + return if_else(imzlt0, h2, eve::sign_alternate(u_t(n))*h2); + } + } + else + { + return cayley_extend_rev(sph_bessel_h2n, n, z); + } + } + + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return sph_bessel_h2n(0, z); + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return sph_bessel_h2n(1, z); + } + + +} diff --git a/include/kyosu/types/impl/detail/bessel_i.hpp b/include/kyosu/types/impl/detail/bessel_i.hpp new file mode 100644 index 00000000..e431ccdd --- /dev/null +++ b/include/kyosu/types/impl/detail/bessel_i.hpp @@ -0,0 +1,182 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_i0 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return cyl_bessel_j0(complex(-ipart(z), real(z))); + } + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_i1 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + auto [r, i] = cyl_bessel_j1(complex(-ipart(z), real(z))); + return complex(i, -r); + } + else + { + return cayley_extend(cyl_bessel_i1, z); + } + } + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_in + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, N n, Z z) noexcept + { + if constexpr(concepts::complex ) + { + using e_t = as_real_type_t; + auto miton = [](N n){ + if (n%4 == 0) return complex(eve::one(eve::as())); + else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::mone(eve::as())); + else if (n%4 == 2) return complex(eve::mone(eve::as())); + else return complex(eve::zero(eve::as()), eve::one(eve::as())); + }; + auto an = eve::abs(n); + return miton(an)*cyl_bessel_jn(an,complex(-ipart(z), real(z))); + } + else + { + return cayley_extend_rev(cyl_bessel_in, n, z); + } + } + + //===------------------------------------------------------------------------------------------- + // spherical besseli + //===------------------------------------------------------------------------------------------- + // spheical i split i,n i1_n and i2_n + //===------------------------------------------------------------------------------------------- + + //===------------------------------------------------------------------------------------------- + // sph_bessel_i1_0 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return sinhc(z); + } + else + { + return cayley_extend(sph_bessel_i1_0, z); + } + } + + //===------------------------------------------------------------------------------------------- + // sph_bessel_i1_1 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + auto [sh, ch] = sinhcosh(z); + auto rz = rec(z); + return fnma(sh, rz, ch)*rz; + } + else + { + return cayley_extend(sph_bessel_i1_1, z); + } + } + + //===------------------------------------------------------------------------------------------- + // sph_bessel_i1n + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, N n, Z z) noexcept + { + if constexpr(concepts::complex ) + { + using e_t = as_real_type_t; + auto iton = [](N n){ + if (n%4 == 0) return complex(eve::one(eve::as())); + else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::one(eve::as())); + else if (n%4 == 2) return complex(eve::mone(eve::as())); + else return complex(eve::zero(eve::as()), eve::mone(eve::as())); + }; + return iton(-n)*sph_bessel_jn(n,complex(-ipart(z), real(z))); + } + else + { + return cayley_extend_rev(sph_bessel_i1n, n, z); + } + } + + //===------------------------------------------------------------------------------------------- + // sph_bessel_i2_0 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return cosh(z)/z; + } + else + { + return cayley_extend(sph_bessel_i2_0, z); + } + } + + //===------------------------------------------------------------------------------------------- + // sph_bessel_i2 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + auto [sh, ch] = sinhcosh(z); + auto rz = rec(z); + return eve::fnma(ch, rz, sh)*rz; + } + else + { + return cayley_extend(sph_bessel_i2_1, z); + } + } + + //===------------------------------------------------------------------------------------------- + // sph_bessel_i2n + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, N n, Z z) noexcept + { + if constexpr(concepts::complex ) + { + using e_t = as_real_type_t; + auto iton = [](N n){ + if (n%4 == 0) return complex(eve::one(eve::as())); + else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::one(eve::as())); + else if (n%4 == 2) return complex(eve::mone(eve::as())); + else return complex(eve::zero(eve::as()), eve::mone(eve::as())); + }; + return iton(-n-1)*sph_bessel_yn(n,complex(-ipart(z), real(z))); + } + else + { + return cayley_extend_rev(sph_bessel_i2n, n, z); + } + } +} diff --git a/include/kyosu/types/impl/detail/bessel_j.hpp b/include/kyosu/types/impl/detail/bessel_j.hpp new file mode 100644 index 00000000..4acb0c93 --- /dev/null +++ b/include/kyosu/types/impl/detail/bessel_j.hpp @@ -0,0 +1,465 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include + +namespace kyosu::_ +{ + //===------------------------------------------------------------------------------------------- + // cyl_bessel_j0 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + using e_t = as_real_type_t; + using u_t = eve::underlying_type_t; + auto saz = kyosu::sqr_abs(z); + + auto ascending_series_cyl_j0 = [](auto z) + { + // Ascending Series from G. N. Watson 'A treatise on the + // theory of Bessel functions', 2ed, Cambridge, 1996, + // Chapter II, equation (3); or from Equation 9.1.12 of + // M. Abramowitz, I. A. Stegun 'Handbook of Mathematical + // Functions'. + // good for abs(z) < 12 + auto eps2 = eve::sqr(eve::eps(eve::as())); + auto j0 = complex(eve::one((eve::as()))); + auto sm = j0; + auto test = sqr_abs(sm) >= eps2*sqr_abs(j0); + auto m(eve::one(eve::as())); + auto qz2 = - sqr(z)*u_t(0.25); + size_t im = 1; + while(eve::any(test)) + { + sm *= qz2/eve::sqr(im); + j0 += sm; + test = sqr_abs(sm) >= eps2*sqr_abs(j0); + ++im; + } + return j0; + }; + + auto semiconvergent_series_cyl_j0 = [saz](auto z) { + auto bound_compute = [saz]() + { + auto bds = eve::if_else(saz < 50*50, e_t(10), e_t(8)); + bds = eve::if_else(saz < 35*35, e_t(12), bds); + bds = eve::if_else(saz < 12*12, e_t(0), bds); + return bds; + }; + + // DLMF 10.17.1 Hankel expansion tabulated p[i] = a[2*i], q[i] = a{2*i+1] + std::array pim{1.0000000000000000e+00, + 7.0312500000000000e-02, 1.5950520833333335e+00, 5.1046874999999998e+00, 1.0609654017857142e+01, + 1.8112673611111113e+01, 2.7614701704545457e+01, 3.9116157280219781e+01, 5.2617252604166673e+01, + 6.8118106617647058e+01, 8.5618791118421058e+01, 1.0511935200216450e+02, 1.2661981997282609e+02, + 1.5012021634615385e+02, 1.7562055638227514e+02 + }; + std::array xim{1.2500000000000000e-01, + 5.8593750000000000e-01, 3.1007812499999998e+00, 7.6075148809523814e+00, 1.4111328125000002e+01, + 2.2613778409090912e+01, 3.3115484775641029e+01, 4.5616741071428571e+01, 6.0117704503676471e+01, + 7.6618466739766092e+01, 9.5119084821428572e+01, 1.1561959609683794e+02, 1.3812002604166668e+02, + 1.6262039262820511e+02, 1.8912070889778326e+02 + }; + auto rz = rec(z); + auto rz2= sqr(rz); + auto Pm = complex(e_t(pim[0])); + auto Qm = rz*xim[0]; + auto P = Pm; + auto Q = Qm; + auto bds = bound_compute(); + + size_t im = 1; + auto bound_not_reached = u_t(im) <= bds; + while (eve::any(bound_not_reached)) + { + Pm *= -pim[im]*rz2; + Qm *= -xim[im]*rz2; + if constexpr(eve::scalar_value) + { + P += Pm; + Q += Qm; + } + else + { + P += if_else( bound_not_reached, Pm, eve::zero); + Q += if_else( bound_not_reached, Qm, eve::zero); + } + bound_not_reached = u_t(++im) <= bds; + } + auto [s, c] = kyosu::sincos(z); + u_t rsqrtpi = eve::rsqrt_pi(eve::as()); + auto r = rsqrtpi*(c*(P-Q) + s*(P+Q))/(kyosu::sqrt(z)); + return r; + }; + auto izneg = eve::is_ltz(imag(z)); + z = if_else(izneg, conj(z), z); + auto rzneg = eve::is_ltz(real(z)); + z = if_else(rzneg, -z, z); + + auto r = kyosu::if_else(is_eqz(saz), Z(1), eve::nan(eve::as(saz))); + auto notdone = kyosu::is_nan(r); + + if( eve::any(notdone) ) + { + notdone = next_interval(ascending_series_cyl_j0, notdone, saz <= as_real_type_t(144), r, z); + if( eve::any(notdone) ) + { + last_interval(semiconvergent_series_cyl_j0, notdone, r, z); + } + } + imag(r) = eve::if_else(is_real(z), eve::zero, imag(r)); + r = if_else(izneg, conj(r), r); + return r; + } + else + { + return cayley_extend(cyl_bessel_j0, z); + } + } + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_j1 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + using e_t = as_real_type_t; + using u_t = eve::underlying_type_t; + auto saz = kyosu::sqr_abs(z); + + auto ascending_series_cyl_j1 = [](auto z) { + // Ascending Series from G. N. Watson 'A treatise on the + // theory of Bessel functions', 2ed, Cambridge, 1996, + // Chapter II, equation (3); or from Equation 9.1.12 of + // M. Abramowitz, I. A. Stegun 'Handbook of Mathematical + // Functions'. + // good for abs(z) < 12 + auto eps2 = sqr(eve::eps(eve::as())); + auto j1 = complex(eve::one((eve::as()))); + auto sm = j1; + auto test = sqr_abs(sm) >= eps2*sqr_abs(j1); + auto m(eve::one(eve::as())); + auto z2o_4 = - sqr(z)*e_t(0.25); + while(eve::any(test)) + { + sm *= z2o_4/(m*inc(m)); + j1 += sm; + test = sqr_abs(sm) >= eps2*sqr_abs(j1); + m = inc(m); + } + return j1*eve::half(eve::as())*z; + }; + + auto semiconvergent_series_cyl_j1 = [saz](auto z) + { + auto bound_compute = [saz]() + { + auto bds = eve::if_else(saz < 50*50, e_t(10), e_t(8)); + bds = eve::if_else(saz < 35*35, e_t(12), bds); + bds = eve::if_else(saz < 12*12, e_t(0), bds); + return bds; + }; + + // DLMF 10.17.1 Hankel expansion tabulated p[i] = a[2*i], q[i] = a{2*i+1] + std::array pim{1.0000000000000000e+00, + -1.1718750000000000e-01, 1.2304687500000000e+00, 4.6921875000000002e+00, 1.0174386160714286e+01, + 1.7664062500000004e+01, 2.7157315340909090e+01, 3.8652558379120876e+01, 5.2149023437499999e+01, + 6.7646292892156850e+01, 8.5144120065789465e+01, 1.0464234983766234e+02, 1.2614087975543477e+02, + 1.4963963942307691e+02, 1.7513857886904762e+02 + }; + std::array xim{-3.7500000000000000e-01, + 2.7343750000000000e-01, 2.7070312500000000e+00, 7.1819196428571432e+00, 1.3668619791666668e+01, + 2.2160369318181818e+01, 3.2654747596153847e+01, 4.5150669642857139e+01, 5.9647575827205877e+01, + 7.6145148026315780e+01, 9.4643191964285705e+01, 1.1514158226284583e+02, 1.3764023437500001e+02, + 1.6213908920940170e+02, 1.8863810421798030e+02 + }; + auto rz = rec(z); + auto rz2= sqr(rz); + auto Pm = complex(e_t(pim[0])); + auto Qm = complex(xim[0])*rz; + auto P = Pm; + auto Q = Qm; + auto bds = bound_compute(); + size_t im = 1; + auto bound_not_reached = u_t(im) <= bds; + + while (eve::any(bound_not_reached)) + { + auto m = u_t(im); + Pm *= -pim[im]*rz2; + Qm *= -xim[im]*rz2; + if constexpr(eve::scalar_value) + { + P += Pm; + Q += Qm; + } + else + { + P += kyosu::if_else(bound_not_reached, Pm, eve::zero); + Q += kyosu::if_else(bound_not_reached, Qm, eve::zero); + } + bound_not_reached = u_t(++im) <= bds; + } + auto [s, c] = kyosu::sincos(z); + u_t rsqrtpi = eve::rsqrt_pi(eve::as()); + auto r = rsqrtpi*(s*(P-Q) - c*(P+Q))/(kyosu::sqrt(z)); + return r; + }; + + auto izneg = eve::is_ltz(imag(z)); + z = if_else(izneg, conj(z), z); + auto rzneg = eve::is_ltz(real(z)); + z = if_else(rzneg, -z, z); + + auto r = kyosu::if_else(is_eqz(saz), Z(1), eve::nan(eve::as(saz))); + auto notdone = kyosu::is_nan(r); + + if( eve::any(notdone) ) + { + notdone = next_interval(ascending_series_cyl_j1, notdone, saz <= as_real_type_t(144), r, z); + if( eve::any(notdone) ) + { + last_interval(semiconvergent_series_cyl_j1, notdone, r, z); + } + } + real(r) = eve::if_else(is_pure(z), eve::zero, real(r)); + r = if_else(rzneg, -r, r); + return if_else(izneg, conj(r), r); + } + else + { + return cayley_extend(cyl_bessel_j1, z); + } + } + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_jn + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, N nn, Z z) noexcept + { + if constexpr(concepts::complex ) + { + if ( is_eqz(nn) ) + { + return cyl_bessel_j0(z); + } + else if (nn == 1) + { + return cyl_bessel_j1(z); + } + else if ( nn == -1 ) + { + return -cyl_bessel_j1(z); + } + else + { + using e_t = as_real_type_t; + using u_t = eve::underlying_type_t; + auto n = u_t(eve::abs(nn)); + auto az = kyosu::abs(z); + + auto forward = [n](auto z){ + auto b0 = cyl_bessel_j0(z); + auto b1 = cyl_bessel_j1(z); + Z bn; + auto rz = rec(z); + for ( int k=1; kINF) ). + return u_t(0.5)*eve::log10(u_t(6.28)*n) - n*eve::log10(u_t(1.36)*az/n); + }; + + auto ini_for_br_1 = [minus_log10_cyl_j_at_infinity](auto az, auto mg) + { + // Starting point for backward recurrence + // for when |Jn(x)|~10e-mg + // using the secant method. + 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; + auto nn = n1 - (n1 - n0)/oneminus(f0/f1); + auto f = minus_log10_cyl_j_at_infinity(nn, az) - mg; + auto test = eve::abs(nn - n1) > 1; + while ( eve::any(test) ) + { + n0 = n1; + f0 = f1; + n1 = nn; + f1 = f; + nn = eve::if_else(test, n1 - (n1 - n0)/oneminus(f0/f1), nn); + f = eve::if_else(test, minus_log10_cyl_j_at_infinity(nn, az) - mg, f); + test = eve::abs(nn - n1) > 1; + } + return eve::trunc(nn); + }; + + auto ini_for_br_2 = [minus_log10_cyl_j_at_infinity](auto n, auto az, auto sd){ + // Starting point for backward recurrence + // for when Jn(x) has sd significant digits + // using the secant method. + auto hmp = eve::half(eve::as())*sd; + 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::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; + auto nn = n1 - (n1-n0)/oneminus(f0/f1); + auto f = minus_log10_cyl_j_at_infinity(nn, az) - obj; + auto test = eve::abs(nn - n1) >= 1; + while ( eve::any(test)) + { + n0 = n1; + f0 = f1; + n1 = nn; + f1 = f; + nn = eve::if_else(test, n1 - (n1-n0)/(oneminus(f0/f1)), nn); + f = eve::if_else(test, minus_log10_cyl_j_at_infinity(nn, az) - obj, f); + test = eve::abs(nn - n1) >= 1; + } + return eve::trunc(nn + 10); + }; + + auto backward = [az, n, ini_for_br_1, ini_for_br_2](auto z){ + 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); + Z bn(cf); + auto k = m; + auto kgez = eve::is_gez(k); + while (eve::any(kgez)) + { + cf = kyosu::if_else(kgez, 2*inc(k)*cf1*rec(z)-cf2, cf); + bn = kyosu::if_else ( k == n, cf, bn); + cf2 = kyosu::if_else(kgez, cf1, cf2); + cf1 = kyosu::if_else(kgez, cf, cf1); + k = dec(k); + kgez = eve::is_gez(k); + } + auto j0 = cyl_bessel_j0(z); + auto j1 = cyl_bessel_j1(z); + bn *= eve::if_else ( sqr_abs(j0) > sqr_abs(j1), j0/cf, j1/cf2); + + return bn; + }; + + auto srz = eve::signnz(real(z)); + z *= srz; + + auto r = kyosu::if_else(is_eqz(az), Z(0), eve::nan(eve::as(az))); + auto notdone = kyosu::is_nan(r); + if( eve::any(notdone) ) + { + notdone = next_interval(forward, notdone, 4*n < az, r, z); + if( eve::any(notdone) ) + { + last_interval(backward, notdone, r, z); + } + } + auto sgnaltern = [n](auto x){return eve::if_else(eve::is_ltz(x), eve::one, eve::sign_alternate(n));}; + r = sgnaltern(srz)*sgnaltern(n)*r; + return nn < 0 ? r*eve::sign_alternate(u_t(nn)) : r; + } + } + else + { + return cayley_extend_rev(cyl_bessel_jn, nn, z); + } + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return sinc(z); + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + auto rz = rec(z); + return if_else(kyosu::abs(z) < eve::eps(eve::as(real(z))) + , eve::zero + , (sinc(z)-cos(z))*rz + ); + } + else + { + return cayley_extend(sph_bessel_j0, z); + } + } + + template + auto dispatch(eve::tag_of, N n, Z z) noexcept + { + if constexpr(concepts::complex ) + { + using u_t = eve::underlying_type_t; + auto bd = [z](int n){ + auto st = eve::abs(eve::sin(eve::abs(arg(z)))); + auto r = kyosu::abs(z); + auto m = eve::maximum(eve::ceil((u_t(1.83)+u_t(4.1)*eve::pow(st, u_t(0.33)))* + eve::pow(r, (u_t(0.91)-u_t(0.43)*eve::pow(st, u_t(0.36))))+9*(1-eve::sqrt(st)))); + auto nn = (eve::any(is_real(z))) ? n+5 : n; + return eve::max(nn, int(inc(m))); + }; + + auto j0 = kyosu::sph_bessel_j0(z); + if (n == 0) return j0; + auto j1 = kyosu::sph_bessel_j1(z); + if (n == 1) return j1; + + auto rz = kyosu::rec(z); + auto nn = bd(n); + Z jnext(kyosu::complex(u_t(0))); + Z j(kyosu::complex(eve::sqrtsmallestposval(eve::as()))); + auto init = j; + auto jcur = jnext; + auto res = j; + for(int i=nn-1; i > 0 ; --i) + { + jcur = (2*i+1)*rz*j-jnext; + if(i == n) res = j; + jnext = j; + j = jcur; + } + auto j0ltj1 = kyosu::abs(j0) <= kyosu::abs(j1); + auto scalej0 = (j0/jcur); + auto scalej1 = (j1/jnext); + res *= if_else(j0ltj1, scalej0, scalej1); + return if_else(is_eqz(z), eve::zero, res); + } + else + { + return cayley_extend_rev(sph_bessel_jn, n, z); + } + } + +} diff --git a/include/kyosu/types/impl/detail/besselkn.hpp b/include/kyosu/types/impl/detail/bessel_k.hpp similarity index 53% rename from include/kyosu/types/impl/detail/besselkn.hpp rename to include/kyosu/types/impl/detail/bessel_k.hpp index f9bf7b4c..043a2138 100644 --- a/include/kyosu/types/impl/detail/besselkn.hpp +++ b/include/kyosu/types/impl/detail/bessel_k.hpp @@ -6,7 +6,7 @@ */ //====================================================================================================================== #pragma once -#include +#include namespace kyosu::_ { @@ -28,7 +28,8 @@ namespace kyosu::_ 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) + auto argzlt0 = eve::is_ltz(argz); + auto r = if_else(argzlt0 , 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); @@ -50,4 +51,46 @@ namespace kyosu::_ { return cyl_bessel_kn(1, z); } + + + //===------------------------------------------------------------------------------------------- + // sph_bessel_kn + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, int n, Z z) noexcept + { + using u_t = eve::underlying_type_t; + if constexpr(concepts::complex) + { + auto iton = [](int n){ + if (n%4 == 0) return complex(eve::one(eve::as())); + else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::one(eve::as())); + else if (n%4 == 2) return complex(eve::mone(eve::as())); + else return complex(eve::zero(eve::as()), eve::mone(eve::as())); + }; + return -eve::pio_2(eve::as())*iton(n)*sph_bessel_h1n(n, complex(-imag(z), real(z))); + } + else + { + return cayley_extend_rev(sph_bessel_kn, n, z); + } + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + using u_t = eve::underlying_type_t; + return eve::pio_2(eve::as())*exp(-z)/z; + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + using u_t = eve::underlying_type_t; + auto rz = rec(z); + return eve::pio_2(eve::as())*exp(-z)*fma(rz, rz, rz); + } + + + } diff --git a/include/kyosu/types/impl/detail/bessel_y.hpp b/include/kyosu/types/impl/detail/bessel_y.hpp new file mode 100644 index 00000000..f3e135dd --- /dev/null +++ b/include/kyosu/types/impl/detail/bessel_y.hpp @@ -0,0 +1,169 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include + +namespace kyosu::_ +{ + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + // utilities + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + + template < typename Z> struct mkjs + { + mkjs(size_t n, Z z): + rz(kyosu::rec(z)), + rs(kyosu::_:: R(n, z)), + j(kyosu::cyl_bessel_jn(n, z)), + i(n-1){} + + auto operator()(){ + auto pj = j; + j = pj*rs; + rs = 2*(i--)*rz-kyosu::rec(rs); + return pj; + } + + Z rz, rs, j, pj; + int i; + }; + + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + // cylindrical bessel + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + + //===------------------------------------------------------------------------------------------- + // 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 bd = bound(z); + Z s{}; + mkjs jj(2*bd-2, z); + auto sgn = eve::sign_alternate(u_t(bd+1)); + for (int k = bd-1; k >= 1; --k) + { + auto sk = sgn*jj()/k; + jj(); + sgn = -sgn; + s+= sk; + } + return if_else(is_eqz(z), complex(eve::minf(eve::as())), twoopi*((log(z/2)+egamma)*jj()-2*s)); + } + + + //===------------------------------------------------------------------------------------------- + // 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); + } + + + //===------------------------------------------------------------------------------------------- + // 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) + { + int n = eve::abs(nn); + auto do_it = [n, z, &y](auto js){ + mkjs jj(n, z); + for(int i=n; i >= 0 ; --i) js[i] = jj(); + using u_t = eve::underlying_type_t; + auto twoopi = eve::two_o_pi(eve::as()); + auto b = twoopi*rec(z); + for(int i=1; i <= n ; ++i) y = fms(js[i], y, b)/js[i-1]; + }; + kyosu::_::with_alloca(n+1, do_it); + 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; + } + + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + // spherical bessel + //===------------------------------------------------------------------------------------------- + //===------------------------------------------------------------------------------------------- + + //===------------------------------------------------------------------------------------------- + // sph_bessel_y0 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return -cos(z)/z; + } + + //===------------------------------------------------------------------------------------------- + // sph_bessel_y1 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return -cos(z)/sqr(z)-sinc(z); + } + + + //===------------------------------------------------------------------------------------------- + // sph_bessel_yn + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, int n, Z z) noexcept + { + if constexpr(concepts::complex ) + { + using u_t = eve::underlying_type_t; + auto imzge0 = eve::is_gez(imag(z)); + auto jn = sph_bessel_jn(n, z); + auto i = complex(u_t(0), u_t(1)); + Z res; + if (eve::any(imzge0)) + { + auto z1 = if_else(imzge0, z, conj(z)); + res = jn-sph_bessel_h1n(n, z1); + } + if (!eve::all(imzge0)) + { + auto z2 = if_else(imzge0, conj(z), z); + res = if_else(imzge0, res, sph_bessel_h2n(n, z2)-jn); + } + return complex(-imag(res), real(res)); + } + else + { + return cayley_extend_rev(sph_bessel_yn, n, z); + } + } + + +} diff --git a/include/kyosu/types/impl/detail/besselhn.hpp b/include/kyosu/types/impl/detail/besselhn.hpp deleted file mode 100644 index 95ddd8be..00000000 --- a/include/kyosu/types/impl/detail/besselhn.hpp +++ /dev/null @@ -1,46 +0,0 @@ -//====================================================================================================================== -/* - 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/besseli1.hpp b/include/kyosu/types/impl/detail/besseli1.hpp deleted file mode 100644 index e6e1523f..00000000 --- a/include/kyosu/types/impl/detail/besseli1.hpp +++ /dev/null @@ -1,30 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include - -namespace kyosu::_ -{ - - //===------------------------------------------------------------------------------------------- - // cyl_bessel_i1 - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, Z z) noexcept - { - if constexpr(concepts::complex ) - { - auto [r, i] = cyl_bessel_j1(complex(-ipart(z), real(z))); - return complex(i, -r); - } - else - { - return cayley_extend(cyl_bessel_i1, z); - } - } -} diff --git a/include/kyosu/types/impl/detail/besselin.hpp b/include/kyosu/types/impl/detail/besselin.hpp deleted file mode 100644 index 925b7595..00000000 --- a/include/kyosu/types/impl/detail/besselin.hpp +++ /dev/null @@ -1,37 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include - -namespace kyosu::_ -{ - - //===------------------------------------------------------------------------------------------- - // cyl_bessel_in - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, N n, Z z) noexcept - { - if constexpr(concepts::complex ) - { - using e_t = as_real_type_t; - auto miton = [](N n){ - if (n%4 == 0) return complex(eve::one(eve::as())); - else if (n%4 == 1) return complex(eve::zero(eve::as()), eve::mone(eve::as())); - else if (n%4 == 2) return complex(eve::mone(eve::as())); - else return complex(eve::zero(eve::as()), eve::one(eve::as())); - }; - auto an = eve::abs(n); - return miton(an)*cyl_bessel_jn(an,complex(-ipart(z), real(z))); - } - else - { - return cayley_extend_rev(cyl_bessel_in, n, z); - } - } -} diff --git a/include/kyosu/types/impl/detail/besselj0.hpp b/include/kyosu/types/impl/detail/besselj0.hpp index 0fc16f93..5c499637 100644 --- a/include/kyosu/types/impl/detail/besselj0.hpp +++ b/include/kyosu/types/impl/detail/besselj0.hpp @@ -76,7 +76,6 @@ namespace kyosu::_ auto Q = Qm; auto bds = bound_compute(); - Z zero{}; size_t im = 1; auto bound_not_reached = u_t(im) <= bds; while (eve::any(bound_not_reached)) @@ -90,8 +89,8 @@ namespace kyosu::_ } else { - P += if_else( bound_not_reached, Pm, zero); - Q += if_else( bound_not_reached, Qm, zero); + P += if_else( bound_not_reached, Pm, eve::zero); + Q += if_else( bound_not_reached, Qm, eve::zero); } bound_not_reached = u_t(++im) <= bds; } diff --git a/include/kyosu/types/impl/detail/besselj1.hpp b/include/kyosu/types/impl/detail/besselj1.hpp index 599d0aaf..36f48aa5 100644 --- a/include/kyosu/types/impl/detail/besselj1.hpp +++ b/include/kyosu/types/impl/detail/besselj1.hpp @@ -74,7 +74,6 @@ namespace kyosu::_ auto P = Pm; auto Q = Qm; auto bds = bound_compute(); - auto zero = Z{}; size_t im = 1; auto bound_not_reached = u_t(im) <= bds; @@ -90,8 +89,8 @@ namespace kyosu::_ } else { - P += kyosu::if_else(bound_not_reached, Pm, zero); - Q += kyosu::if_else(bound_not_reached, Qm, zero); + P += kyosu::if_else(bound_not_reached, Pm, eve::zero); + Q += kyosu::if_else(bound_not_reached, Qm, eve::zero); } bound_not_reached = u_t(++im) <= bds; } diff --git a/include/kyosu/types/impl/detail/besseljn.hpp b/include/kyosu/types/impl/detail/besseljn.hpp index fd9b484d..a7758a8e 100644 --- a/include/kyosu/types/impl/detail/besseljn.hpp +++ b/include/kyosu/types/impl/detail/besseljn.hpp @@ -6,8 +6,6 @@ */ //====================================================================================================================== #pragma once -#include -#include namespace kyosu::_ { diff --git a/include/kyosu/types/impl/detail/bessely0.hpp b/include/kyosu/types/impl/detail/bessely0.hpp deleted file mode 100644 index 17ef8ad4..00000000 --- a/include/kyosu/types/impl/detail/bessely0.hpp +++ /dev/null @@ -1,38 +0,0 @@ -//====================================================================================================================== -/* - 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 deleted file mode 100644 index d2705a28..00000000 --- a/include/kyosu/types/impl/detail/bessely1.hpp +++ /dev/null @@ -1,28 +0,0 @@ -//====================================================================================================================== -/* - 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 deleted file mode 100644 index 7b484b4b..00000000 --- a/include/kyosu/types/impl/detail/besselyn.hpp +++ /dev/null @@ -1,40 +0,0 @@ -//====================================================================================================================== -/* - 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 index 6b5b023a..c55b8f62 100644 --- a/include/kyosu/types/impl/detail/cyl_bessely1.hpp +++ b/include/kyosu/types/impl/detail/cyl_bessely1.hpp @@ -6,7 +6,7 @@ */ //====================================================================================================================== #pragma once -#include +#include namespace kyosu::_ { diff --git a/include/kyosu/types/impl/detail/cylseq.hpp b/include/kyosu/types/impl/detail/cylseq.hpp new file mode 100644 index 00000000..a9ae7bc2 --- /dev/null +++ b/include/kyosu/types/impl/detail/cylseq.hpp @@ -0,0 +1,51 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +namespace kyosu::_ +{ + + template + void dispatch(eve::tag_of, F f, size_t n, Z z, Ra & cbs) noexcept + { + std::cout << "icitte seq" << std::endl; + if constexpr(concepts::complex) + { + if constexpr(concepts::cayley_dickson) + *std::begin(cbs) = f; + else + *std::begin(cbs) = f(n, z); + if (n == 0) + return; + else + { + auto rz = kyosu::rec(z); + EVE_ASSERT(cbs.size() >= n+1, "cbs has not room enough"); + auto rs =kyosu::_:: R(n, z); + int i = n-1; + + auto cur = std::rbegin(cbs); + *cur = f(n, z); + auto pred = cur++; + while (i >= 0) + { + *cur++ = *pred++*rs; + rs = 2*(i--)*rz-kyosu::rec(rs); + } + } + } + else + { + auto p = kyosu::pure(z); + auto az = kyosu::abs(p); + cyl_bessel_cbs(n, kyosu::complex(kyosu::real(z), az), cbs); + for(size_t i = 0; i <= n; ++i) cbs[i] = kyosu::real(cbs[i]) + kyosu::ipart(cbs[i])*kyosu::sign(p); + } + } + +} diff --git a/include/kyosu/types/impl/trigo.hpp b/include/kyosu/types/impl/trigo.hpp index dd566c56..4db715fb 100644 --- a/include/kyosu/types/impl/trigo.hpp +++ b/include/kyosu/types/impl/trigo.hpp @@ -264,4 +264,21 @@ namespace kyosu::_ return cayley_extend(sinc, z); } } + + template + KYOSU_FORCEINLINE constexpr + auto dispatch(eve::tag_of const&, C const& z) noexcept + { + if constexpr(concepts::complex ) + { + auto s = kyosu::sinh(z); + using u_t = eve::underlying_type_t; + return kyosu::if_else(kyosu::abs(z) < eve::eps(eve::as(u_t())), eve::one(eve::as(u_t())), s/z); + } + else + { + return cayley_extend(sinhc, z); + } + } + } diff --git a/test/doc/cyl_bessel_yn.cpp b/test/doc/cyl_bessel_yn.cpp index bc8bf246..8ec23aca 100644 --- a/test/doc/cyl_bessel_yn.cpp +++ b/test/doc/cyl_bessel_yn.cpp @@ -30,5 +30,11 @@ 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_yn(n,z0) << "\n"; + + + for(int i=75; i < 100; ++i) + std::cout << "i " << i << " -> " << cyl_bessel_yn(i,kyosu::complex(1.0, 1.0)) << "\n"; + + return 0; } diff --git a/test/doc/sinhc.cpp b/test/doc/sinhc.cpp new file mode 100644 index 00000000..382ce095 --- /dev/null +++ b/test/doc/sinhc.cpp @@ -0,0 +1,25 @@ +#include +#include +#include + +int main() +{ + using kyosu::sinhc; + using kyosu::complex_t; + using kyosu::quaternion_t; + + std::cout << "Real: "; + std::cout << 12.9f << " -> " << sinhc(12.9f) << "\n"; + + std::cout << "Complex: "; + std::cout << kyosu::complex_t(3.5f,-2.9f) << " -> " << sinhc(kyosu::complex_t(3.5f,-2.9f)) << "\n"; + + std::cout << "Quaternion: "; + std::cout << kyosu::quaternion_t(1.,2.,3.,4.) << " -> " << sinhc(kyosu::quaternion_t(1.,2.,3.,4.)) << "\n"; + + std::cout << "SIMD: "; + using wc_t = eve::wide, eve::fixed<2>>; + std::cout << wc_t(kyosu::complex_t(1.3,-3.7)) << " -> " << sinhc(wc_t(kyosu::complex_t(1.3,-3.7))) << "\n"; + + return 0; +} diff --git a/test/doc/sph_bessel_h1_0.cpp b/test/doc/sph_bessel_h1_0.cpp new file mode 100644 index 00000000..32cdebc3 --- /dev/null +++ b/test/doc/sph_bessel_h1_0.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_h1_0; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_h1_0(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_h1_0(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_h1_0(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_h1_0(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_h1_0(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_h1_0(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_h1_0(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_h1_0(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 -> " << sph_bessel_h1_0(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_h1_0(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_h1_0(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_h1_0(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_h1_0(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-> " << sph_bessel_h1_0(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_h1_1.cpp b/test/doc/sph_bessel_h1_1.cpp new file mode 100644 index 00000000..e2aaf44c --- /dev/null +++ b/test/doc/sph_bessel_h1_1.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_h1_1; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_h1_1(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_h1_1(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_h1_1(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_h1_1(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_h1_1(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_h1_1(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_h1_1(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_h1_1(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 -> " << sph_bessel_h1_1(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_h1_1(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_h1_1(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_h1_1(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_h1_1(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-> " << sph_bessel_h1_1(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_h1n.cpp b/test/doc/sph_bessel_h1n.cpp new file mode 100644 index 00000000..dac409f4 --- /dev/null +++ b/test/doc/sph_bessel_h1n.cpp @@ -0,0 +1,41 @@ +#include +#include +#include + +int main() +{ + + size_t n = 12; + using kyosu::sph_bessel_h1n; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_h1n(n,1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_h1n(n,kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_h1n(n,15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_h1n(n,kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_h1n(n,40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_h1n(n,kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_h1n(n,60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_h1n(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 -> " << sph_bessel_h1n(n,z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_h1n(n,zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_h1n(n,1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_h1n(n,kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_h1n(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-> " << sph_bessel_h1n(n,z0) << "\n"; + + + for(int i=75; i < 100; ++i) + std::cout << "i " << i << " -> " << kyosu::sph_bessel_h1n(i,kyosu::complex(1.0, 1.0)) << "\n"; + + + return 0; +} diff --git a/test/doc/sph_bessel_h2_0.cpp b/test/doc/sph_bessel_h2_0.cpp new file mode 100644 index 00000000..1614fa4b --- /dev/null +++ b/test/doc/sph_bessel_h2_0.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_h2_0; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_h2_0(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_h2_0(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_h2_0(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_h2_0(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_h2_0(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_h2_0(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_h2_0(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_h2_0(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 -> " << sph_bessel_h2_0(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_h2_0(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_h2_0(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_h2_0(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_h2_0(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-> " << sph_bessel_h2_0(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_h2_1.cpp b/test/doc/sph_bessel_h2_1.cpp new file mode 100644 index 00000000..a5d45979 --- /dev/null +++ b/test/doc/sph_bessel_h2_1.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_h2_1; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_h2_1(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_h2_1(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_h2_1(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_h2_1(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_h2_1(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_h2_1(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_h2_1(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_h2_1(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 -> " << sph_bessel_h2_1(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_h2_1(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_h2_1(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_h2_1(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_h2_1(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-> " << sph_bessel_h2_1(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_h2n.cpp b/test/doc/sph_bessel_h2n.cpp new file mode 100644 index 00000000..fc9de767 --- /dev/null +++ b/test/doc/sph_bessel_h2n.cpp @@ -0,0 +1,40 @@ +#include +#include +#include + +int main() +{ + size_t n = 12; + using kyosu::sph_bessel_h2n; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_h2n(n,1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_h2n(n,kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_h2n(n,15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_h2n(n,kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_h2n(n,40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_h2n(n,kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_h2n(n,60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_h2n(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 -> " << sph_bessel_h2n(n,z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_h2n(n,zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_h2n(n,1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_h2n(n,kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_h2n(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-> " << sph_bessel_h2n(n,z0) << "\n"; + + + for(int i=75; i < 100; ++i) + std::cout << "i " << i << " -> " << kyosu::sph_bessel_h2n(i,kyosu::complex(1.0, 1.0)) << "\n"; + + + return 0; +} diff --git a/test/doc/sph_bessel_i1_0.cpp b/test/doc/sph_bessel_i1_0.cpp new file mode 100644 index 00000000..b3deb974 --- /dev/null +++ b/test/doc/sph_bessel_i1_0.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_i1_0; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_i1_0(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_i1_0(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_i1_0(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_i1_0(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_i1_0(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_i1_0(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_i1_0(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_i1_0(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 -> " << sph_bessel_i1_0(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_i1_0(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_i1_0(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_i1_0(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_i1_0(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-> " << sph_bessel_i1_0(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_i1_1.cpp b/test/doc/sph_bessel_i1_1.cpp new file mode 100644 index 00000000..5f14cb86 --- /dev/null +++ b/test/doc/sph_bessel_i1_1.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_i1_1; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_i1_1(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_i1_1(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_i1_1(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_i1_1(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_i1_1(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_i1_1(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_i1_1(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_i1_1(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 -> " << sph_bessel_i1_1(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_i1_1(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_i1_1(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_i1_1(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_i1_1(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-> " << sph_bessel_i1_1(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_i1n.cpp b/test/doc/sph_bessel_i1n.cpp new file mode 100644 index 00000000..cc782739 --- /dev/null +++ b/test/doc/sph_bessel_i1n.cpp @@ -0,0 +1,40 @@ +#include +#include +#include + +int main() +{ + size_t n = 12; + using kyosu::sph_bessel_i1n; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_i1n(n,1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_i1n(n,kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_i1n(n,15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_i1n(n,kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_i1n(n,40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_i1n(n,kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_i1n(n,60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_i1n(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 -> " << sph_bessel_i1n(n,z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_i1n(n,zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_i1n(n,1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_i1n(n,kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_i1n(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-> " << sph_bessel_i1n(n,z0) << "\n"; + + + for(int i=75; i < 100; ++i) + std::cout << "i " << i << " -> " << kyosu::sph_bessel_i1n(i,kyosu::complex(1.0, 1.0)) << "\n"; + + + return 0; +} diff --git a/test/doc/sph_bessel_i2_0.cpp b/test/doc/sph_bessel_i2_0.cpp new file mode 100644 index 00000000..6e227fa3 --- /dev/null +++ b/test/doc/sph_bessel_i2_0.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_i2_0; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_i2_0(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_i2_0(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_i2_0(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_i2_0(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_i2_0(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_i2_0(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_i2_0(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_i2_0(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 -> " << sph_bessel_i2_0(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_i2_0(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_i2_0(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_i2_0(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_i2_0(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-> " << sph_bessel_i2_0(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_i2_1.cpp b/test/doc/sph_bessel_i2_1.cpp new file mode 100644 index 00000000..66a4d537 --- /dev/null +++ b/test/doc/sph_bessel_i2_1.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_i2_1; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_i2_1(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_i2_1(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_i2_1(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_i2_1(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_i2_1(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_i2_1(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_i2_1(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_i2_1(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 -> " << sph_bessel_i2_1(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_i2_1(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_i2_1(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_i2_1(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_i2_1(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-> " << sph_bessel_i2_1(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_i2n.cpp b/test/doc/sph_bessel_i2n.cpp new file mode 100644 index 00000000..a31c3eae --- /dev/null +++ b/test/doc/sph_bessel_i2n.cpp @@ -0,0 +1,40 @@ +#include +#include +#include + +int main() +{ + size_t n = 12; + using kyosu::sph_bessel_i2n; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_i2n(n,1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_i2n(n,kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_i2n(n,15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_i2n(n,kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_i2n(n,40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_i2n(n,kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_i2n(n,60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_i2n(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 -> " << sph_bessel_i2n(n,z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_i2n(n,zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_i2n(n,1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_i2n(n,kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_i2n(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-> " << sph_bessel_i2n(n,z0) << "\n"; + + + for(int i=75; i < 100; ++i) + std::cout << "i " << i << " -> " << kyosu::sph_bessel_i2n(i,kyosu::complex(1.0, 1.0)) << "\n"; + + + return 0; +} diff --git a/test/doc/sph_bessel_j0.cpp b/test/doc/sph_bessel_j0.cpp new file mode 100644 index 00000000..0c065941 --- /dev/null +++ b/test/doc/sph_bessel_j0.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_j0; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_j0(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_j0(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_j0(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_j0(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_j0(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_j0(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_j0(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_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 -> " << sph_bessel_j0(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_j0(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_j0(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_j0(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_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-> " << sph_bessel_j0(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_j1.cpp b/test/doc/sph_bessel_j1.cpp new file mode 100644 index 00000000..b5edd07f --- /dev/null +++ b/test/doc/sph_bessel_j1.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_j1; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_j1(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_j1(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_j1(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_j1(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_j1(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_j1(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_j1(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_j1(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 -> " << sph_bessel_j1(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_j1(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_j1(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_j1(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_j1(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-> " << sph_bessel_j1(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_jn.cpp b/test/doc/sph_bessel_jn.cpp new file mode 100644 index 00000000..2084fdc0 --- /dev/null +++ b/test/doc/sph_bessel_jn.cpp @@ -0,0 +1,41 @@ +#include +#include +#include + +int main() +{ + + size_t n = 12; + using kyosu::sph_bessel_jn; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_jn(n,1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_jn(n,kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_jn(n,15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_jn(n,kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_jn(n,40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_jn(n,kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_jn(n,60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_jn(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 -> " << sph_bessel_jn(n,z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_jn(n,zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_jn(n,1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_jn(n,kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_jn(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-> " << sph_bessel_jn(n,z0) << "\n"; + + + for(int i=75; i < 100; ++i) + std::cout << "i " << i << " -> " << kyosu::sph_bessel_jn(i,kyosu::complex(1.0, 1.0)) << "\n"; + + + return 0; +} diff --git a/test/doc/sph_bessel_k0.cpp b/test/doc/sph_bessel_k0.cpp new file mode 100644 index 00000000..5690e586 --- /dev/null +++ b/test/doc/sph_bessel_k0.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_k0; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_k0(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_k0(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_k0(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_k0(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_k0(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_k0(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_k0(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_k0(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 -> " << sph_bessel_k0(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_k0(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_k0(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_k0(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_k0(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-> " << sph_bessel_k0(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_k1.cpp b/test/doc/sph_bessel_k1.cpp new file mode 100644 index 00000000..66c8710e --- /dev/null +++ b/test/doc/sph_bessel_k1.cpp @@ -0,0 +1,34 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_k1; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_k1(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_k1(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_k1(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_k1(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_k1(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_k1(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_k1(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_k1(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 -> " << sph_bessel_k1(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_k1(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_k1(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_k1(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0j " << " -> " << sph_bessel_k1(kyosu::quaternion(1.0, 0.0, 36.0, 0.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_k1(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-> " << sph_bessel_k1(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_kn.cpp b/test/doc/sph_bessel_kn.cpp new file mode 100644 index 00000000..9643e45a --- /dev/null +++ b/test/doc/sph_bessel_kn.cpp @@ -0,0 +1,41 @@ +#include +#include +#include + +int main() +{ + + size_t n = 12; + using kyosu::sph_bessel_kn; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_kn(n,1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_kn(n,kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_kn(n,15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_kn(n,kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_kn(n,40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_kn(n,kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_kn(n,60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_bessel_kn(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 -> " << sph_bessel_kn(n,z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_kn(n,zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_kn(n,1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_kn(n,kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_bessel_kn(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-> " << sph_bessel_kn(n,z0) << "\n"; + + + for(int i=75; i < 100; ++i) + std::cout << "i " << i << " -> " << kyosu::sph_bessel_kn(i,kyosu::complex(1.0, 1.0)) << "\n"; + + + return 0; +} diff --git a/test/doc/sph_bessel_y0.cpp b/test/doc/sph_bessel_y0.cpp new file mode 100644 index 00000000..a95f4c7b --- /dev/null +++ b/test/doc/sph_bessel_y0.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_y0; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_y0(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_y0(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_y0(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_y0(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_y0(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_y0(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_y0(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_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 -> " << sph_bessel_y0(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_y0(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_y0(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_y0(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_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-> " << sph_bessel_y0(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_y1.cpp b/test/doc/sph_bessel_y1.cpp new file mode 100644 index 00000000..e1428175 --- /dev/null +++ b/test/doc/sph_bessel_y1.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::sph_bessel_y1; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << sph_bessel_y1(1.0) << "\n"; + std::cout << "1+0i " << " -> " << sph_bessel_y1(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << sph_bessel_y1(15.0) << "\n"; + std::cout << "15+0i " << " -> " << sph_bessel_y1(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << sph_bessel_y1(40.0) << "\n"; + std::cout << "40+0i " << " -> " << sph_bessel_y1(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << sph_bessel_y1(60.0) << "\n"; + std::cout << "60+0i " << " -> " << sph_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 -> " << sph_bessel_y1(z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(sph_bessel_y1(zz)) << "\n"; + + std::cout << "1.0 " << " -> " << sph_bessel_y1(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << sph_bessel_y1(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << sph_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-> " << sph_bessel_y1(z0) << "\n"; + return 0; +} diff --git a/test/doc/sph_bessel_yn.cpp b/test/doc/sph_bessel_yn.cpp new file mode 100644 index 00000000..46b23843 --- /dev/null +++ b/test/doc/sph_bessel_yn.cpp @@ -0,0 +1,40 @@ +#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"; + + + for(int i=75; i < 100; ++i) + std::cout << "i " << i << " -> " << kyosu::sph_bessel_yn(i,kyosu::complex(1.0, 1.0)) << "\n"; + + + return 0; +} diff --git a/test/unit/function/sinhc.cpp b/test/unit/function/sinhc.cpp new file mode 100644 index 00000000..93db6f84 --- /dev/null +++ b/test/unit/function/sinhc.cpp @@ -0,0 +1,36 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +# if __has_include () +# include +# define HAS_BOOST +# endif + +TTS_CASE_WITH ( "Check kyosu::sinhc over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10)) + ) +(auto data) +{ + TTS_ULP_EQUAL(kyosu::sinhc(data), eve::sinhc(data), 0.5); +}; + +TTS_CASE_WITH ( "Check kyosu::sinhc over quaternion" + , kyosu::simd_real_types + , tts::generate ( tts::randoms(-10,10), tts::randoms(-10,10) + , tts::randoms(-10,10), tts::randoms(-10,10) + ) + ) +(T r, T i, T j, T k) +{ + using ke_t = kyosu::quaternion_t; + auto q = ke_t(r,i, j, k); + TTS_RELATIVE_EQUAL(kyosu::sinhc(q), kyosu::sinh(q)/q, 1e-4); +}; diff --git a/test/unit/function/sph_bessel_h0.cpp b/test/unit/function/sph_bessel_h0.cpp new file mode 100644 index 00000000..f90f52cd --- /dev/null +++ b/test/unit/function/sph_bessel_h0.cpp @@ -0,0 +1,46 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_h1_0 over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10)) + ) +(T a0, T a1) +{ + using u_t = eve::underlying_type_t; + auto i = kyosu::complex(u_t(0), u_t(1)); + auto h1_0 = [i](auto z){ return kyosu::sph_bessel_j0(z)+i*kyosu::sph_bessel_y0(z) ; }; + auto h2_0 = [i](auto z){ return kyosu::sph_bessel_j0(z)-i*kyosu::sph_bessel_y0(z) ; }; + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(a0); + auto im = i*a1; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h1_0(re), h1_0(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h1_0(im), h1_0(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h1_0(z) , h1_0(z) , 1.0e-4) << i << " <- " << z << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h2_0(re), h2_0(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h2_0(im), h2_0(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h2_0(z) , h2_0(z) , 1.0e-4) << i << " <- " << z << '\n'; +}; + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_h1_0 over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T , T ) +{ + using kyosu::sph_bessel_h1_0; + using kyosu::sph_bessel_h1_1; + auto z = kyosu::complex(T(0), T(0)); + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_h1_1(z)))); + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_h1_0(z)))); +}; diff --git a/test/unit/function/sph_bessel_h1.cpp b/test/unit/function/sph_bessel_h1.cpp new file mode 100644 index 00000000..efe5b20d --- /dev/null +++ b/test/unit/function/sph_bessel_h1.cpp @@ -0,0 +1,46 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_h1_1 over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10)) + ) +(T a0, T a1) +{ + using u_t = eve::underlying_type_t; + auto i = kyosu::complex(u_t(0), u_t(1)); + auto h1_1 = [i](auto z){ return kyosu::sph_bessel_j1(z)+i*kyosu::sph_bessel_y1(z) ; }; + auto h2_1 = [i](auto z){ return kyosu::sph_bessel_j1(z)-i*kyosu::sph_bessel_y1(z) ; }; + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(a0); + auto im = i*a1; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h1_1(re), h1_1(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h1_1(im), h1_1(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h1_1(z) , h1_1(z) , 1.0e-4) << i << " <- " << z << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h2_1(re), h2_1(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h2_1(im), h2_1(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h2_1(z) , h2_1(z) , 1.0e-4) << i << " <- " << z << '\n'; +}; + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_h1_1 over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T , T ) +{ + using kyosu::sph_bessel_h1_1; + using kyosu::sph_bessel_h2_1; + auto z = kyosu::complex(T(0), T(0)); + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_h1_1(z)))); + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_h2_1(z)))); +}; diff --git a/test/unit/function/sph_bessel_hn.cpp b/test/unit/function/sph_bessel_hn.cpp new file mode 100644 index 00000000..0eae5cee --- /dev/null +++ b/test/unit/function/sph_bessel_hn.cpp @@ -0,0 +1,77 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_h1nn over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ +// using u_t = eve::underlying_type_t; + if constexpr(sizeof(T) == 8) + { + { + std::array re{-3.0741715027083147e+00, 4.3426648638416063e+00, -2.7981941356021967e+00, -8.8600087717455200e-01, -2.1297833758572637e+00, -4.9849851779139147e+00, 1.5223768684011907e+00, -8.4260002884283591e-01}; + + std::array im{-2.0307188187984737e+00, -4.0705193405198852e+00, 2.3547162061032791e+00, -2.5738011925258233e+00, 2.6849682520111706e+00, 2.5683492208093019e+00, 3.7050613060981217e+00, -1.9582614234846174e+00}; + + std::array reres{-1.8987639369507804e+01, -9.7846582033372931e-02, 7.6790116788295499e+00, -2.0361029682674200e+02, -2.7927440764619799e+01, 1.3454938584765222e-01, 3.9201689873947365e+00, -1.6300834687965348e+03}; + + std::array imres{2.1905468729751316e+00, -8.1937022872644705e-02, 1.6738406475084535e+01, 2.0527437676205267e+00, 3.5791027983986723e+00, -6.3121489926765517e-01, 3.4600793328396318e+00, -1.1728487363689339e+03}; + for(int i=0; i < 8; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h1n(8, c), res, 1.0e-4) << i << " <- " << c << '\n'; + } + } + { + std::array re{4.8023241784273516e+00, -7.3624994664291843e-01, 2.7845895150886490e+00, -4.2439315522389700e+00, -3.2700626063635507e+00, 3.2777149501976064e+00, -1.7350443377902092e+00, 3.4240826382686604e+00}; + std::array im{-1.0943783015244868e+00, 4.2508564240326341e-01, -1.8583059717139760e+00, -4.8549618376057779e+00, 2.1934429650892029e+00, 2.7593426313960578e+00, -2.0573580642750766e+00, 8.2956711439982622e-01}; + std::array reres{ -2.537018910232239 ,-8841141.702494428,42.74928523513201,-0.006411322923734045,-10.58900986724366,-2.513588322804079,-252.9552712623419,32.93818432005056}; + std::array imres{-0.1086979192666555,+183900.2111111085,+10.09288763529341,-0.07774596762539486,-1.210214973107798,+3.719306160697736,-66.91213603662391,-12.65422016639942}; + for(int i=0; i < 8; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_h2n(8, c), res, 1.0e-4); + } + } + } +}; + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_jn over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + if constexpr(sizeof(eve::underlying_type_t) == 8) + { + using u_t = eve::underlying_type_t; + auto c = kyosu::complex(a0, a1); + auto cm = -c; + auto zer = kyosu::complex(T(0), T(0)); + using kyosu::sph_bessel_h1n; + using kyosu::sph_bessel_h2n; + + for(int i=1; i < 4 ; ++i) + { + auto h1nc = sph_bessel_h1n(i, c); + auto h2nc = sph_bessel_h2n(i, c); + TTS_IEEE_EQUAL(h1nc, eve::sign_alternate(u_t(i))*sph_bessel_h2n(i, cm)) << "c " << c << "\n"; + TTS_IEEE_EQUAL(h2nc, eve::sign_alternate(u_t(i))*sph_bessel_h1n(i, cm)) << "c " << c << "\n"; + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_h1n(i, zer)))) << i << sph_bessel_h1n(i, zer) << '\n'; + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_h2n(i, zer)))) << i << sph_bessel_h2n(i, zer) << '\n'; + } + } +}; diff --git a/test/unit/function/sph_bessel_i0.cpp b/test/unit/function/sph_bessel_i0.cpp new file mode 100644 index 00000000..81b8a252 --- /dev/null +++ b/test/unit/function/sph_bessel_i0.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_i1_0 over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10)) + ) +(T a0, T a1) +{ + auto i1_0 = [](auto z){ return kyosu::sinh(z)/z; }; + auto i2_0 = [](auto z){ return kyosu::cosh(z)/z; }; + using u_t = eve::underlying_type_t; + auto i = kyosu::complex(u_t(0), u_t(1)); + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(a0); + auto im = i*a1; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1_0(re), i1_0(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1_0(im), i1_0(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1_0(z) , i1_0(z) , 1.0e-4) << i << " <- " << z << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2_0(re), i2_0(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2_0(im), i2_0(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2_0(z) , i2_0(z) , 1.0e-4) << i << " <- " << z << '\n'; +}; + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_i1_0 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 cm= -c; + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + using kyosu::sph_bessel_i1_0; + using kyosu::sph_bessel_i1_1; + auto i1_0c = sph_bessel_i1_0(c); + auto i1_1c = sph_bessel_i1_1(c); + TTS_EQUAL(i1_0c, sph_bessel_i1_0(cm)); + TTS_EQUAL(i1_1c, -sph_bessel_i1_1(cm)); + TTS_EQUAL(i1_0c, kyosu::conj(sph_bessel_i1_0(cb))); + TTS_EQUAL(i1_1c, kyosu::conj(sph_bessel_i1_1(cb))); + TTS_EXPECT(eve::all(kyosu::is_real(sph_bessel_i1_0(cr)))); + TTS_EXPECT(eve::all(kyosu::is_real(sph_bessel_i1_1(cr)))); + TTS_EXPECT(eve::all(kyosu::is_real(sph_bessel_i1_0(ci)))); + TTS_EXPECT(eve::all(kyosu::is_pure(sph_bessel_i1_1(ci)))); + auto z = kyosu::complex(T(0), T(0)); + auto o = kyosu::complex(T(1), T(0)); + TTS_EQUAL(sph_bessel_i1_0(z), o); + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_i1_1(z)))); +}; diff --git a/test/unit/function/sph_bessel_i1.cpp b/test/unit/function/sph_bessel_i1.cpp new file mode 100644 index 00000000..a5e8384c --- /dev/null +++ b/test/unit/function/sph_bessel_i1.cpp @@ -0,0 +1,70 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i1 over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10)) + ) +(T a0, T a1) +{ + auto i1_1 = [](auto z){ + auto [sh, ch] = kyosu::sinhcosh(z); + auto rz = kyosu::rec(z); + return kyosu::fnma(sh, rz, ch)*rz; + }; + auto i2_1 = [](auto z){ + auto [sh, ch] = kyosu::sinhcosh(z); + auto rz = kyosu::rec(z); + return kyosu::fnma(ch, rz, sh)*rz; + }; + using u_t = eve::underlying_type_t; + auto i = kyosu::complex(u_t(0), u_t(1)); + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(a0); + auto im = i*a1; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1_1(re), i1_1(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1_1(im), i1_1(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1_1(z) , i1_1(z) , 1.0e-4) << i << " <- " << z << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2_1(re), i2_1(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2_1(im), i2_1(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2_1(z) , i2_1(z) , 1.0e-4) << i << " <- " << z << '\n'; +}; + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_i1 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 cm= -c; + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + using kyosu::sph_bessel_i1_1; + using kyosu::sph_bessel_i2_1; + auto i1c = sph_bessel_i1_1(c); + auto i2c = sph_bessel_i2_1(c); + TTS_EQUAL(i1c, -sph_bessel_i1_1(cm)); + TTS_EQUAL(i2c, sph_bessel_i2_1(cm)); + TTS_EQUAL(i1c, kyosu::conj(sph_bessel_i1_1(cb))); + TTS_EQUAL(i2c, kyosu::conj(sph_bessel_i2_1(cb))); + TTS_EXPECT(eve::all(kyosu::is_real(sph_bessel_i1_1(cr)))); + TTS_EXPECT(eve::all(kyosu::is_real(sph_bessel_i2_1(cr)))); + TTS_EXPECT(eve::all(kyosu::is_pure(sph_bessel_i1_1(ci)))); + TTS_EXPECT(eve::all(kyosu::is_real(sph_bessel_i2_1(ci)))); + auto z = kyosu::complex(T(0), T(0)); + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_i1_1(z)))); + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_i2_1(z)))); + +}; diff --git a/test/unit/function/sph_bessel_in.cpp b/test/unit/function/sph_bessel_in.cpp new file mode 100644 index 00000000..0b50d34d --- /dev/null +++ b/test/unit/function/sph_bessel_in.cpp @@ -0,0 +1,76 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_i1nn over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + if constexpr(sizeof(T) == 8) + { + { + std::array re{-1.3102964777073500e+00, 2.6429797151271517e+00, -1.2991425964097203e+00, -2.0616201216776595e+00, -3.5692327918927735e-01, 4.5409169416450617e+00, 1.3915394762511490e+00, -1.5742664543118334e+00}; + std::array im{-4.9978344205200456e+00, 4.6966388867930240e+00, 4.5575687436111947e+00, -4.3856054937978231e+00, 1.7496198325295331e+00, -4.8976206253867760e+00, -2.2292456733273172e+00, -3.6547854593577123e+00}; + std::array reres{-8.9448633564799367e-04, -1.3550451109077518e-02, -1.4079674805043841e-03, -5.9344800455695375e-03, -1.6670926554323534e-08, 1.0879026036746578e-02, -2.4460612590389117e-05, -1.3445907340862448e-03}; + std::array imres{-7.8433567383746679e-03, 3.9123036962636004e-03, 4.2071250786199152e-03, -7.1820574884934665e-04, 2.7760904475739103e-06, -1.1171678180080313e-01, -5.5883425389321553e-05, -2.7135427306845860e-04}; + for(int i=0; i < 8; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i1n(8, c), res, 1.0e-7) << i << " <- " << c << '\n'; + } + } + { + std::array re{-1.0353298912411590e+00, -4.0590492982139348e-01, -3.7533860841866939e+00, 2.5026965372926147e-01, 1.3027543743434511e-01,-3.9640689770100854e+00, -1.0995476278801330e+00, -3.8763377787749520e-01}; + std::array im{-4.2831047668722677e+00, 7.7061553129385785e-01, 1.3657001488056097e+00, 4.1210952543771051e-01, 4.4818969533279871e+00, -3.7525013370467541e+00, 4.8289663184328813e-01, -7.0094140409912353e-01}; + std::array reres{-5.7006149471406644e+00, 6.6457154244604167e+06, 4.9722413633652538e+00, -1.4179514732481084e+09, 1.2132650006681511e+00, -4.7090187937440880e-02, 3.0782381172903924e+05, 1.4851214398026228e+07}; + std::array imres{-1.3655417114292883e+00, 2.5782447371366373e+06, 1.6560730551132823e+00, -2.7650888434981310e+08, -5.5042068213386841e+00, 4.0322394841215375e-01, 2.1882447301295877e+05, -2.7656213800463108e+06}; + for(int i=0; i < 8; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_i2n(8, c), res, 1.0e-7) << i << " <- " << c << '\n'; + } + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_jn over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + if constexpr(sizeof(eve::underlying_type_t) == 8) + { + using u_t = eve::underlying_type_t; + auto c = kyosu::complex(a0, a1); + auto cm = -c; + auto zer = kyosu::complex(T(0), T(0)); + // auto one = kyosu::complex(T(1), T(0)); + + using kyosu::sph_bessel_i1n; + using kyosu::sph_bessel_i2n; + + for(int i=1; i < 4 ; ++i) + { + auto i1nc = sph_bessel_i1n(i, c); + auto i2nc = sph_bessel_i2n(i, c); + TTS_IEEE_EQUAL(i1nc, eve::sign_alternate(u_t(i))*sph_bessel_i1n(i, cm)) << "c " << c << "\n"; + TTS_IEEE_EQUAL(i2nc, eve::sign_alternate(u_t(i+1))*sph_bessel_i2n(i, cm)) << "c " << c << "\n"; + TTS_EXPECT(eve::all(kyosu::is_eqz(sph_bessel_i1n(i, zer)))) << i << sph_bessel_i1n(i, zer) << '\n'; + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_i2n(i, zer)))) << i << sph_bessel_i2n(i, zer) << '\n'; + } + } +}; diff --git a/test/unit/function/sph_bessel_j0.cpp b/test/unit/function/sph_bessel_j0.cpp new file mode 100644 index 00000000..cd8fa9cf --- /dev/null +++ b/test/unit/function/sph_bessel_j0.cpp @@ -0,0 +1,50 @@ +//====================================================================================================================== +/* + 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), + tts::randoms(-10,10)) + ) +(T a0, T a1) +{ + using u_t = eve::underlying_type_t; + auto i = kyosu::complex(u_t(0), u_t(1)); + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(a0); + auto im = i*a1; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_j0(re), kyosu::sinc(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_j0(im), kyosu::sinc(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_j0(z) , kyosu::sinc(z) , 1.0e-4) << i << " <- " << z << '\n'; +}; + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_j0 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 cm= -c; + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + using kyosu::sph_bessel_j0; + auto j0c = sph_bessel_j0(c); + TTS_EQUAL(j0c, sph_bessel_j0(cm)); + TTS_EQUAL(j0c, kyosu::conj(sph_bessel_j0(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 o = kyosu::complex(T(1), T(0)); + TTS_IEEE_EQUAL(sph_bessel_j0(z), o); +}; diff --git a/test/unit/function/sph_bessel_j1.cpp b/test/unit/function/sph_bessel_j1.cpp new file mode 100644 index 00000000..fcdd8de1 --- /dev/null +++ b/test/unit/function/sph_bessel_j1.cpp @@ -0,0 +1,50 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j1 over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10)) + ) +(T a0, T a1) +{ + auto j1 = [](auto z){ return (kyosu::sinc(z)-kyosu::cos(z))/z; }; + using u_t = eve::underlying_type_t; + auto i = kyosu::complex(u_t(0), u_t(1)); + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(a0); + auto im = i*a1; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_j1(re), j1(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_j1(im), j1(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_j1(z) , j1(z) , 1.0e-4) << i << " <- " << z << '\n'; +}; + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_j1 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 cm= -c; + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + using kyosu::sph_bessel_j1; + auto j1c = sph_bessel_j1(c); + TTS_EQUAL(j1c, -sph_bessel_j1(cm)); + TTS_EQUAL(j1c, kyosu::conj(sph_bessel_j1(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(sph_bessel_j1(z), z); +}; diff --git a/test/unit/function/sph_bessel_jn.cpp b/test/unit/function/sph_bessel_jn.cpp new file mode 100644 index 00000000..97707021 --- /dev/null +++ b/test/unit/function/sph_bessel_jn.cpp @@ -0,0 +1,67 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_jn over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + if constexpr(sizeof(T) == 8) + { + std::array re{-3.5864413497599945e-01, 1.8305506745510136e+00, 4.8145177703097399e+00, 4.7301901125586436e+00, -4.3947566261794710e+00, 1.0212518380533187e+00, -4.1514744061148807e+00, 1.2766581423615309e-01}; + + std::array im{-1.4900727734412866e+00, 4.3733972760463349e+00, 6.3220189518343628e-01, 5.1967694618977767e-01, 6.2761895523580580e-01, 3.8772322566505766e+00, -4.6074082935612806e+00, -4.4849065338089478e+00}; + + std::array reres{2.6695000315635153e-02, -2.2215759788276843e+00, 2.5805202705605512e-01, 2.5337598654371535e-01, -2.6051834312740757e-01, -1.1816979617769454e+00, -3.2042363997219185e+00, -3.3242357291991742e-01}; + + std::array imres{2.7709520397028151e-02, 9.1259760936216550e-01, -1.9344122386863866e-02, -1.1252144284238647e-02, 8.5011995341820815e-03, -5.1276579958119139e-01, -2.1027562948226137e+00, 2.3711156404785037e+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::sph_bessel_jn(3, c), res, 1.0e-7) << i << " <- " << c << '\n'; + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_jn over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + if constexpr(sizeof(eve::underlying_type_t) == 8) + { + using u_t = eve::underlying_type_t; + auto c = kyosu::complex(a0, a1); + auto 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::sph_bessel_jn; + + for(int i=0; i < 10 ; ++i) + { + auto jnc = sph_bessel_jn(i, c); + TTS_IEEE_EQUAL(jnc, eve::sign_alternate(u_t(i))*sph_bessel_jn(i, cm)) << "c " << c << "\n"; + TTS_IEEE_EQUAL(jnc, kyosu::conj(sph_bessel_jn(i, cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + TTS_IEEE_EQUAL(sph_bessel_jn(i, zer), i ? zer : one) << i << '\n'; + } + } +}; diff --git a/test/unit/function/sph_bessel_k0.cpp b/test/unit/function/sph_bessel_k0.cpp new file mode 100644 index 00000000..5c890df3 --- /dev/null +++ b/test/unit/function/sph_bessel_k0.cpp @@ -0,0 +1,48 @@ +//====================================================================================================================== +/* + 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), + tts::randoms(-10,10)) + ) +(T a0, T a1) +{ + using u_t = eve::underlying_type_t; + auto i = kyosu::complex(u_t(0), u_t(1)); + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(a0); + auto im = i*a1; + auto k0 = [](auto z){ return eve::pio_2(eve::as())*kyosu::exp(-z)/z; }; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_k0(re), k0(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_k0(im), k0(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_k0(z) , k0(z) , 1.0e-4) << i << " <- " << z << '\n'; +}; + +TTS_CASE_WITH ( "Check kyosu::sph_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); + auto ci= kyosu::complex(T(0), a1); + using kyosu::sph_bessel_k0; + auto k0c = sph_bessel_k0(c); + TTS_EQUAL(k0c, kyosu::conj(sph_bessel_k0(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_EXPECT(eve::all(kyosu::is_nan(sph_bessel_k0(z)))); +}; diff --git a/test/unit/function/sph_bessel_k1.cpp b/test/unit/function/sph_bessel_k1.cpp new file mode 100644 index 00000000..b0d1c560 --- /dev/null +++ b/test/unit/function/sph_bessel_k1.cpp @@ -0,0 +1,51 @@ +//====================================================================================================================== +/* + 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), + tts::randoms(-10,10)) + ) +(T a0, T a1) +{ + using u_t = eve::underlying_type_t; + auto i = kyosu::complex(u_t(0), u_t(1)); + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(a0); + auto im = i*a1; + auto k1 = [](auto z){ + auto rz = kyosu::rec(z); + return eve::pio_2(eve::as())*kyosu::exp(-z)*kyosu::fma(rz, rz, rz); + }; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_k1(re), k1(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_k1(im), k1(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_k1(z) , k1(z) , 1.0e-4) << i << " <- " << z << '\n'; +}; + +TTS_CASE_WITH ( "Check kyosu::sph_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); + auto ci= kyosu::complex(T(0), a1); + using kyosu::sph_bessel_k1; + auto k1c = sph_bessel_k1(c); + TTS_EQUAL(k1c, kyosu::conj(sph_bessel_k1(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_EXPECT(eve::all(kyosu::is_nan(sph_bessel_k1(z)))); +}; diff --git a/test/unit/function/sph_bessel_kn.cpp b/test/unit/function/sph_bessel_kn.cpp new file mode 100644 index 00000000..ba556045 --- /dev/null +++ b/test/unit/function/sph_bessel_kn.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::sph_bessel_kn over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + if constexpr(sizeof(T) == 8) + { + std::array re{2.0344637669563861e+00, -4.3813469489459145e+00, -2.6657638099374692e+00, 3.3889708858072520e+00, -3.7647998208558220e+00, 1.7491848302215496e+00, 6.6959832811814834e-02, -1.0359975448622161e+00}; + + std::array im{9.1599183409523333e-01, 3.8563462289706973e-01, 3.1031208865787474e-01, -3.1600822008299359e+00, 4.0144397073512161e+00, -3.6311009559927134e+00, -2.8259820882536459e+00, 4.1712992187897164e+00}; + + std::array reres{-1.4678177240255734e+03, -1.7772703403589662e+00, -1.6005659389435411e+02, 1.1469395680075878e+00, 2.9114121290916373e-01, -6.8942371077974842e+00, 7.2045172227253374e+01, -1.0631979003342050e+01}; + + std::array imres{1.4756336825094897e+03, -2.1777133097986439e+00, -3.1159733154561093e+02, 2.7355335304581043e+00, -5.3062640683955642e-01, -1.4181875854319991e+01, 3.5537641050392568e+02, 3.3112069122251264e+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::sph_bessel_kn(8, c), res, 1.0e-7) << i << " <- " << c << '\n'; + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_kn over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + if constexpr(sizeof(eve::underlying_type_t) == 8) + { + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + auto zer = kyosu::complex(T(0), T(0)); + + using kyosu::sph_bessel_kn; + + for(int i=2; i < 10 ; ++i) + { + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_kn(i, zer)))) << i << '\n'; + } + } +}; diff --git a/test/unit/function/sph_bessel_y0.cpp b/test/unit/function/sph_bessel_y0.cpp new file mode 100644 index 00000000..b2250910 --- /dev/null +++ b/test/unit/function/sph_bessel_y0.cpp @@ -0,0 +1,50 @@ +//====================================================================================================================== +/* + 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), + tts::randoms(-10,10)) + ) +(T a0, T a1) +{ + using u_t = eve::underlying_type_t; + auto i = kyosu::complex(u_t(0), u_t(1)); + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(a0); + auto im = i*a1; + auto y0 = [](auto z){return -kyosu::cos(z)/z; }; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_y0(re), y0(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_y0(im), y0(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_y0(z) , y0(z) , 1.0e-4) << i << " <- " << z << '\n'; +}; + +TTS_CASE_WITH ( "Check kyosu::sph_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 cm= -c; + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + using kyosu::sph_bessel_y0; + auto y0c = sph_bessel_y0(c); + TTS_EQUAL(y0c, -sph_bessel_y0(cm)); + TTS_EQUAL(y0c, kyosu::conj(sph_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)); + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_y0(z)))); +}; diff --git a/test/unit/function/sph_bessel_y1.cpp b/test/unit/function/sph_bessel_y1.cpp new file mode 100644 index 00000000..1306f57f --- /dev/null +++ b/test/unit/function/sph_bessel_y1.cpp @@ -0,0 +1,50 @@ +//====================================================================================================================== +/* + 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), + tts::randoms(-10,10)) + ) +(T a0, T a1) +{ + using u_t = eve::underlying_type_t; + auto i = kyosu::complex(u_t(0), u_t(1)); + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(a0); + auto im = i*a1; + auto y1 = [](auto z){return -kyosu::cos(z)/kyosu::sqr(z)-kyosu::sinc(z); }; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_y1(re), y1(re), 1.0e-4) << i << " <- " << re << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_y1(im), y1(im), 1.0e-4) << i << " <- " << im << '\n'; + TTS_RELATIVE_EQUAL(kyosu::sph_bessel_y1(z) , y1(z) , 1.0e-4) << i << " <- " << z << '\n'; +}; + +TTS_CASE_WITH ( "Check kyosu::sph_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 cm= -c; + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + using kyosu::sph_bessel_y1; + auto y1c = sph_bessel_y1(c); + TTS_EQUAL(y1c, sph_bessel_y1(cm)); + TTS_EQUAL(y1c, kyosu::conj(sph_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)); + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_y1(z)))); +}; diff --git a/test/unit/function/sph_bessel_yn.cpp b/test/unit/function/sph_bessel_yn.cpp new file mode 100644 index 00000000..00c4f73e --- /dev/null +++ b/test/unit/function/sph_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::sph_bessel_yn over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + if constexpr(sizeof(T) == 8) + { + std::array re{2.9751292610110780e+00, 1.3139296459466121e+00, 3.2268716431231459e+00, 1.6502599567093656e+00, 4.6332810203103314e+00, 2.8212863760434983e+00, -1.1810230326333682e+00, 4.4699044019977441e+00}; + + std::array im{-2.0222848019757578e+00, 4.6155627764756613e+00, 2.1866177359975039e+00, -1.5599498003423562e+00, 3.2640938797973753e+00, 1.4714931919159735e+00, -3.8260183135371872e+00, -2.1131147779384407e-01}; + + std::array reres{-3.5158398033526106e-01, 4.9066840081936158e-01, -3.9771160672274486e-01, 3.4727405740898953e-01, -1.0049148303947214e-01, -2.9662326625355484e-01, 2.8364317746458689e-01, -1.1622279259142240e-01}; + + std::array imres{-2.5850292463184665e-01, -2.8908255599320762e+00, 3.1145881726375613e-01, -2.6570747897125668e-01, 1.2559393288754663e+00, 3.0517784454448199e-01, -1.2282698054615229e+00, -4.3334436047326737e-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::sph_bessel_yn(3, c), res, 1.0e-7) << i << " <- " << c << '\n'; + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::sph_bessel_yn over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + if constexpr(sizeof(eve::underlying_type_t) == 8) + { + using u_t = eve::underlying_type_t; + auto c = kyosu::complex(a0, a1); + auto cm= -c; + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + auto zer = kyosu::complex(T(0), T(0)); + + using kyosu::sph_bessel_yn; + + for(int i=2; i < 10 ; ++i) + { + auto ync = sph_bessel_yn(i, c); + TTS_IEEE_EQUAL(ync, eve::sign_alternate(u_t(i+1))*sph_bessel_yn(i, cm)) << "i " << i << " c " << c << "\n"; + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + TTS_EXPECT(eve::all(kyosu::is_nan(sph_bessel_yn(i, zer)))) << i << '\n'; + } + } +};