From 98e4b3b955f8dc915be13ffc710d87073b78d859 Mon Sep 17 00:00:00 2001 From: jtlap Date: Thu, 26 Oct 2023 16:55:39 +0200 Subject: [PATCH 01/17] lj0 lj1 ljn still incorrect --- include/kyosu/functions.hpp | 4 +- include/kyosu/functions/cyl_lbessel_j0.hpp | 80 +++++++ include/kyosu/functions/cyl_lbessel_j1.hpp | 80 +++++++ include/kyosu/functions/cyl_lbessel_jn.hpp | 81 +++++++ include/kyosu/types/impl/bessel.hpp | 4 + .../kyosu/types/impl/detail/bessel_utils.hpp | 226 ++++++++++++++++++ include/kyosu/types/impl/detail/lj0.hpp | 29 +++ include/kyosu/types/impl/detail/lj1.hpp | 29 +++ include/kyosu/types/impl/detail/ljn.hpp | 51 ++++ test/doc/cyl_lbessel_jn.cpp | 27 +++ 10 files changed, 610 insertions(+), 1 deletion(-) create mode 100644 include/kyosu/functions/cyl_lbessel_j0.hpp create mode 100644 include/kyosu/functions/cyl_lbessel_j1.hpp create mode 100644 include/kyosu/functions/cyl_lbessel_jn.hpp create mode 100644 include/kyosu/types/impl/detail/bessel_utils.hpp create mode 100644 include/kyosu/types/impl/detail/lj0.hpp create mode 100644 include/kyosu/types/impl/detail/lj1.hpp create mode 100644 include/kyosu/types/impl/detail/ljn.hpp create mode 100644 test/doc/cyl_lbessel_jn.cpp diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index 70bdc83c..537526b3 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -53,7 +53,9 @@ #include #include #include -#include +#include +#include +#include #include #include #include diff --git a/include/kyosu/functions/cyl_lbessel_j0.hpp b/include/kyosu/functions/cyl_lbessel_j0.hpp new file mode 100644 index 00000000..24ed12c9 --- /dev/null +++ b/include/kyosu/functions/cyl_lbessel_j0.hpp @@ -0,0 +1,80 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_lbessel_j0: eve::elementwise + { + using callable_tag_type = callable_cyl_lbessel_j0; + + KYOSU_DEFERS_CALLABLE(cyl_lbessel_j0_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept + { + + return log(cyl_bessel_j0(complex(v))); + } + + template + KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) + { + return eve::tag_invoke(*this, target); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_lbessel_j0 +//! @brief Computes the logarithm of the bessel function of the first kind, +//! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) +//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. +//! +//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_lbessel_j0(T z) noexcept; +//! template constexpr T cyl_lbessel_j0(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return \f$\log(J_0(z))\f$. If T is a floating point value (real), a complex value is returned. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_lbessel_j0.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_lbessel_j0 cyl_lbessel_j0 = {}; +} diff --git a/include/kyosu/functions/cyl_lbessel_j1.hpp b/include/kyosu/functions/cyl_lbessel_j1.hpp new file mode 100644 index 00000000..88cc8569 --- /dev/null +++ b/include/kyosu/functions/cyl_lbessel_j1.hpp @@ -0,0 +1,80 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_lbessel_j1: eve::elementwise + { + using callable_tag_type = callable_cyl_lbessel_j1; + + KYOSU_DEFERS_CALLABLE(cyl_lbessel_j1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept + { + callable_cyl_lbessel_j1 fn{}; + return fn(complex(v)); + } + + template + KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) + { + return eve::tag_invoke(*this, target); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_lbessel_j1 +//! @brief Computes the logarithm of the Bessel function of the first kind, +//! \f$ J_1(x)=\frac1{\pi }\int _{0}^{\pi}\cos(\tau-x\sin \tau )\,\mathrm {d} \tau \f$ +//! extended to the complex plane and cayley_dickson values. +//! +//! \f$ J_1(\f$ is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_lbessel_j1(T z) noexcept; +//! template constexpr T cyl_lbessel_j1(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return the cylindrical \f$J_0(z)\f$. If T is a floating point value (real), a complex value is returned. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_lbessel_j1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_lbessel_j1 cyl_lbessel_j1 = {}; +} diff --git a/include/kyosu/functions/cyl_lbessel_jn.hpp b/include/kyosu/functions/cyl_lbessel_jn.hpp new file mode 100644 index 00000000..5cd51c65 --- /dev/null +++ b/include/kyosu/functions/cyl_lbessel_jn.hpp @@ -0,0 +1,81 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_lbessel_jn: eve::elementwise + { + using callable_tag_type = callable_cyl_lbessel_jn; + + KYOSU_DEFERS_CALLABLE(cyl_lbessel_jn_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& v) noexcept + { + callable_cyl_lbessel_jn fn{}; + return 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 +//! @{ +//! @brief Computes the logarithm of the Bessel functions of the first kind, +//! \f$ J_{n}(x)=\sum_{p=0}^{\infty}{\frac{(-1)^p}{p!\,\Gamma (p+n +1)}} +//! {\left({x \over 2}\right)}^{2p+n }\f$. +//! +//! \f$ J_{n}\f$ is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which +//! \f$ y(0) = 0\f$ if \f$n \ne 0\f$ else \f$1\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_lbessel_jn(int n, T z) noexcept; +//! template constexpr T cyl_lbessel_jn(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return the cylindrical \f$J_n(z)\f$. If T is a floating point value (real), a complex value is returned. +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_lbessel_jn.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_lbessel_jn cyl_lbessel_jn = {}; +} diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp index 71257611..b443efd4 100644 --- a/include/kyosu/types/impl/bessel.hpp +++ b/include/kyosu/types/impl/bessel.hpp @@ -12,9 +12,13 @@ #include #include +#include #include #include #include #include #include #include +#include +#include +#include diff --git a/include/kyosu/types/impl/detail/bessel_utils.hpp b/include/kyosu/types/impl/detail/bessel_utils.hpp new file mode 100644 index 00000000..a0b2bf31 --- /dev/null +++ b/include/kyosu/types/impl/detail/bessel_utils.hpp @@ -0,0 +1,226 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // bound + //===------------------------------------------------------------------------------------------- + template + auto bound(Z const & z) noexcept + { + using u_t = eve::underlying_type_t; + return abs(z)*u_t(1.1)+2; + } + + //===------------------------------------------------------------------------------------------- + // R use lentz + //===------------------------------------------------------------------------------------------- + template + struct fraction_traits_simple + { + using value_type = typename Gen::result_type; + using result_type = typename Gen::result_type; + using scalar_type = typename Gen::scalar_type; + + static result_type a(const value_type&) noexcept + { + return result_type(1); + } + static result_type b(const value_type& v) noexcept + { + return v; + } + + }; + + template + struct fraction_traits_pair + { + using value_type = typename Gen::value_type; + using result_type = typename kumi::element<0, value_type>; + using scalar_type = typename Gen::scalar_type; + + static result_type a(const value_type& v) noexcept + { + return kumi::get<1>(v); + } + static result_type b(const value_type& v) noexcept + { + return kumi::get<0>(v); + } + }; + + template + struct fraction_traits + : std::conditional_t< + std::same_as, + fraction_traits_simple, + fraction_traits_pair> + { + }; + + template + inline auto lentz_b(Gen& g, const U& eps, size_t const& max_terms = 1000) noexcept + { + using traits = fraction_traits; + using value_type = typename traits::value_type; + using result_type = typename traits::result_type; + using scalar_type = typename traits::scalar_type; + using u_t = scalar_type; + + u_t tiny = 16*eve::smallestposval(eve::as()) ; + u_t terminator(eps); + constexpr auto b_only = std::same_as; + + auto v = g(); + + result_type f, C, delta; + f = traits::b(v); + // if constexpr(b_only) +// { +// f = kyosu::if_else(kyosu::is_eqz(v), tiny, v) ; +// } +// else +// { +// f = kyosu::if_else(kyosu::is_eqz(kumi::get<0>(v)), tiny, kumi::get<0>(v)) ; +// } + C = f; + result_type D{}; //(u_t(0)); + + std::uintmax_t counter(max_terms); + do{ + v = g(); + if constexpr(b_only) + { + D = v + D; + C = v + kyosu::rec(C); + } + else + { + D = kumi::get<0>(v) + kumi::get<1>(v) * D; + C = kumi::get<0>(v) + kumi::get<1>(v) / C; + } + D = kyosu::if_else(kyosu::is_eqz(D), tiny, D); + C = kyosu::if_else(kyosu::is_eqz(C), tiny, C); + D = kyosu::rec(D); + delta = C*D; + f = f * delta; + } while(eve::any((kyosu::abs(kyosu::dec(delta)) > terminator) && --counter)); + g.reset(); + return f; + } + + template + struct R_estimate + { + using result_type = T; + using value_type = T; + using scalar_type = eve::underlying_type_t; + + int n; + T rz; + int i; + int sgn; + + R_estimate(size_t n_, T z_) + : n(n_), rz(kyosu::rec(z_)) + { + reset(); + } + + void reset() + { + i = 0; sgn = -1; + } + + value_type operator()() + { + sgn = -sgn; + ++i; + auto r = sgn*2*(n+i-1)*rz; + return r; + } + }; + + template + inline auto R(size_t n, Z z) noexcept + // compute the ratio Jn(n, z)/Jn(n-1, z)J + { + using u_t = eve::underlying_type_t; + R_estimate re(n, z); + auto r = lentz_b(re, eve::eps(eve::as())); + return r; + } + + template + inline auto bjn(size_t n, Z z) noexcept + // compute log(Jn(z)) + { + Z r{}; + if (n > 0) + { + auto rz = rec(z); + auto rn = R(n, z); + auto pri = rn; + for(int i=n-1; i > 0; --i) + { + rn = 2*(i)*rz-rec(rn); +// std::cout << "i " << i << " ri " << rn << std::endl; +// std::cout << "i " << i << " Ri " << R(i,z) << std::endl; + pri *= rn; + } + auto s = kyosu::cyl_bessel_j0(z)/pri; + return s; + } + else + return cyl_bessel_j0(z); + + } + + template + inline auto arg_adjust(Z z) noexcept + // modify ipart of z modulo 2*pi fit in ]-pi +pi] + { + if constexpr(kyosu::concepts::complex) + { + ipart(z) = eve::to_nearest(eve::rem)(ipart(z), 2*eve::pi(eve::as(ipart(z)))); + return z; + } + else return z; + } + +// template +// inline auto lbesseljn(size_t n, Z z) noexcept +// // compute log(Jn(z)) +// { +// Z r{}, s{}; +// if (n > 0) +// { +// auto rz = rec(z); +// auto rn = R(n, z); +// auto lri = log(rn); +// for(int i=n-1; i > 0; --i) +// { +// rn = 2*(i)*rz-rec(rn); +// // std::cout << "i " << i << " ri " << rn << std::endl; +// // std::cout << "i " << i << " Ri " << R(i,z) << std::endl; +// lri += log(rn); +// } +// s = log(kyosu::cyl_bessel_j0(z))-lri; +// } +// else +// s = kyosu::log(cyl_bessel_j0(z)); +// // ipart(s) = eve::to_nearest(eve::rem)(ipart(s), 2*eve::pi(eve::as(ipart(s)))); +// return arg_adjust(s); +// } + +} diff --git a/include/kyosu/types/impl/detail/lj0.hpp b/include/kyosu/types/impl/detail/lj0.hpp new file mode 100644 index 00000000..bedb8c5e --- /dev/null +++ b/include/kyosu/types/impl/detail/lj0.hpp @@ -0,0 +1,29 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +//#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_lbessel_j0 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return log(cyl_bessel_j0(z)); + } + else + { + return cayley_extend(cyl_bessel_j0, z); + } + } +} diff --git a/include/kyosu/types/impl/detail/lj1.hpp b/include/kyosu/types/impl/detail/lj1.hpp new file mode 100644 index 00000000..bdba02fa --- /dev/null +++ b/include/kyosu/types/impl/detail/lj1.hpp @@ -0,0 +1,29 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_lbessel_j1 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return cyl_bessel_jn(1, z); + } + else + { + return cayley_extend(cyl_bessel_j0, z); + } + } +} diff --git a/include/kyosu/types/impl/detail/ljn.hpp b/include/kyosu/types/impl/detail/ljn.hpp new file mode 100644 index 00000000..0655caaf --- /dev/null +++ b/include/kyosu/types/impl/detail/ljn.hpp @@ -0,0 +1,51 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_lbessel_jn + //===------------------------------------------------------------------------------------------- + + template + auto dispatch(eve::tag_of, N n, Z z) noexcept + { + if constexpr(concepts::complex ) + { + if ( is_eqz(n) ) + { + return cyl_lbessel_j0(z); + } + else + { + Z r{}, s{}; + auto rz = rec(z); + auto rn = R(n, z); + auto lri = log(rn); + for(int i=n-1; i > 0; --i) + { + rn = 2*(i)*rz-rec(rn); +// std::cout << "i " << i << " ri " << rn << std::endl; +// std::cout << "i " << i << " Ri " << R(i,z) << std::endl; + lri += log(rn); + } + s = log(kyosu::cyl_lbessel_j0(z))-lri; + return arg_adjust(s); + } + } + else + { + return cayley_extend_rev(cyl_lbessel_jn, n, z); + } + } + +} diff --git a/test/doc/cyl_lbessel_jn.cpp b/test/doc/cyl_lbessel_jn.cpp new file mode 100644 index 00000000..d09bed00 --- /dev/null +++ b/test/doc/cyl_lbessel_jn.cpp @@ -0,0 +1,27 @@ +#include +#include +#include +//#include + + +int main() +{ + std::cout.precision(16); + // using w_t = eve::wide>; + //auto z = kyosu::complex(w_t(20.0, 1.0)); + auto z = kyosu::complex(1.0); + for(int n=0; n < 10; ++n) + { +// auto jj = kyosu::log(kyosu::cyl_bessel_jn(n, z)); +// std::cout << jj << std::endl; +// auto val1 = kyosu::_::lbesseljn(n, z); +// std::cout << val1 << std::endl; +// } +// auto n = 2; + // auto j0 = kyosu::cyl_bessel_j0(z); + auto jn = kyosu::cyl_lbessel_jn(n, z); + std::cout << jn << std::endl; + std::cout << kyosu::log(kyosu::cyl_bessel_jn(n, z)) << std::endl; + } + return 0; +} From 305296c173093f6969624d0ff47cc70a7e36a3bf Mon Sep 17 00:00:00 2001 From: jtlap Date: Thu, 26 Oct 2023 22:27:42 +0200 Subject: [PATCH 02/17] cayley/log_bessel --- .../kyosu/types/impl/detail/bessel_utils.hpp | 48 +++++++++---------- include/kyosu/types/impl/detail/ljn.hpp | 42 +++++++++++----- test/doc/cyl_lbessel_jn.cpp | 21 ++++---- 3 files changed, 63 insertions(+), 48 deletions(-) diff --git a/include/kyosu/types/impl/detail/bessel_utils.hpp b/include/kyosu/types/impl/detail/bessel_utils.hpp index a0b2bf31..3bf7dc9a 100644 --- a/include/kyosu/types/impl/detail/bessel_utils.hpp +++ b/include/kyosu/types/impl/detail/bessel_utils.hpp @@ -198,29 +198,29 @@ namespace kyosu::_ else return z; } -// template -// inline auto lbesseljn(size_t n, Z z) noexcept -// // compute log(Jn(z)) -// { -// Z r{}, s{}; -// if (n > 0) -// { -// auto rz = rec(z); -// auto rn = R(n, z); -// auto lri = log(rn); -// for(int i=n-1; i > 0; --i) -// { -// rn = 2*(i)*rz-rec(rn); -// // std::cout << "i " << i << " ri " << rn << std::endl; -// // std::cout << "i " << i << " Ri " << R(i,z) << std::endl; -// lri += log(rn); -// } -// s = log(kyosu::cyl_bessel_j0(z))-lri; -// } -// else -// s = kyosu::log(cyl_bessel_j0(z)); -// // ipart(s) = eve::to_nearest(eve::rem)(ipart(s), 2*eve::pi(eve::as(ipart(s)))); -// return arg_adjust(s); -// } + template + inline auto lbesseljn(size_t n, Z z) noexcept + // compute log(Jn(z)) + { + Z r{}, s{}; + if (n > 0) + { + auto rz = rec(z); + auto rn = R(n, z); + auto lri = log(rn); + for(int i=n-1; i > 0; --i) + { + rn = 2*(i)*rz-rec(rn); + std::cout << "i " << i << " ri " << rn << std::endl; + std::cout << "i " << i << " Ri " << R(i,z) << std::endl; + lri += log(rn); + } + s = log(kyosu::cyl_bessel_j0(z))-lri; + } + else + s = kyosu::log(cyl_bessel_j0(z)); + // ipart(s) = eve::to_nearest(eve::rem)(ipart(s), 2*eve::pi(eve::as(ipart(s)))); + return arg_adjust(s); + } } diff --git a/include/kyosu/types/impl/detail/ljn.hpp b/include/kyosu/types/impl/detail/ljn.hpp index 0655caaf..422fdb26 100644 --- a/include/kyosu/types/impl/detail/ljn.hpp +++ b/include/kyosu/types/impl/detail/ljn.hpp @@ -21,26 +21,46 @@ namespace kyosu::_ { if constexpr(concepts::complex ) { - if ( is_eqz(n) ) + Z r{}, s{}; + if (n > 0) { - return cyl_lbessel_j0(z); - } - else - { - Z r{}, s{}; auto rz = rec(z); auto rn = R(n, z); auto lri = log(rn); for(int i=n-1; i > 0; --i) { - rn = 2*(i)*rz-rec(rn); -// std::cout << "i " << i << " ri " << rn << std::endl; -// std::cout << "i " << i << " Ri " << R(i,z) << std::endl; + rn = 2*i*rz-rec(rn); +// std::cout << "i " << i << " ri " << rn << std::endl; +// std::cout << "i " << i << " Ri " << R(i,z) << std::endl; lri += log(rn); } - s = log(kyosu::cyl_lbessel_j0(z))-lri; - return arg_adjust(s); + s = log(kyosu::cyl_bessel_j0(z))-lri; } + else s = kyosu::log(cyl_bessel_j0(z)); + return arg_adjust(s); + + +// return lbesseljn(n, z); +// if ( is_eqz(n) ) +// { +// return cyl_lbessel_j0(z); +// } +// else +// { +// Z r{}, s{}; +// auto rz = rec(z); +// auto rn = R(n, z); +// auto lri = log(rn); +// for(int i=n-1; i > 0; --i) +// { +// rn = 2*(i)*rz-rec(rn); +// std::cout << "i " << i << " ri " << rn << std::endl; +// std::cout << "i " << i << " Ri " << R(i,z) << std::endl; +// lri += log(rn); +// } +// s = log(kyosu::cyl_lbessel_j0(z))-lri; +// return arg_adjust(s); +// } } else { diff --git a/test/doc/cyl_lbessel_jn.cpp b/test/doc/cyl_lbessel_jn.cpp index d09bed00..0205c20b 100644 --- a/test/doc/cyl_lbessel_jn.cpp +++ b/test/doc/cyl_lbessel_jn.cpp @@ -7,21 +7,16 @@ int main() { std::cout.precision(16); - // using w_t = eve::wide>; - //auto z = kyosu::complex(w_t(20.0, 1.0)); - auto z = kyosu::complex(1.0); - for(int n=0; n < 10; ++n) + using w_t = eve::wide>; + auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); +// auto z = kyosu::complex(1.0); + for(int n=0; n <= 10; ++n) { -// auto jj = kyosu::log(kyosu::cyl_bessel_jn(n, z)); -// std::cout << jj << std::endl; -// auto val1 = kyosu::_::lbesseljn(n, z); -// std::cout << val1 << std::endl; -// } -// auto n = 2; - // auto j0 = kyosu::cyl_bessel_j0(z); auto jn = kyosu::cyl_lbessel_jn(n, z); - std::cout << jn << std::endl; - std::cout << kyosu::log(kyosu::cyl_bessel_jn(n, z)) << std::endl; + std::cout << "n " << n << " z " << z << std::endl; + std::cout << jn << std::endl; + std::cout << kyosu::log(kyosu::cyl_bessel_jn(n, z)) << std::endl; +// std::cout << kyosu::_::lbesseljn(n, z) << std::endl; } return 0; } From 4557daedbb8a8215d212c2c263152993edaeea27 Mon Sep 17 00:00:00 2001 From: jtlap Date: Sat, 28 Oct 2023 16:40:14 +0200 Subject: [PATCH 03/17] tests for ljn lj0 and lj1 --- include/kyosu/functions/cyl_lbessel_j0.hpp | 6 +- include/kyosu/types/impl/detail/lj1.hpp | 4 +- include/kyosu/types/impl/detail/ljn.hpp | 28 +------ test/doc/cyl_lbessel_j0.cpp | 33 ++++++++ test/unit/function/cyl_lbessel_j0.cpp | 85 ++++++++++++++++++++ test/unit/function/cyl_lbessel_j1.cpp | 83 ++++++++++++++++++++ test/unit/function/cyl_lbessel_jn.cpp | 90 ++++++++++++++++++++++ 7 files changed, 298 insertions(+), 31 deletions(-) create mode 100644 test/doc/cyl_lbessel_j0.cpp create mode 100644 test/unit/function/cyl_lbessel_j0.cpp create mode 100644 test/unit/function/cyl_lbessel_j1.cpp create mode 100644 test/unit/function/cyl_lbessel_jn.cpp diff --git a/include/kyosu/functions/cyl_lbessel_j0.hpp b/include/kyosu/functions/cyl_lbessel_j0.hpp index 24ed12c9..0e1727de 100644 --- a/include/kyosu/functions/cyl_lbessel_j0.hpp +++ b/include/kyosu/functions/cyl_lbessel_j0.hpp @@ -8,7 +8,7 @@ #pragma once #include -#include + namespace kyosu::tags { @@ -21,8 +21,8 @@ namespace kyosu::tags template static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { - - return log(cyl_bessel_j0(complex(v))); + auto fn = callable_cyl_lbessel_j0{}; + return fn(complex(v)); } template diff --git a/include/kyosu/types/impl/detail/lj1.hpp b/include/kyosu/types/impl/detail/lj1.hpp index bdba02fa..aa0db7b3 100644 --- a/include/kyosu/types/impl/detail/lj1.hpp +++ b/include/kyosu/types/impl/detail/lj1.hpp @@ -19,11 +19,11 @@ namespace kyosu::_ { if constexpr(concepts::complex ) { - return cyl_bessel_jn(1, z); + return cyl_lbessel_jn(1, z); } else { - return cayley_extend(cyl_bessel_j0, z); + return cayley_extend(cyl_bessel_j1, z); } } } diff --git a/include/kyosu/types/impl/detail/ljn.hpp b/include/kyosu/types/impl/detail/ljn.hpp index 422fdb26..2b8f560a 100644 --- a/include/kyosu/types/impl/detail/ljn.hpp +++ b/include/kyosu/types/impl/detail/ljn.hpp @@ -30,37 +30,13 @@ namespace kyosu::_ for(int i=n-1; i > 0; --i) { rn = 2*i*rz-rec(rn); -// std::cout << "i " << i << " ri " << rn << std::endl; -// std::cout << "i " << i << " Ri " << R(i,z) << std::endl; lri += log(rn); } s = log(kyosu::cyl_bessel_j0(z))-lri; } else s = kyosu::log(cyl_bessel_j0(z)); - return arg_adjust(s); - - -// return lbesseljn(n, z); -// if ( is_eqz(n) ) -// { -// return cyl_lbessel_j0(z); -// } -// else -// { -// Z r{}, s{}; -// auto rz = rec(z); -// auto rn = R(n, z); -// auto lri = log(rn); -// for(int i=n-1; i > 0; --i) -// { -// rn = 2*(i)*rz-rec(rn); -// std::cout << "i " << i << " ri " << rn << std::endl; -// std::cout << "i " << i << " Ri " << R(i,z) << std::endl; -// lri += log(rn); -// } -// s = log(kyosu::cyl_lbessel_j0(z))-lri; -// return arg_adjust(s); -// } + using e_t = eve::element_type_t; + return if_else(is_nez(z), arg_adjust(s), kyosu::complex(eve::minf(eve::as()))); } else { diff --git a/test/doc/cyl_lbessel_j0.cpp b/test/doc/cyl_lbessel_j0.cpp new file mode 100644 index 00000000..bc1b376d --- /dev/null +++ b/test/doc/cyl_lbessel_j0.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::cyl_lbessel_j0; + std::cout.precision(16); + + std::cout << "1.0 " << " -> " << cyl_lbessel_j0(1.0) << "\n"; + std::cout << "1+0i " << " -> " << cyl_lbessel_j0(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << cyl_lbessel_j0(15.0) << "\n"; + std::cout << "15+0i " << " -> " << cyl_lbessel_j0(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << cyl_lbessel_j0(40.0) << "\n"; + std::cout << "40+0i " << " -> " << cyl_lbessel_j0(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << cyl_lbessel_j0(60.0) << "\n"; + std::cout << "60+0i " << " -> " << cyl_lbessel_j0(kyosu::complex(60.0, 0.0)) << "\n"; + + eve::wide> z(1.0, 15.0, 40.0, 60.0); + auto zz = kyosu::complex(z); + std::cout << z << " -> \n" << cyl_lbessel_j0(z) << "\n"; + std::cout << zz << " -> \n" << cyl_lbessel_j0(zz) << "\n"; + + std::cout << "1.0 " << " -> " << cyl_lbessel_j0(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << cyl_lbessel_j0(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_lbessel_j0(kyosu::quaternion(1.0, 36.0, 2.0, 1.5)) << "\n"; + + + eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); + auto z0 = kyosu::complex(z1, z2); + std::cout << z0 << "\n -> " << cyl_lbessel_j0(z0) << "\n"; + return 0; +} diff --git a/test/unit/function/cyl_lbessel_j0.cpp b/test/unit/function/cyl_lbessel_j0.cpp new file mode 100644 index 00000000..8f67b136 --- /dev/null +++ b/test/unit/function/cyl_lbessel_j0.cpp @@ -0,0 +1,85 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::abs over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + using e_t = eve::element_type_t; + if constexpr(sizeof(e_t) == 8) + { + std::cout.precision(17); + std::array re{ e_t(-9.8403264736233908e+01), e_t(8.6287324687322851e+01), e_t(2.8677183208357349e+01), e_t(7.5642931882481321e+01), e_t(9.6351793077594323e+01) + , e_t(-7.8716695947540288e+01), e_t(-5.0019457963533974e+01), e_t(5.9030075032843811e+01), e_t(-7.5596124248257235e+00), e_t(8.5810875418338028e+01)}; + + std::array im{ e_t(-1.9271037994246587e+01), e_t(7.9043790847200569e+01), e_t(-3.0599555605285602e+01), e_t(-7.1720588923445590e+01), e_t(-6.1110206702938363e+01) + , e_t(9.0167523214068439e+01), e_t(-7.0131934820406244e+01), e_t(-9.3405497312943496e+01), e_t(8.3431692318213322e+01), e_t(-4.1557438097347173e+01)}; + + std::array reres{e_t(-8.8414975251226071e+06), e_t(-3.9177291392384149e+32), e_t(-1.2012865579185671e+12), e_t(5.4246749994831845e+29), e_t(-3.8201314596721725e+23) + , e_t(-5.1783207689008547e+37), e_t(1.0487044863432896e+29), e_t(-8.2091088236413212e+38), e_t(2.4924320903438798e+34), e_t(-4.1577568667728512e+16)}; + std::array imres{e_t(2.9689225084136692e+06), e_t(6.8130876628410661e+32), e_t(-2.9179495183191563e+10), e_t(-8.8752530111744507e+28), e_t(1.2946830695756804e+25) + , e_t(9.5653949361454081e+36), e_t(6.5259614844547437e+28), e_t(1.1301644890532119e+39), e_t(7.0530420820017420e+34), e_t(-1.8907103618909208e+16)}; + for(int i=0; i < 10; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::log(kyosu::complex(reres[i], imres[i])); + TTS_RELATIVE_EQUAL(kyosu::cyl_lbessel_j0(c), res, 1.0e-4) << i << " <- " << c << '\n'; + } + } + std::cout.precision(17); + std::array re{2.175102683907268e+00, 3.007738475083643e-01, 3.571791197760752e+00, 1.916834892035037e+00, 2.309329011574615e+00, 1.116302757926106e+00, 3.454188251414261e+00, 2.504516854443353e+00}; + std::array im{ 3.008362577802520e+00, 3.733804275914350e+00, 2.988312632706299e+00, 4.411110052057489e+00, 3.055249958745497e+00, 2.508994160634028e+00, 3.100976745146938e-01, 3.854880352108083e+00}; + + std::array reres{-1.1271605877499498e+00, 8.6819450415271469e+00, -3.7342482747514563e+00, -1.9594143015933545e+00, -1.6842573645694947e+00, 1.9975730077018730e+00, -3.9366090510718693e-01, -5.2795374955460659e+00}; + std::array imres{-4.1516002202378823e+00, -2.2772457219968798e+00, -1.0734738341492034e-01, -1.5246293957752114e+01, -4.0999678002136202e+00, -2.3878228304877860e+00, -4.8664505025944993e-02, -7.2799291393660077e+00}; + for(int i=0; i < 8; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::log(kyosu::complex(reres[i], imres[i])); + double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; + TTS_RELATIVE_EQUAL(kyosu::cyl_lbessel_j0(c), res, tol) << i << " <- " << c << '\n'; + } +}; + + +TTS_CASE_WITH ( "Check kyosu::abs over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + if constexpr(sizeof(T) == 8) + { + using e_t = eve::element_type_t; + double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; + auto c = kyosu::complex(a0, a1); + auto cb= kyosu::conj(c); + auto cm= -c; + using kyosu::cyl_lbessel_j0; + using kyosu::cyl_bessel_j0; + auto lj0c = cyl_lbessel_j0(c); + auto logj0c = kyosu::log(cyl_bessel_j0(c)); + TTS_RELATIVE_EQUAL(lj0c, logj0c, tol) << "c " << c << "\n"; + + auto lj0cb = cyl_lbessel_j0(cb); + auto logj0cb = kyosu::log(cyl_bessel_j0(cb)); + TTS_RELATIVE_EQUAL(lj0cb, logj0cb, tol) << "c " << c << "\n"; + + auto lj0cm = cyl_lbessel_j0(cm); + auto logj0cm = kyosu::log(cyl_bessel_j0(cm)); + TTS_RELATIVE_EQUAL(lj0cm, logj0cm, tol) << "c " << c << "\n"; + auto z = kyosu::complex(T(0), T(0)); + TTS_IEEE_EQUAL(cyl_lbessel_j0(z), kyosu::complex(eve::zero(eve::as()))); + } +}; diff --git a/test/unit/function/cyl_lbessel_j1.cpp b/test/unit/function/cyl_lbessel_j1.cpp new file mode 100644 index 00000000..2c71a74a --- /dev/null +++ b/test/unit/function/cyl_lbessel_j1.cpp @@ -0,0 +1,83 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::abs over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + using e_t = eve::element_type_t; + if constexpr(sizeof(e_t) == 8) + { + std::cout.precision(17); + std::array re{ -9.8403264736233908e+01, 8.6287324687322851e+01, 2.8677183208357349e+01, 7.5642931882481321e+01, 9.6351793077594323e+01 + , -7.8716695947540288e+01, -5.0019457963533974e+01, 5.9030075032843811e+01, -7.5596124248257235e+00, 8.5810875418338028e+01}; + std::array im{-1.9271037994246587e+01, 7.9043790847200569e+01, -3.0599555605285602e+01, -7.1720588923445590e+01, -6.1110206702938363e+01 + , 9.0167523214068439e+01, -7.0131934820406244e+01, -9.3405497312943496e+01, 8.3431692318213322e+01, -4.1557438097347173e+01}; + + std::array reres{3.0094177278540167e+06, -6.8058089697014938e+32, -3.8806283458664719e+10, -8.6565181213118240e+28, 1.2915079982169630e+25, + -9.3925711342553668e+36, 6.4595726071319036e+28, 1.1238442153777222e+39, -7.0123433800633681e+34, -1.9060735303763808e+16}; + std::array imres{8.8185850069141667e+06, -3.8848970122848989e+32, 1.1905919032156870e+12, -5.4098751798106477e+29, 4.2914274806719309e+23, + -5.1646555066948327e+37, -1.0459530606430648e+29, 8.2050979033727447e+38, 2.4737508957685513e+34, 4.1393451863209432e+16}; + for(int i=0; i < 10; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::log(kyosu::complex(reres[i], imres[i])); + TTS_RELATIVE_EQUAL(kyosu::cyl_lbessel_j1(c), res, 1.0e-7) << i << " <- " << c << '\n'; + } + } + std::cout.precision(17); + std::array re{2.175102683907268e+00, 3.007738475083643e-01, 3.571791197760752e+00, 1.916834892035037e+00, 2.309329011574615e+00, 1.116302757926106e+00, 3.454188251414261e+00, 2.504516854443353e+00}; + std::array im{ 3.008362577802520e+00, 3.733804275914350e+00, 2.988312632706299e+00, 4.411110052057489e+00, 3.055249958745497e+00, 2.508994160634028e+00, 3.100976745146938e-01, 3.854880352108083e+00}; + + std::array reres{3.6097545599299043e+00, 2.0607036224436954e+00, -2.1724057594208968e-01, 1.3644664632775582e+01, 3.5380568884368788e+00, 2.1507545829381751e+00, 1.5764037998430128e-01, 6.2619292962097814e+00}; + std::array imres{-1.3891357467781533e+00, 7.3810832876707710e+00, -3.4818290031753225e+00, -2.4845686655182946e+00, -1.8826359832772270e+00, 1.3851339474357725e+00, -1.3157891968242472e-01, -5.2738074786282327e+00}; + for(int i=0; i < 8; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::log(kyosu::complex(reres[i], imres[i])); + double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; + TTS_RELATIVE_EQUAL(kyosu::cyl_lbessel_j1(c), res, tol) << i << " <- " << c << '\n'; + } +}; + +// TTS_CASE_WITH ( "Check kyosu::cyl_lbessel_j1 over real" +// , kyosu::real_types +// , tts::generate(tts::randoms(-10,10), +// tts::randoms(-10,10) +// ) +// ) +// (T a0, T a1) +// { +// if constexpr(sizeof(T) == 8) +// { +// using e_t = eve::element_type_t; +// double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; +// auto c = kyosu::complex(a0, a1); +// auto cb= kyosu::conj(c); +// auto cm= -c; +// using kyosu::cyl_lbessel_j1; +// using kyosu::cyl_bessel_j1; +// auto lj1c = cyl_lbessel_j1(c); +// auto logj1c = kyosu::log(cyl_bessel_j1(c)); +// TTS_RELATIVE_EQUAL(lj1c, logj1c, tol) << "c " << c << "\n"; + +// auto lj1cb = cyl_lbessel_j1(cb); +// auto logj1cb = kyosu::log(cyl_bessel_j1(cb)); +// TTS_RELATIVE_EQUAL(lj1cb, logj1cb, tol) << "c " << c << "\n"; + +// auto lj1cm = cyl_lbessel_j1(cm); +// auto logj1cm = kyosu::log(cyl_bessel_j1(cm)); +// TTS_RELATIVE_EQUAL(lj1cm, logj1cm, tol) << "c " << c << "\n"; +// auto z = kyosu::complex(T(0), T(0)); +// TTS_IEEE_EQUAL(cyl_lbessel_j1(z), kyosu::complex(eve::minf(eve::as()))); +// } +// }; diff --git a/test/unit/function/cyl_lbessel_jn.cpp b/test/unit/function/cyl_lbessel_jn.cpp new file mode 100644 index 00000000..1696861a --- /dev/null +++ b/test/unit/function/cyl_lbessel_jn.cpp @@ -0,0 +1,90 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::abs over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + if constexpr(sizeof(T) == 8) + { + std::cout.precision(17); + std::array re{ -9.8403264736233908e+01, 8.6287324687322851e+01, 2.8677183208357349e+01, 7.5642931882481321e+01, 9.6351793077594323e+01 + , -7.8716695947540288e+01, -5.0019457963533974e+01, 5.9030075032843811e+01, -7.5596124248257235e+00, 8.5810875418338028e+01}; + std::array im{-1.9271037994246587e+01, 7.9043790847200569e+01, -3.0599555605285602e+01, -7.1720588923445590e+01, -6.1110206702938363e+01 + , 9.0167523214068439e+01, -7.0131934820406244e+01, -9.3405497312943496e+01, 8.3431692318213322e+01, -4.1557438097347173e+01}; + + std::array reres{-3.3279158735087640e+06, 6.7446373092398792e+32, 1.0973422267168294e+11, 6.9510412682541761e+28, -1.2657523599731811e+25, + 8.0429582124849094e+36, -5.9430267579521807e+28, -1.0741846420111567e+39, 6.6950456879149885e+34, 2.0245919743884628e+16}; + std::array imres{-8.6289796115820184e+06, 3.6264735895350172e+32, -1.1056118662675706e+12, 5.2905159458975805e+29, -7.9797425680943093e+23, + 5.0551834834322169e+37, 1.0236118318065309e+29, -8.1681100953492296e+38, -2.3290495380992445e+34, -3.9909835121233848e+16}; + for(int i=0; i < 10; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_jn(3, c), res, 1.0e-7) << i << " <- " << c << '\n'; + } + } + + std::array re{-2.2013892240657134e+00, -2.8382367527759522e-01, -7.0737598522551182e-01, 8.6178612999047166e-02, 1.8217799341954277e-01, -7.6718326053284047e-01, -5.9167207480492257e-01, 7.5911892713725948e-01}; + + std::array im{8.7228967955690240e-01, -2.1140894389723703e+00, 2.2905192466959483e+00, -1.5275122978034243e+00, 5.5435596674685006e-01, -1.1418929858153533e+00, 1.1666506779144503e+00, -1.9198575359478032e+00}; + + std::array reres{-1.3672173596226031e-01, 1.2004993008474076e-01, 3.4070032613198098e-01, -1.5845109266378627e-02, -3.4617874091075179e-03, 5.4024365099246512e-02, 4.9541047773116346e-02, -2.1467277569289250e-01}; + + std::array imres{1.6497274394928027e-01, 2.3540711137379494e-01, -1.7769121638069371e-01, 8.4636178354763933e-02, -2.3973187185149024e-03, -1.7452629167850800e-02, -3.8485807736362864e-03, 5.8435910432745555e-02}; + using e_t = eve::element_type_t; + + for(size_t i=0; i < re.size(); ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::log(kyosu::complex(reres[i], imres[i])); + double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; + TTS_RELATIVE_EQUAL(kyosu::cyl_lbessel_jn(3, c), res, tol) << i << " <- " << c << '\n'; + } +}; + + +TTS_CASE_WITH ( "Check kyosu::abs over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + if constexpr(sizeof(T) == 8) + { + using e_t = eve::element_type_t; + double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; + auto c = kyosu::complex(a0, a1); + auto cb= kyosu::conj(c); + auto cm= -c; + auto cr = kyosu::complex(eve::abs(a0)); + using kyosu::cyl_lbessel_jn; + using kyosu::cyl_bessel_jn; + auto ljnc = cyl_lbessel_jn(3, c); + auto logjnc = kyosu::log(cyl_bessel_jn(3, c)); + TTS_RELATIVE_EQUAL(ljnc, logjnc, tol) << "c " << c << "\n"; + + auto ljncb = cyl_lbessel_jn(3, cb); + auto logjncb = kyosu::log(cyl_bessel_jn(3, cb)); + TTS_RELATIVE_EQUAL(ljncb, logjncb, tol) << "c " << c << "\n"; + + auto ljncm = cyl_lbessel_jn(3, cm); + auto logjncm = kyosu::log(cyl_bessel_jn(3, cm)); + TTS_RELATIVE_EQUAL(ljncm, logjncm, tol) << "c " << c << "\n"; + TTS_EXPECT(eve::all(kyosu::is_real(kyosu::cyl_lbessel_jn(3, cr)))); +// TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + auto z = kyosu::complex(T(0), T(0)); + TTS_IEEE_EQUAL(cyl_lbessel_jn(3, z), kyosu::complex(eve::minf(eve::as()))); + } +}; From 210923821332d2c7b8bb14e2b64d6265b0c1d0e6 Mon Sep 17 00:00:00 2001 From: jtlap Date: Sat, 28 Oct 2023 23:37:57 +0200 Subject: [PATCH 04/17] i0 i1 in --- include/kyosu/functions.hpp | 3 + include/kyosu/functions/cyl_lbessel_i0.hpp | 80 ++++++++++++++++++++++ include/kyosu/functions/cyl_lbessel_i1.hpp | 80 ++++++++++++++++++++++ include/kyosu/functions/cyl_lbessel_in.hpp | 79 +++++++++++++++++++++ include/kyosu/types/impl/bessel.hpp | 3 + include/kyosu/types/impl/detail/li0.hpp | 29 ++++++++ include/kyosu/types/impl/detail/li1.hpp | 29 ++++++++ include/kyosu/types/impl/detail/lin.hpp | 38 ++++++++++ test/doc/cyl_lbessel_i0.cpp | 18 +++++ test/doc/cyl_lbessel_i1.cpp | 18 +++++ test/doc/cyl_lbessel_in.cpp | 27 ++++++++ test/doc/cyl_lbessel_j1.cpp | 19 +++++ 12 files changed, 423 insertions(+) create mode 100644 include/kyosu/functions/cyl_lbessel_i0.hpp create mode 100644 include/kyosu/functions/cyl_lbessel_i1.hpp create mode 100644 include/kyosu/functions/cyl_lbessel_in.hpp create mode 100644 include/kyosu/types/impl/detail/li0.hpp create mode 100644 include/kyosu/types/impl/detail/li1.hpp create mode 100644 include/kyosu/types/impl/detail/lin.hpp create mode 100644 test/doc/cyl_lbessel_i0.cpp create mode 100644 test/doc/cyl_lbessel_i1.cpp create mode 100644 test/doc/cyl_lbessel_in.cpp create mode 100644 test/doc/cyl_lbessel_j1.cpp diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index 537526b3..fb179499 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -56,6 +56,9 @@ #include #include #include +#include +#include +#include #include #include #include diff --git a/include/kyosu/functions/cyl_lbessel_i0.hpp b/include/kyosu/functions/cyl_lbessel_i0.hpp new file mode 100644 index 00000000..74fcc6ff --- /dev/null +++ b/include/kyosu/functions/cyl_lbessel_i0.hpp @@ -0,0 +1,80 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_lbessel_i0: eve::elementwise + { + using callable_tag_type = callable_cyl_lbessel_i0; + + KYOSU_DEFERS_CALLABLE(cyl_lbessel_i0_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept + { + callable_cyl_lbessel_i0 fn{}; + return fn(complex(v)); + } + + template + KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) + { + return eve::tag_invoke(*this, target); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_lbessel_i0 +//! @brief Computes the logarithm of the Bessel function of the first kind, +//! \f$ J_1(x)=\frac1{\pi }\int _{0}^{\pi}\cos(\tau-x\sin \tau )\,\mathrm {d} \tau \f$ +//! extended to the complex plane and cayley_dickson values. +//! +//! \f$ J_1(\f$ is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_lbessel_i0(T z) noexcept; +//! template constexpr T cyl_lbessel_i0(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return the cylindrical \f$J_0(z)\f$. If T is a floating point value (real), a complex value is returned. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_lbessel_i0.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_lbessel_i0 cyl_lbessel_i0 = {}; +} diff --git a/include/kyosu/functions/cyl_lbessel_i1.hpp b/include/kyosu/functions/cyl_lbessel_i1.hpp new file mode 100644 index 00000000..636de4e9 --- /dev/null +++ b/include/kyosu/functions/cyl_lbessel_i1.hpp @@ -0,0 +1,80 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_lbessel_i1: eve::elementwise + { + using callable_tag_type = callable_cyl_lbessel_i1; + + KYOSU_DEFERS_CALLABLE(cyl_lbessel_i1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept + { + callable_cyl_lbessel_i1 fn{}; + return fn(complex(v)); + } + + template + KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) + { + return eve::tag_invoke(*this, target); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_lbessel_i1 +//! @brief Computes the logarithm of the Bessel function of the first kind, +//! \f$ J_1(x)=\frac1{\pi }\int _{0}^{\pi}\cos(\tau-x\sin \tau )\,\mathrm {d} \tau \f$ +//! extended to the complex plane and cayley_dickson values. +//! +//! \f$ J_1(\f$ is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_lbessel_i1(T z) noexcept; +//! template constexpr T cyl_lbessel_i1(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return the cylindrical \f$J_0(z)\f$. If T is a floating point value (real), a complex value is returned. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_lbessel_i1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_lbessel_i1 cyl_lbessel_i1 = {}; +} diff --git a/include/kyosu/functions/cyl_lbessel_in.hpp b/include/kyosu/functions/cyl_lbessel_in.hpp new file mode 100644 index 00000000..748d6741 --- /dev/null +++ b/include/kyosu/functions/cyl_lbessel_in.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_cyl_lbessel_in: eve::elementwise + { + using callable_tag_type = callable_cyl_lbessel_in; + + KYOSU_DEFERS_CALLABLE(cyl_lbessel_in_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& v) noexcept { + callable_cyl_lbessel_in fn{}; + return 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 +//! @{ +//! @brief Computes the logarithm of the Bessel functions of the first kind, +//! \f$ I_{n}(z)\f$. +//! +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_lbessel_in(int n, T z) noexcept; +//! template constexpr T cyl_lbessel_in(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return the logarithm of \f$I_n(z)\f$. +//! +//! @warning Up to now the cayley_dickson versions have only been implemented for scalar int values of n. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_lbessel_in.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_lbessel_in cyl_lbessel_in = {}; +} diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp index b443efd4..a6f54575 100644 --- a/include/kyosu/types/impl/bessel.hpp +++ b/include/kyosu/types/impl/bessel.hpp @@ -22,3 +22,6 @@ #include #include #include +#include +#include +#include diff --git a/include/kyosu/types/impl/detail/li0.hpp b/include/kyosu/types/impl/detail/li0.hpp new file mode 100644 index 00000000..fea51dad --- /dev/null +++ b/include/kyosu/types/impl/detail/li0.hpp @@ -0,0 +1,29 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_lbessel_i0 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return cyl_lbessel_in(0, z); + } + else + { + return cayley_extend(cyl_lbessel_i0, z); + } + } +} diff --git a/include/kyosu/types/impl/detail/li1.hpp b/include/kyosu/types/impl/detail/li1.hpp new file mode 100644 index 00000000..a8d586b6 --- /dev/null +++ b/include/kyosu/types/impl/detail/li1.hpp @@ -0,0 +1,29 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_lbessel_j1 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + if constexpr(concepts::complex ) + { + return cyl_lbessel_in(1, z); + } + else + { + return cayley_extend(cyl_bessel_i1, z); + } + } +} diff --git a/include/kyosu/types/impl/detail/lin.hpp b/include/kyosu/types/impl/detail/lin.hpp new file mode 100644 index 00000000..320fd544 --- /dev/null +++ b/include/kyosu/types/impl/detail/lin.hpp @@ -0,0 +1,38 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_in + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, N n, Z z) noexcept + { + if constexpr(concepts::complex ) + { + using e_t = as_real_t; + auto logmiton = [](N n){ + auto nmod4 = n%4; + if (nmod4 == 0) return complex(eve::zero(eve::as())); + else if (nmod4 == 1) return complex(eve::zero(eve::as()), -eve::pio_2(eve::as())); + else if (nmod4 == 2) return complex(eve::zero(eve::as()), eve::pi(eve::as())); + else return complex(eve::zero(eve::as()), eve::pio_2(eve::as())); + }; + + return arg_adjust(logmiton(n)+cyl_lbessel_jn(n,complex(-ipart(z), real(z)))); + } + else + { + return cayley_extend_rev(cyl_bessel_in, n, z); + } + } +} diff --git a/test/doc/cyl_lbessel_i0.cpp b/test/doc/cyl_lbessel_i0.cpp new file mode 100644 index 00000000..44d6c29b --- /dev/null +++ b/test/doc/cyl_lbessel_i0.cpp @@ -0,0 +1,18 @@ +#include +#include +#include + +int main() +{ + std::cout.precision(16); + using w_t = eve::wide>; + auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); +// auto z = kyosu::complex(30.0, 0.0); +// auto zi= kyosu::complex(00.0, 30.0); + // auto z = kyosu::complex(1.0); + auto in = kyosu::cyl_lbessel_i0(z); + std::cout << "0 " << in << std::endl; + auto login = kyosu::log(kyosu::cyl_bessel_i0(z)); + std::cout << "1 " << login << std::endl; + return 0; +} diff --git a/test/doc/cyl_lbessel_i1.cpp b/test/doc/cyl_lbessel_i1.cpp new file mode 100644 index 00000000..180ff649 --- /dev/null +++ b/test/doc/cyl_lbessel_i1.cpp @@ -0,0 +1,18 @@ +#include +#include +#include + +int main() +{ + std::cout.precision(16); + using w_t = eve::wide>; + auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); +// auto z = kyosu::complex(30.0, 0.0); +// auto zi= kyosu::complex(00.0, 30.0); + // auto z = kyosu::complex(1.0); + auto in = kyosu::cyl_lbessel_i1(z); + std::cout << "0 " << in << std::endl; + auto login = kyosu::log(kyosu::cyl_bessel_i1(z)); + std::cout << "1 " << login << std::endl; + return 0; +} diff --git a/test/doc/cyl_lbessel_in.cpp b/test/doc/cyl_lbessel_in.cpp new file mode 100644 index 00000000..c414fb83 --- /dev/null +++ b/test/doc/cyl_lbessel_in.cpp @@ -0,0 +1,27 @@ +#include +#include +#include + +int main() +{ + std::cout.precision(16); + // using w_t = eve::wide>; + //auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); + auto z = kyosu::complex(30.0, 0.0); + auto zi= kyosu::complex(00.0, 30.0); +// auto z = kyosu::complex(1.0); + for(int n=2; n <= 2; ++n) + { + auto in = kyosu::cyl_lbessel_in(n, z); + std::cout << "n " << n << " z " << z << std::endl; + std::cout << "0 " << in << std::endl; + std::cout << "1 " << kyosu::cyl_bessel_in(n, z) << std::endl; + std::cout << "2 " << kyosu::cyl_bessel_jn(n, zi) << std::endl; + std::cout << "3 " << kyosu::cyl_bessel_in(n, 20.0) << std::endl; + auto login = kyosu::log(kyosu::cyl_bessel_in(n, z)); + std::cout << "4 " << login << std::endl; + std::cout << "5 " << kyosu::exp(in) << std::endl; + std::cout << "6 " << kyosu::exp(login) << std::endl; + } + return 0; +} diff --git a/test/doc/cyl_lbessel_j1.cpp b/test/doc/cyl_lbessel_j1.cpp new file mode 100644 index 00000000..5ec32440 --- /dev/null +++ b/test/doc/cyl_lbessel_j1.cpp @@ -0,0 +1,19 @@ +#include +#include +#include + +int main() +{ + std::cout.precision(16); + using w_t = eve::wide>; + auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); +// auto z = kyosu::complex(30.0, 0.0); +// auto zi= kyosu::complex(00.0, 30.0); + // auto z = kyosu::complex(1.0); + auto in = kyosu::cyl_lbessel_i1(z); + std::cout << "0 " << in << std::endl; + std::cout << "1 " << kyosu::cyl_bessel_i1(z) << std::endl; + auto login = kyosu::log(kyosu::cyl_bessel_i1(z)); + std::cout << "2 " << login << std::endl; + return 0; +} From c11830ad5eba942ae854186e0d62ad77dada65bb Mon Sep 17 00:00:00 2001 From: jtlap Date: Mon, 30 Oct 2023 15:21:58 +0100 Subject: [PATCH 05/17] y0 --- include/kyosu/functions.hpp | 2 + include/kyosu/functions/cyl_bessel_y0.hpp | 5 +- include/kyosu/types/impl/bessel.hpp | 1 + .../kyosu/types/impl/detail/bessel_utils.hpp | 58 +++++-------------- include/kyosu/types/impl/detail/bessely0.hpp | 36 ++++++++++++ test/doc/cyl_bessel_y0.cpp | 33 +++++++++++ 6 files changed, 89 insertions(+), 46 deletions(-) create mode 100644 include/kyosu/types/impl/detail/bessely0.hpp create mode 100644 test/doc/cyl_bessel_y0.cpp diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index fb179499..fca40d1e 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -53,6 +53,8 @@ #include #include #include +#include + #include #include #include diff --git a/include/kyosu/functions/cyl_bessel_y0.hpp b/include/kyosu/functions/cyl_bessel_y0.hpp index b597c144..d59f5481 100644 --- a/include/kyosu/functions/cyl_bessel_y0.hpp +++ b/include/kyosu/functions/cyl_bessel_y0.hpp @@ -19,7 +19,10 @@ namespace kyosu::tags KYOSU_DEFERS_CALLABLE(cyl_bessel_y0_); template - static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::cyl_bessel_y0(v); } + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept + { + return eve::cyl_bessel_y0(v); + } template KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp index a6f54575..79c0fafe 100644 --- a/include/kyosu/types/impl/bessel.hpp +++ b/include/kyosu/types/impl/bessel.hpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include diff --git a/include/kyosu/types/impl/detail/bessel_utils.hpp b/include/kyosu/types/impl/detail/bessel_utils.hpp index 3bf7dc9a..4c069178 100644 --- a/include/kyosu/types/impl/detail/bessel_utils.hpp +++ b/include/kyosu/types/impl/detail/bessel_utils.hpp @@ -161,29 +161,22 @@ namespace kyosu::_ return r; } - template - inline auto bjn(size_t n, Z z) noexcept - // compute log(Jn(z)) + template + inline auto R(size_t n, size_t p, Z z) noexcept + // compute the ratio Jn(n, z)/Jp(n-1, z)J { - Z r{}; - if (n > 0) + using u_t = eve::underlying_type_t; + if (is_eqz(n)) return Z(eve::nan(as())); + if (n == p) return Z(eve::one(as())); + Z r{}, s{}; + auto rz = rec(z); + auto rn = R(n, z); + for(int i=n-1; i > p; --i) { - auto rz = rec(z); - auto rn = R(n, z); - auto pri = rn; - for(int i=n-1; i > 0; --i) - { - rn = 2*(i)*rz-rec(rn); -// std::cout << "i " << i << " ri " << rn << std::endl; -// std::cout << "i " << i << " Ri " << R(i,z) << std::endl; - pri *= rn; - } - auto s = kyosu::cyl_bessel_j0(z)/pri; - return s; + rn = 2*i*rz-rec(rn); + pri *= rn; } - else - return cyl_bessel_j0(z); - + return lri; } template @@ -198,29 +191,4 @@ namespace kyosu::_ else return z; } - template - inline auto lbesseljn(size_t n, Z z) noexcept - // compute log(Jn(z)) - { - Z r{}, s{}; - if (n > 0) - { - auto rz = rec(z); - auto rn = R(n, z); - auto lri = log(rn); - for(int i=n-1; i > 0; --i) - { - rn = 2*(i)*rz-rec(rn); - std::cout << "i " << i << " ri " << rn << std::endl; - std::cout << "i " << i << " Ri " << R(i,z) << std::endl; - lri += log(rn); - } - s = log(kyosu::cyl_bessel_j0(z))-lri; - } - else - s = kyosu::log(cyl_bessel_j0(z)); - // ipart(s) = eve::to_nearest(eve::rem)(ipart(s), 2*eve::pi(eve::as(ipart(s)))); - return arg_adjust(s); - } - } diff --git a/include/kyosu/types/impl/detail/bessely0.hpp b/include/kyosu/types/impl/detail/bessely0.hpp new file mode 100644 index 00000000..a8443e38 --- /dev/null +++ b/include/kyosu/types/impl/detail/bessely0.hpp @@ -0,0 +1,36 @@ +//====================================================================================================================== +/* + 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); + Z s{}, sk{}; + auto sgn = -rec(j0z); + int k = 1; + do { + sk = sgn*cyl_bessel_jn(k+k, z)/k; + ++k; + sgn = -sgn; + s+= sk; + } while (eve::any(kyosu::abs(sk) > abs(s)*eps)); + return twoopi*((log(z/2)+egamma)-2*s)*cyl_bessel_j0(z); + } +} diff --git a/test/doc/cyl_bessel_y0.cpp b/test/doc/cyl_bessel_y0.cpp new file mode 100644 index 00000000..b4350f33 --- /dev/null +++ b/test/doc/cyl_bessel_y0.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::cyl_bessel_y0; + std::cout.precision(16); + std::cout << std::scientific; + std::cout << "1.0 " << " -> " << cyl_bessel_y0(1.0) << "\n"; + std::cout << "1+0i " << " -> " << cyl_bessel_y0(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << cyl_bessel_y0(15.0) << "\n"; + std::cout << "15+0i " << " -> " << cyl_bessel_y0(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << cyl_bessel_y0(40.0) << "\n"; + std::cout << "40+0i " << " -> " << cyl_bessel_y0(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << cyl_bessel_y0(60.0) << "\n"; + std::cout << "60+0i " << " -> " << cyl_bessel_y0(kyosu::complex(60.0, 0.0)) << "\n"; + + eve::wide> z(1.0, 15.0, 40.0, 60.0); + auto zz = kyosu::complex(z); + std::cout << z << " -> \n" << cyl_bessel_y0(z) << "\n"; + std::cout << zz << " -> \n" << cyl_bessel_y0(zz) << "\n"; + + std::cout << "1.0 " << " -> " << cyl_bessel_y0(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << cyl_bessel_y0(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_bessel_y0(kyosu::quaternion(1.0, 36.0, 2.0, 1.5)) << "\n"; + + + eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); + auto z0 = kyosu::complex(z1, z2); + std::cout << z0 << "\n -> " << cyl_bessel_y0(z0) << "\n"; + return 0; +} From 2e08de4fb0e6dcecde4f0eb88c4d0f6c6d2d3059 Mon Sep 17 00:00:00 2001 From: jtlap Date: Tue, 31 Oct 2023 15:17:23 +0100 Subject: [PATCH 06/17] y1 yn --- include/kyosu/functions.hpp | 2 + include/kyosu/functions/cyl_bessel_y1.hpp | 79 +++++++++++++++++++ include/kyosu/functions/cyl_bessel_yn.hpp | 79 +++++++++++++++++++ include/kyosu/functions/if_else.hpp | 5 ++ include/kyosu/types/impl/bessel.hpp | 2 + .../kyosu/types/impl/detail/bessel_utils.hpp | 41 +++++++--- include/kyosu/types/impl/detail/bessely0.hpp | 4 +- include/kyosu/types/impl/detail/bessely1.hpp | 29 +++++++ include/kyosu/types/impl/detail/besselyn.hpp | 37 +++++++++ .../kyosu/types/impl/detail/cyl_bessely1.hpp | 28 +++++++ test/doc/cyl_bessel_y1.cpp | 33 ++++++++ test/doc/cyl_bessel_yn.cpp | 34 ++++++++ 12 files changed, 361 insertions(+), 12 deletions(-) create mode 100644 include/kyosu/functions/cyl_bessel_y1.hpp create mode 100644 include/kyosu/functions/cyl_bessel_yn.hpp create mode 100644 include/kyosu/types/impl/detail/bessely1.hpp create mode 100644 include/kyosu/types/impl/detail/besselyn.hpp create mode 100644 include/kyosu/types/impl/detail/cyl_bessely1.hpp create mode 100644 test/doc/cyl_bessel_y1.cpp create mode 100644 test/doc/cyl_bessel_yn.cpp diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index fca40d1e..603829cb 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -54,6 +54,8 @@ #include #include #include +#include +#include #include #include diff --git a/include/kyosu/functions/cyl_bessel_y1.hpp b/include/kyosu/functions/cyl_bessel_y1.hpp new file mode 100644 index 00000000..3af150ab --- /dev/null +++ b/include/kyosu/functions/cyl_bessel_y1.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_cyl_bessel_y1: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_y1; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_y1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept + { + return eve::cyl_bessel_y1(v); + } + + template + KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) + { + return eve::tag_invoke(*this, target); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_bessel_y1 +//! @brief Computes the Bessel function of the second kind, +//! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) +//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. +//! +//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_y1(T z) noexcept; +//! template constexpr T cyl_bessel_y1(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return the cylindrical \f$J_0(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_y1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_y1 cyl_bessel_y1 = {}; +} diff --git a/include/kyosu/functions/cyl_bessel_yn.hpp b/include/kyosu/functions/cyl_bessel_yn.hpp new file mode 100644 index 00000000..703ca4c5 --- /dev/null +++ b/include/kyosu/functions/cyl_bessel_yn.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_cyl_bessel_yn: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_yn; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_yn_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& v) noexcept { return eve::cyl_bessel_yn(n, v); } + + template + KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target0, target1)) + { + return eve::tag_invoke(*this, target0, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @brief Computes the Bessel functions of the first kind, +//! \f$ J_{n}(x)=\sum_{p=0}^{\infty}{\frac{(-1)^p}{p!\,\Gamma (p+n +1)}} +//! {\left({x \over 2}\right)}^{2p+n }\f$. +//! +//! In the real field, it is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which +//! \f$ y(0) = 0\f$ if \f$n \ne 0\f$ else \f$1\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_yn(int n, T z) noexcept; +//! template constexpr T cyl_bessel_yn(N n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return the cylindrical \f$J_n(z)\f$. +//! +//! @warning Up to now the cayley_dickson versions have only been imlemented dor scalar int values of n. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_yn.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_yn cyl_bessel_yn = {}; +} diff --git a/include/kyosu/functions/if_else.hpp b/include/kyosu/functions/if_else.hpp index 3361cdda..7708eb05 100644 --- a/include/kyosu/functions/if_else.hpp +++ b/include/kyosu/functions/if_else.hpp @@ -15,6 +15,11 @@ namespace kyosu::tags { using callable_tag_type = callable_if_else; + KYOSU_DEFERS_CALLABLE(if_else_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, auto c, T const& t, U u) noexcept { return eve::if_else(c, t, u); } + KYOSU_FORCEINLINE auto operator()(auto const& m, auto const& t, auto const& f) const noexcept -> decltype(eve::tag_invoke(*this,m,t,f)) { diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp index 79c0fafe..40caa155 100644 --- a/include/kyosu/types/impl/bessel.hpp +++ b/include/kyosu/types/impl/bessel.hpp @@ -20,6 +20,8 @@ #include #include #include +#include +#include #include #include #include diff --git a/include/kyosu/types/impl/detail/bessel_utils.hpp b/include/kyosu/types/impl/detail/bessel_utils.hpp index 4c069178..c746c45f 100644 --- a/include/kyosu/types/impl/detail/bessel_utils.hpp +++ b/include/kyosu/types/impl/detail/bessel_utils.hpp @@ -7,7 +7,7 @@ //====================================================================================================================== #pragma once #include - +#include namespace kyosu::_ { @@ -18,7 +18,7 @@ namespace kyosu::_ auto bound(Z const & z) noexcept { using u_t = eve::underlying_type_t; - return abs(z)*u_t(1.1)+2; + return int(eve::ceil(eve::maximum(abs(z)*u_t(1.1)+2))); } //===------------------------------------------------------------------------------------------- @@ -161,24 +161,43 @@ namespace kyosu::_ return r; } +// template +// inline auto R(size_t n, size_t p, Z z) noexcept +// // compute the ratio Jn(n, z)/Jp(n-1, z)J +// { +// using u_t = eve::underlying_type_t; +// if (is_eqz(n)) return Z(eve::nan(as())); +// if (n == p) return Z(eve::one(as())); +// Z r{}, s{}; +// auto rz = rec(z); +// auto rn = R(n, z); +// for(int i=n-1; i > p; --i) +// { +// rn = 2*i*rz-rec(rn); +// pri *= rn; +// } +// return lri; +// } + + template - inline auto R(size_t n, size_t p, Z z) noexcept - // compute the ratio Jn(n, z)/Jp(n-1, z)J + inline auto Rs(size_t n, Z z) noexcept + // compute the ratio Jn(i, z)/Jn(i-1, z)J for i = n:-1:1 { + std::vector rs(n); using u_t = eve::underlying_type_t; - if (is_eqz(n)) return Z(eve::nan(as())); - if (n == p) return Z(eve::one(as())); + if (is_eqz(n)) return rs; Z r{}, s{}; auto rz = rec(z); - auto rn = R(n, z); - for(int i=n-1; i > p; --i) + rs[n-1] = R(n, z); + for(int i=n-1; i >= 1; --i) { - rn = 2*i*rz-rec(rn); - pri *= rn; + rs[i-1] = 2*i*rz-rec(rs[i]); } - return lri; + return rs; } + template inline auto arg_adjust(Z z) noexcept // modify ipart of z modulo 2*pi fit in ]-pi +pi] diff --git a/include/kyosu/types/impl/detail/bessely0.hpp b/include/kyosu/types/impl/detail/bessely0.hpp index a8443e38..e1011658 100644 --- a/include/kyosu/types/impl/detail/bessely0.hpp +++ b/include/kyosu/types/impl/detail/bessely0.hpp @@ -22,6 +22,7 @@ namespace kyosu::_ auto egamma = eve::egamma(eve::as()); auto eps = eve::eps(eve::as()); auto j0z = cyl_bessel_j0(z); + auto bd = bound(z); Z s{}, sk{}; auto sgn = -rec(j0z); int k = 1; @@ -30,7 +31,8 @@ namespace kyosu::_ ++k; sgn = -sgn; s+= sk; - } while (eve::any(kyosu::abs(sk) > abs(s)*eps)); + } while (k < bd && eve::any(kyosu::abs(sk) > abs(s)*eps)); + std::cout << "k " << k << " bound(z) "<< bound(z) << std::endl; return twoopi*((log(z/2)+egamma)-2*s)*cyl_bessel_j0(z); } } diff --git a/include/kyosu/types/impl/detail/bessely1.hpp b/include/kyosu/types/impl/detail/bessely1.hpp new file mode 100644 index 00000000..bb1b574a --- /dev/null +++ b/include/kyosu/types/impl/detail/bessely1.hpp @@ -0,0 +1,29 @@ +//====================================================================================================================== +/* + 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 y0*recs1; + } +} diff --git a/include/kyosu/types/impl/detail/besselyn.hpp b/include/kyosu/types/impl/detail/besselyn.hpp new file mode 100644 index 00000000..f9909ff3 --- /dev/null +++ b/include/kyosu/types/impl/detail/besselyn.hpp @@ -0,0 +1,37 @@ +//====================================================================================================================== +/* + 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 n, Z z) noexcept + { + std::cout << "icitte " << n << std::endl; + if (n == 0) return cyl_bessel_y0(z); + using u_t = eve::underlying_type_t; + auto twoopi = eve::two_o_pi(eve::as()); + auto rs = Rs(n, z); + auto y1 = cyl_bessel_y1(z); + auto ypred = y0; + Z pr(one(eve::as())); ; + for(int i=n-1; i > 0; --i) + { + auto recsi = rec(rs[i-1])-twoopi/(cyl_bessel_jn(i-1, z)*ypred); + pr *= recsi; + } + return pr; + } +} diff --git a/include/kyosu/types/impl/detail/cyl_bessely1.hpp b/include/kyosu/types/impl/detail/cyl_bessely1.hpp new file mode 100644 index 00000000..6b5b023a --- /dev/null +++ b/include/kyosu/types/impl/detail/cyl_bessely1.hpp @@ -0,0 +1,28 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_y1 + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, Z z) noexcept + { + using u_t = eve::underlying_type_t; + auto twoopi = eve::two_o_pi(eve::as()); + auto r1 = R(1, z); + auto y0 = cyl_bessel_y0(z); + auto recs1 = rec(r1)-twoopi/(z*cyl_bessel_j0(z)*y0) + + return y0*recs1; + } +} diff --git a/test/doc/cyl_bessel_y1.cpp b/test/doc/cyl_bessel_y1.cpp new file mode 100644 index 00000000..0ca817f2 --- /dev/null +++ b/test/doc/cyl_bessel_y1.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +int main() +{ + using kyosu::cyl_bessel_y1; + std::cout.precision(16); + std::cout << std::scientific; + std::cout << "1.0 " << " -> " << cyl_bessel_y1(1.0) << "\n"; + std::cout << "1+0i " << " -> " << cyl_bessel_y1(kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << cyl_bessel_y1(15.0) << "\n"; + std::cout << "15+0i " << " -> " << cyl_bessel_y1(kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << cyl_bessel_y1(40.0) << "\n"; + std::cout << "40+0i " << " -> " << cyl_bessel_y1(kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << cyl_bessel_y1(60.0) << "\n"; + std::cout << "60+0i " << " -> " << cyl_bessel_y1(kyosu::complex(60.0, 0.0)) << "\n"; + + eve::wide> z(1.0, 15.0, 40.0, 60.0); + auto zz = kyosu::complex(z); + std::cout << z << " -> \n" << cyl_bessel_y1(z) << "\n"; + std::cout << zz << " -> \n" << cyl_bessel_y1(zz) << "\n"; + + std::cout << "1.0 " << " -> " << cyl_bessel_y1(1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << cyl_bessel_y1(kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_bessel_y1(kyosu::quaternion(1.0, 36.0, 2.0, 1.5)) << "\n"; + + + eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); + auto z0 = kyosu::complex(z1, z2); + std::cout << z0 << "\n -> " << cyl_bessel_y1(z0) << "\n"; + return 0; +} diff --git a/test/doc/cyl_bessel_yn.cpp b/test/doc/cyl_bessel_yn.cpp new file mode 100644 index 00000000..6a6863d2 --- /dev/null +++ b/test/doc/cyl_bessel_yn.cpp @@ -0,0 +1,34 @@ +#include +#include +#include + +int main() +{ + size_t n = 3; + using kyosu::cyl_bessel_jn; + std::cout.precision(16); + std::cout << std::scientific << std::endl; + std::cout << "1.0 " << " -> " << cyl_bessel_jn(n,1.0) << "\n"; + std::cout << "1+0i " << " -> " << cyl_bessel_jn(n,kyosu::complex(1.0, 0.0)) << "\n"; + std::cout << "15.0 " << " -> " << cyl_bessel_jn(n,15.0) << "\n"; + std::cout << "15+0i " << " -> " << cyl_bessel_jn(n,kyosu::complex(15.0, 0.0)) << "\n"; + std::cout << "40.0 " << " -> " << cyl_bessel_jn(n,40.0) << "\n"; + std::cout << "40+0i " << " -> " << cyl_bessel_jn(n,kyosu::complex(40.0, 0.0)) << "\n"; + std::cout << "60.0 " << " -> " << cyl_bessel_jn(n,60.0) << "\n"; + std::cout << "60+0i " << " -> " << cyl_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 -> " << cyl_bessel_jn(n,z) << "\n"; + std::cout << zz << "\n -> " << kyosu::real(cyl_bessel_jn(n,zz)) << "\n"; + + std::cout << "1.0 " << " -> " << cyl_bessel_jn(n,1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << cyl_bessel_jn(n,kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_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-> " << cyl_bessel_jn(n,z0) << "\n"; + return 0; +} From 8afd9f22dc0ede55377707f74400e2ed962147fa Mon Sep 17 00:00:00 2001 From: jtlap Date: Tue, 31 Oct 2023 21:28:00 +0100 Subject: [PATCH 07/17] yn h1n h2n --- include/kyosu/functions.hpp | 2 + include/kyosu/functions/cyl_bessel_h1n.hpp | 77 ++++++++++++++++++++ include/kyosu/functions/cyl_bessel_h2n.hpp | 77 ++++++++++++++++++++ include/kyosu/types/impl/bessel.hpp | 2 + include/kyosu/types/impl/detail/besselhn.hpp | 46 ++++++++++++ include/kyosu/types/impl/detail/bessely0.hpp | 2 +- include/kyosu/types/impl/detail/bessely1.hpp | 5 ++ include/kyosu/types/impl/detail/besselyn.hpp | 36 ++++++--- test/doc/cyl_bessel_h1n.cpp | 20 +++++ test/doc/cyl_bessel_h2n.cpp | 20 +++++ test/doc/cyl_bessel_yn.cpp | 42 +++++------ 11 files changed, 297 insertions(+), 32 deletions(-) create mode 100644 include/kyosu/functions/cyl_bessel_h1n.hpp create mode 100644 include/kyosu/functions/cyl_bessel_h2n.hpp create mode 100644 include/kyosu/types/impl/detail/besselhn.hpp create mode 100644 test/doc/cyl_bessel_h1n.cpp create mode 100644 test/doc/cyl_bessel_h2n.cpp diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index 603829cb..f2efb65e 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -56,6 +56,8 @@ #include #include #include +#include +#include #include #include diff --git a/include/kyosu/functions/cyl_bessel_h1n.hpp b/include/kyosu/functions/cyl_bessel_h1n.hpp new file mode 100644 index 00000000..963f5950 --- /dev/null +++ b/include/kyosu/functions/cyl_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_cyl_bessel_h1n: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_h1n; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_h1n_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& z) noexcept + { + using e_t = eve::element_type_t; + return complex(cyl_bessel_jn(n, z), cyl_bessel_yn(n, z)); + } + + template + KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target0, target1)) + { + return eve::tag_invoke(*this, target0, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @brief Computes the Bessel/Hankel functions of the third kind , +//! \f$ H_n^{(1)}\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_h1n(int n, T z) noexcept; +//! template constexpr auto cyl_bessel_h1n(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return \f$H_n^{(1)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_h1n.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_h1n cyl_bessel_h1n = {}; +} diff --git a/include/kyosu/functions/cyl_bessel_h2n.hpp b/include/kyosu/functions/cyl_bessel_h2n.hpp new file mode 100644 index 00000000..64680bcd --- /dev/null +++ b/include/kyosu/functions/cyl_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_cyl_bessel_h2n: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_h1n; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_h2n_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& z) noexcept + { + using e_t = eve::element_type_t; + return complex(cyl_bessel_jn(n, z), -cyl_bessel_yn(n, z)); + } + + template + KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target0, target1)) + { + return eve::tag_invoke(*this, target0, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @brief Computes the Bessel/Hankel functions of the third kind , +//! \f$ H_n^{(2)}\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_h2n(int n, T z) noexcept; +//! template constexpr auto cyl_bessel_h2n(int n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return \f$H_n^{(2)}(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_h2n.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_h2n cyl_bessel_h2n = {}; +} diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp index 40caa155..88b164ce 100644 --- a/include/kyosu/types/impl/bessel.hpp +++ b/include/kyosu/types/impl/bessel.hpp @@ -22,6 +22,8 @@ #include #include #include +#include + #include #include #include diff --git a/include/kyosu/types/impl/detail/besselhn.hpp b/include/kyosu/types/impl/detail/besselhn.hpp new file mode 100644 index 00000000..9e7b5ea2 --- /dev/null +++ b/include/kyosu/types/impl/detail/besselhn.hpp @@ -0,0 +1,46 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include + +namespace kyosu::_ +{ + + //===------------------------------------------------------------------------------------------- + // cyl_bessel_h1n + //===------------------------------------------------------------------------------------------- + template + auto dispatch(eve::tag_of, int n, Z z) noexcept + { + if constexpr(concepts::complex) + { + auto muli = [](auto z){ auto [r, i] = z; return complex(-i, r); }; + return cyl_bessel_jn(n, z)+muli(cyl_bessel_yn(n, z)); + } + else + { + return cayley_extend_rev(cyl_bessel_h1n, n, z); + } + } + + template + auto dispatch(eve::tag_of, int n, Z z) noexcept + { + if constexpr(concepts::complex) + { + auto muli = [](auto z){ auto [r, i] = z; return complex(-i, r); }; + return cyl_bessel_jn(n, z)-muli(cyl_bessel_yn(n, z)); + } + else + { + return cayley_extend_rev(cyl_bessel_h1n, n, z); + } + } +} diff --git a/include/kyosu/types/impl/detail/bessely0.hpp b/include/kyosu/types/impl/detail/bessely0.hpp index e1011658..5e23fd07 100644 --- a/include/kyosu/types/impl/detail/bessely0.hpp +++ b/include/kyosu/types/impl/detail/bessely0.hpp @@ -32,7 +32,7 @@ namespace kyosu::_ sgn = -sgn; s+= sk; } while (k < bd && eve::any(kyosu::abs(sk) > abs(s)*eps)); - std::cout << "k " << k << " bound(z) "<< bound(z) << std::endl; +// std::cout << "k " << k << " bound(z) "<< bound(z) << std::endl; return twoopi*((log(z/2)+egamma)-2*s)*cyl_bessel_j0(z); } } diff --git a/include/kyosu/types/impl/detail/bessely1.hpp b/include/kyosu/types/impl/detail/bessely1.hpp index bb1b574a..808ec187 100644 --- a/include/kyosu/types/impl/detail/bessely1.hpp +++ b/include/kyosu/types/impl/detail/bessely1.hpp @@ -18,11 +18,16 @@ namespace kyosu::_ template auto dispatch(eve::tag_of, Z z) noexcept { + std::cout << std::endl << "latte " << 1 << std::endl; 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); + std::cout << "ypred " << y0 << std::endl; + std::cout << "r " << r1 << std::endl; auto recs1 = rec(r1)-twoopi/(z*cyl_bessel_j0(z)*y0); + std::cout << "recs1 " << recs1 << std::endl; + std::cout << "recs1*y0 " << recs1*y0 << std::endl; return y0*recs1; } diff --git a/include/kyosu/types/impl/detail/besselyn.hpp b/include/kyosu/types/impl/detail/besselyn.hpp index f9909ff3..0fd7999c 100644 --- a/include/kyosu/types/impl/detail/besselyn.hpp +++ b/include/kyosu/types/impl/detail/besselyn.hpp @@ -19,19 +19,35 @@ namespace kyosu::_ template auto dispatch(eve::tag_of, int n, Z z) noexcept { - std::cout << "icitte " << n << std::endl; - if (n == 0) return cyl_bessel_y0(z); +// std::cout << std::endl << "icitte " << n << std::endl; + auto y = cyl_bessel_y0(z); + auto jold = cyl_bessel_j0(z); + if (n == 0) return y; using u_t = eve::underlying_type_t; auto twoopi = eve::two_o_pi(eve::as()); - auto rs = Rs(n, z); - auto y1 = cyl_bessel_y1(z); - auto ypred = y0; - Z pr(one(eve::as())); ; - for(int i=n-1; i > 0; --i) + + for(int i=1; i <= n ; ++i) { - auto recsi = rec(rs[i-1])-twoopi/(cyl_bessel_jn(i-1, z)*ypred); - pr *= recsi; + auto jnew = cyl_bessel_jn(i, z); + y = (jnew*y-twoopi*rec(z))/jold; + jold = jnew; } - return pr; + return y; + +// auto rs = Rs(n, z); +// auto ypred = y0; +// std::cout << "ypred " << ypred << std::endl; +// Z pr(eve::one(eve::as())); ; +// for(int i=n; i > 0; --i) +// { +// std::cout << "r " << rs[i-1] << std::endl; +// auto recsi = rec(rs[i-1])-twoopi/(z*cyl_bessel_jn(i-1, z)*ypred); +// pr *= recsi; +// ypred = pr*y0; +// std::cout << "pr " << pr << std::endl; + +// } +// std::cout << "pr*y0 " << pr*y0 << std::endl; +// return pr*y0; } } diff --git a/test/doc/cyl_bessel_h1n.cpp b/test/doc/cyl_bessel_h1n.cpp new file mode 100644 index 00000000..1abb3ddd --- /dev/null +++ b/test/doc/cyl_bessel_h1n.cpp @@ -0,0 +1,20 @@ +#include +#include +#include + +int main() +{ + std::cout.precision(16); + // using w_t = eve::wide>; + //auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); + auto z = kyosu::complex(30.0, 0.0); +// auto z = kyosu::complex(1.0); + for(int n=2; n <= 2; ++n) + { + auto h1n = kyosu::cyl_bessel_h1n(n, z); + std::cout << "n " << n << " z " << z << std::endl; + std::cout << "1 " << h1n << std::endl; + std::cout << "2 " << kyosu::cyl_bessel_h1n(n, 30.0) << std::endl; + } + return 0; +} diff --git a/test/doc/cyl_bessel_h2n.cpp b/test/doc/cyl_bessel_h2n.cpp new file mode 100644 index 00000000..62c4f254 --- /dev/null +++ b/test/doc/cyl_bessel_h2n.cpp @@ -0,0 +1,20 @@ +#include +#include +#include + +int main() +{ + std::cout.precision(16); + // using w_t = eve::wide>; + //auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); + auto z = kyosu::complex(30.0, 0.0); +// auto z = kyosu::complex(1.0); + for(int n=2; n <= 2; ++n) + { + auto h2n = kyosu::cyl_bessel_h2n(n, z); + std::cout << "n " << n << " z " << z << std::endl; + std::cout << "0 " << h2n << std::endl; + std::cout << "2 " << kyosu::cyl_bessel_h2n(n, 30.0) << std::endl; + } + return 0; +} diff --git a/test/doc/cyl_bessel_yn.cpp b/test/doc/cyl_bessel_yn.cpp index 6a6863d2..4fdcf837 100644 --- a/test/doc/cyl_bessel_yn.cpp +++ b/test/doc/cyl_bessel_yn.cpp @@ -4,31 +4,31 @@ int main() { - size_t n = 3; - using kyosu::cyl_bessel_jn; + 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_jn(n,1.0) << "\n"; - std::cout << "1+0i " << " -> " << cyl_bessel_jn(n,kyosu::complex(1.0, 0.0)) << "\n"; - std::cout << "15.0 " << " -> " << cyl_bessel_jn(n,15.0) << "\n"; - std::cout << "15+0i " << " -> " << cyl_bessel_jn(n,kyosu::complex(15.0, 0.0)) << "\n"; - std::cout << "40.0 " << " -> " << cyl_bessel_jn(n,40.0) << "\n"; - std::cout << "40+0i " << " -> " << cyl_bessel_jn(n,kyosu::complex(40.0, 0.0)) << "\n"; - std::cout << "60.0 " << " -> " << cyl_bessel_jn(n,60.0) << "\n"; - std::cout << "60+0i " << " -> " << cyl_bessel_jn(n,kyosu::complex(60.0, 0.0)) << "\n"; + 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_jn(n,z) << "\n"; - std::cout << zz << "\n -> " << kyosu::real(cyl_bessel_jn(n,zz)) << "\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_jn(n,1.0) << "\n"; - std::cout << "1.0+36.0i " << " -> " << cyl_bessel_jn(n,kyosu::complex(1.0, 36.0)) << "\n"; - std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_bessel_jn(n,kyosu::quaternion(1.0, 36.0, 2.0, 1.5)) << "\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_jn(n,z0) << "\n"; - return 0; +// eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); +// auto z0 = kyosu::complex(z1, z2); +// std::cout << z0 << " \n-> " << cyl_bessel_yn(n,z0) << "\n"; +// return 0; } From ac9b1dec160e92d800f1abe3c1654811806b3e7e Mon Sep 17 00:00:00 2001 From: jtlap Date: Wed, 1 Nov 2023 15:27:10 +0100 Subject: [PATCH 08/17] kn and more unit tests --- include/kyosu/functions.hpp | 1 + include/kyosu/functions/cyl_bessel_kn.hpp | 79 ++++++++++++++++++++ include/kyosu/types/impl/bessel.hpp | 1 + include/kyosu/types/impl/detail/besselkn.hpp | 43 +++++++++++ include/kyosu/types/impl/detail/bessely0.hpp | 3 +- include/kyosu/types/impl/detail/bessely1.hpp | 8 +- include/kyosu/types/impl/detail/besselyn.hpp | 41 ++++------ test/doc/cyl_bessel_kn.cpp | 20 +++++ test/doc/cyl_bessel_yn.cpp | 4 +- test/doc/rstest.cpp | 21 ++++++ test/unit/function/cyl_bessel_i0.cpp | 4 +- test/unit/function/cyl_bessel_i1.cpp | 4 +- test/unit/function/cyl_bessel_in.cpp | 21 ++++-- test/unit/function/cyl_bessel_j0.cpp | 4 +- test/unit/function/cyl_bessel_j1.cpp | 4 +- test/unit/function/cyl_bessel_jn.cpp | 25 ++++--- test/unit/function/cyl_bessel_y0.cpp | 60 +++++++++++++++ test/unit/function/cyl_bessel_y1.cpp | 59 +++++++++++++++ test/unit/function/cyl_bessel_yn.cpp | 63 ++++++++++++++++ test/unit/function/cyl_lbessel_j0.cpp | 4 +- test/unit/function/cyl_lbessel_j1.cpp | 2 +- test/unit/function/cyl_lbessel_jn.cpp | 4 +- 22 files changed, 408 insertions(+), 67 deletions(-) create mode 100644 include/kyosu/functions/cyl_bessel_kn.hpp create mode 100644 include/kyosu/types/impl/detail/besselkn.hpp create mode 100644 test/doc/cyl_bessel_kn.cpp create mode 100644 test/doc/rstest.cpp create mode 100644 test/unit/function/cyl_bessel_y0.cpp create mode 100644 test/unit/function/cyl_bessel_y1.cpp create mode 100644 test/unit/function/cyl_bessel_yn.cpp diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index f2efb65e..ca073505 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -58,6 +58,7 @@ #include #include #include +#include #include #include diff --git a/include/kyosu/functions/cyl_bessel_kn.hpp b/include/kyosu/functions/cyl_bessel_kn.hpp new file mode 100644 index 00000000..caa0c924 --- /dev/null +++ b/include/kyosu/functions/cyl_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_cyl_bessel_kn: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_kn; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_kn_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& v) noexcept { return eve::cyl_bessel_kn(n, v); } + + template + KYOSU_FORCEINLINE auto operator()(N const & target0, T const& target1) const noexcept + -> decltype(eve::tag_invoke(*this, target0, target1)) + { + return eve::tag_invoke(*this, target0, target1); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @brief Computes the Bessel functions of the first kind, +//! \f$ J_{n}(x)=\sum_{p=0}^{\infty}{\frac{(-1)^p}{p!\,\Gamma (p+n +1)}} +//! {\left({x \over 2}\right)}^{2p+n }\f$. +//! +//! In the real field, it is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which +//! \f$ y(0) = 0\f$ if \f$n \ne 0\f$ else \f$1\f$. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_kn(int n, T z) noexcept; +//! template constexpr T cyl_bessel_kn(N n, T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return the cylindrical \f$J_n(z)\f$. +//! +//! @warning Up to now the cayley_dickson versions have only been imlemented dor scalar int values of n. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_kn.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_kn cyl_bessel_kn = {}; +} diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp index 88b164ce..7d73a0fc 100644 --- a/include/kyosu/types/impl/bessel.hpp +++ b/include/kyosu/types/impl/bessel.hpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include diff --git a/include/kyosu/types/impl/detail/besselkn.hpp b/include/kyosu/types/impl/detail/besselkn.hpp new file mode 100644 index 00000000..5f843f9d --- /dev/null +++ b/include/kyosu/types/impl/detail/besselkn.hpp @@ -0,0 +1,43 @@ +//====================================================================================================================== +/* + 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) + { + using u_t = eve::underlying_type_t; + auto argz = arg(z); + auto piotwo = eve::pio_2(eve::as()); + auto i = complex(u_t(0), u_t(1)); + auto cpi = piotwo*i*exp_ipi(u_t(n)/2); + auto cmi = -piotwo*i*exp_ipi(-u_t(n)/2); + auto epio2 = exp_ipi(eve::half(eve::as())); + auto empio2 = exp_ipi(eve::mhalf(eve::as())); + auto r = if_else(eve::is_ltz(argz) + , cpi*cyl_bessel_h1n(n, z*epio2) + , cmi*cyl_bessel_h2n(n, z*empio2)); + + return r; + } + else + { + return cayley_extend_rev(cyl_bessel_kn, n, z); + } + } +} diff --git a/include/kyosu/types/impl/detail/bessely0.hpp b/include/kyosu/types/impl/detail/bessely0.hpp index 5e23fd07..b34d0320 100644 --- a/include/kyosu/types/impl/detail/bessely0.hpp +++ b/include/kyosu/types/impl/detail/bessely0.hpp @@ -32,7 +32,6 @@ namespace kyosu::_ sgn = -sgn; s+= sk; } while (k < bd && eve::any(kyosu::abs(sk) > abs(s)*eps)); -// std::cout << "k " << k << " bound(z) "<< bound(z) << std::endl; - return twoopi*((log(z/2)+egamma)-2*s)*cyl_bessel_j0(z); + return if_else(is_eqz(z), complex(eve::minf(eve::as())), twoopi*((log(z/2)+egamma)-2*s)*cyl_bessel_j0(z)); } } diff --git a/include/kyosu/types/impl/detail/bessely1.hpp b/include/kyosu/types/impl/detail/bessely1.hpp index 808ec187..d2705a28 100644 --- a/include/kyosu/types/impl/detail/bessely1.hpp +++ b/include/kyosu/types/impl/detail/bessely1.hpp @@ -18,17 +18,11 @@ namespace kyosu::_ template auto dispatch(eve::tag_of, Z z) noexcept { - std::cout << std::endl << "latte " << 1 << std::endl; 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); - std::cout << "ypred " << y0 << std::endl; - std::cout << "r " << r1 << std::endl; auto recs1 = rec(r1)-twoopi/(z*cyl_bessel_j0(z)*y0); - std::cout << "recs1 " << recs1 << std::endl; - std::cout << "recs1*y0 " << recs1*y0 << std::endl; - - return y0*recs1; + 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 index 0fd7999c..de3734d0 100644 --- a/include/kyosu/types/impl/detail/besselyn.hpp +++ b/include/kyosu/types/impl/detail/besselyn.hpp @@ -19,35 +19,22 @@ namespace kyosu::_ template auto dispatch(eve::tag_of, int n, Z z) noexcept { -// std::cout << std::endl << "icitte " << n << std::endl; - auto y = cyl_bessel_y0(z); - auto jold = cyl_bessel_j0(z); - if (n == 0) return y; using u_t = eve::underlying_type_t; - auto twoopi = eve::two_o_pi(eve::as()); - - for(int i=1; i <= n ; ++i) + auto y = cyl_bessel_y0(z); + if (n != 0) { - auto jnew = cyl_bessel_jn(i, z); - y = (jnew*y-twoopi*rec(z))/jold; - jold = jnew; - } - return y; - -// auto rs = Rs(n, z); -// auto ypred = y0; -// std::cout << "ypred " << ypred << std::endl; -// Z pr(eve::one(eve::as())); ; -// for(int i=n; i > 0; --i) -// { -// std::cout << "r " << rs[i-1] << std::endl; -// auto recsi = rec(rs[i-1])-twoopi/(z*cyl_bessel_jn(i-1, z)*ypred); -// pr *= recsi; -// ypred = pr*y0; -// std::cout << "pr " << pr << std::endl; + auto jold = cyl_bessel_j0(z); + using u_t = eve::underlying_type_t; + auto twoopi = eve::two_o_pi(eve::as()); -// } -// std::cout << "pr*y0 " << pr*y0 << std::endl; -// return pr*y0; + for(int i=1; i <= n ; ++i) + { + auto jnew = cyl_bessel_jn(i, z); + y = (jnew*y-twoopi*rec(z))/jold; + jold = jnew; + } + } + auto r = if_else(eve::is_gtz(real(z)) && is_real(z) && is_nan(y), complex(real(y)), y); + return if_else(is_eqz(z), complex(eve::minf(eve::as())), r); } } diff --git a/test/doc/cyl_bessel_kn.cpp b/test/doc/cyl_bessel_kn.cpp new file mode 100644 index 00000000..2a6a5193 --- /dev/null +++ b/test/doc/cyl_bessel_kn.cpp @@ -0,0 +1,20 @@ +#include +#include +#include + +int main() +{ + std::cout.precision(16); + // using w_t = eve::wide>; + //auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); + auto z = kyosu::complex(0.0, 3.0); +// auto z = kyosu::complex(1.0); + for(int n=2; n <= 10; ++n) + { + auto kn = kyosu::cyl_bessel_kn(n, z); + std::cout << "n " << n << " z " << z << std::endl; + std::cout << "1 " << kn << std::endl; + std::cout << "2 " << kyosu::cyl_bessel_kn(n, 4.0) << std::endl; + } + return 0; +} diff --git a/test/doc/cyl_bessel_yn.cpp b/test/doc/cyl_bessel_yn.cpp index 4fdcf837..13c147a4 100644 --- a/test/doc/cyl_bessel_yn.cpp +++ b/test/doc/cyl_bessel_yn.cpp @@ -8,8 +8,8 @@ int main() 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 << "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"; diff --git a/test/doc/rstest.cpp b/test/doc/rstest.cpp new file mode 100644 index 00000000..7ba83085 --- /dev/null +++ b/test/doc/rstest.cpp @@ -0,0 +1,21 @@ +#include +#include +#include +#include + +int main() +{ + std::cout.precision(16); + using w_t = eve::wide>; +// auto z = kyosu::complex(1.0, 3.5); + auto z = kyosu::complex(w_t(20.0, 1.0), w_t(0.0, 5.0)); +// auto z = kyosu::complex(1.0, 3.5); + auto rs = kyosu::_::Rs(10, z); + for(int n=1; n < 10; ++n) + { + auto rn = kyosu::_::R(n, z); + std::cout << "n " << n << " R(n, z) " << rn << " == " << rs[n-1] << std::endl; + } + + return 0; +} diff --git a/test/unit/function/cyl_bessel_i0.cpp b/test/unit/function/cyl_bessel_i0.cpp index 78ce839e..0d8778ea 100644 --- a/test/unit/function/cyl_bessel_i0.cpp +++ b/test/unit/function/cyl_bessel_i0.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i0 over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -56,7 +56,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i0 over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) diff --git a/test/unit/function/cyl_bessel_i1.cpp b/test/unit/function/cyl_bessel_i1.cpp index b5db1d5e..f0a20c3c 100644 --- a/test/unit/function/cyl_bessel_i1.cpp +++ b/test/unit/function/cyl_bessel_i1.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i1 over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -53,7 +53,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" } }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i1 over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) diff --git a/test/unit/function/cyl_bessel_in.cpp b/test/unit/function/cyl_bessel_in.cpp index a92ce65c..2fe133b2 100644 --- a/test/unit/function/cyl_bessel_in.cpp +++ b/test/unit/function/cyl_bessel_in.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_in over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -52,7 +52,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_in over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) @@ -60,17 +60,24 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" ) (T a0, T a1) { + using u_t = eve::underlying_type_t; auto c = kyosu::complex(a0, a1); auto cb= kyosu::conj(c); auto cm= -c; auto cr= kyosu::complex(a0); auto ci= kyosu::complex(T(0), a1); + auto zer = kyosu::complex(T(0), T(0)); + auto one = kyosu::complex(T(1), T(0)); + using kyosu::cyl_bessel_in; - auto inc = cyl_bessel_in(3, c); - TTS_IEEE_EQUAL(inc, -cyl_bessel_in(3, cm)); - TTS_IEEE_EQUAL(inc, kyosu::conj(cyl_bessel_in(3, cb))); + + for(int i=0; i < 10; ++i) + { + auto inc = cyl_bessel_in(i, c); + TTS_IEEE_EQUAL(inc, eve::sign_alternate(u_t(i))*cyl_bessel_in(i, cm)); + TTS_IEEE_EQUAL(inc, kyosu::conj(cyl_bessel_in(i, cb))); TTS_EXPECT(eve::all(kyosu::is_real(cr))); TTS_EXPECT(eve::all(kyosu::is_pure(ci))); - auto z = kyosu::complex(T(0), T(0)); - TTS_IEEE_EQUAL(cyl_bessel_in(3, z), z); + TTS_IEEE_EQUAL(cyl_bessel_in(i, zer), i ? zer : one); + } }; diff --git a/test/unit/function/cyl_bessel_j0.cpp b/test/unit/function/cyl_bessel_j0.cpp index ce52eef9..0302271b 100644 --- a/test/unit/function/cyl_bessel_j0.cpp +++ b/test/unit/function/cyl_bessel_j0.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j0 over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -54,7 +54,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" } }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j0 over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) diff --git a/test/unit/function/cyl_bessel_j1.cpp b/test/unit/function/cyl_bessel_j1.cpp index d3730aec..5e3ce470 100644 --- a/test/unit/function/cyl_bessel_j1.cpp +++ b/test/unit/function/cyl_bessel_j1.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j1 over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -52,7 +52,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j1 over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) diff --git a/test/unit/function/cyl_bessel_jn.cpp b/test/unit/function/cyl_bessel_jn.cpp index 02cdd7e4..946fed45 100644 --- a/test/unit/function/cyl_bessel_jn.cpp +++ b/test/unit/function/cyl_bessel_jn.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_jn over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -53,7 +53,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_jn over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) @@ -61,17 +61,24 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" ) (T a0, T a1) { + using u_t = eve::underlying_type_t; auto c = kyosu::complex(a0, a1); auto cb= kyosu::conj(c); auto cm= -c; auto cr= kyosu::complex(a0); auto ci= kyosu::complex(T(0), a1); + auto zer = kyosu::complex(T(0), T(0)); + auto one = kyosu::complex(T(1), T(0)); + using kyosu::cyl_bessel_jn; - auto jnc = cyl_bessel_jn(3, c); - TTS_IEEE_EQUAL(jnc, -cyl_bessel_jn(3, cm)) << "c " << c << "\n"; - TTS_IEEE_EQUAL(jnc, kyosu::conj(cyl_bessel_jn(3, cb))); - TTS_EXPECT(eve::all(kyosu::is_real(cr))); - TTS_EXPECT(eve::all(kyosu::is_pure(ci))); - auto z = kyosu::complex(T(0), T(0)); - TTS_IEEE_EQUAL(cyl_bessel_jn(3, z), z); + + for(int i=0; i < 10 ; ++i) + { + auto jnc = cyl_bessel_jn(i, c); + TTS_IEEE_EQUAL(jnc, eve::sign_alternate(u_t(i))*cyl_bessel_jn(i, cm)) << "c " << c << "\n"; + TTS_IEEE_EQUAL(jnc, kyosu::conj(cyl_bessel_jn(i, cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + TTS_IEEE_EQUAL(cyl_bessel_jn(i, zer), i ? zer : one); + } }; diff --git a/test/unit/function/cyl_bessel_y0.cpp b/test/unit/function/cyl_bessel_y0.cpp new file mode 100644 index 00000000..85085b41 --- /dev/null +++ b/test/unit/function/cyl_bessel_y0.cpp @@ -0,0 +1,60 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_y0 over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + using e_t = eve::element_type_t; + if constexpr(sizeof(e_t) == 8) + { + std::cout.precision(17); + std::array re{2.4955917718518910e+00, -1.6982259765731111e+00, 1.2315182267355251e+00, 2.1258824905764118e+00, -2.1711533910447676e+00, 9.1742698952666013e-01, -2.5479731910545422e-02, 1.4924891723416751e+00}; + + std::array im{-5.0148281002365624e-01, 2.0253551112861641e+00, -5.0175898965530097e-01, -4.0455814438323523e-01, -2.1923016231568171e-01, -1.3185119199506317e+00, 1.0779142409536568e+00, 1.2528910061509635e-01}; + + std::array reres{5.5469839057084858e-01, -1.8165994943959307e+00, 3.3680718591187470e-01, 5.6362426783532593e-01, 2.8645154466364819e-01, 6.0960351781866107e-01, -2.5569311956619334e-01, 3.8450178879111341e-01}; + + std::array imres{7.8480890054575703e-02, 5.8352584787301165e-01, -2.9298142302554275e-01, -1.2803370889276175e-02, -2.4399346361135993e-01, -9.9087866505382638e-01, 1.3205924969979352e+00, 5.2208242352489359e-02}; + + + for(int i=0; i < 8; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_yn(0, c), res, 1.0e-4) << i << " <- " << c << '\n'; + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_y0 over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + auto c = kyosu::complex(a0, a1); + auto cb= kyosu::conj(c); + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + using kyosu::cyl_bessel_y0; + auto y0c = cyl_bessel_y0(c); + TTS_IEEE_EQUAL(y0c, kyosu::conj(cyl_bessel_y0(cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + auto z = kyosu::complex(T(0), T(0)); + auto minf= kyosu::complex(eve::minf(eve::as())); + TTS_IEEE_EQUAL(cyl_bessel_y0(z), minf); +}; diff --git a/test/unit/function/cyl_bessel_y1.cpp b/test/unit/function/cyl_bessel_y1.cpp new file mode 100644 index 00000000..8e5d1a52 --- /dev/null +++ b/test/unit/function/cyl_bessel_y1.cpp @@ -0,0 +1,59 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_y1 over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + using e_t = eve::element_type_t; + if constexpr(sizeof(e_t) == 8) + { + std::cout.precision(17); + std::array re{6.5758573026995004e-01, -2.4649387808542378e+00, 1.9205766369481703e+00, -6.9739731404600436e-01, 4.4532863650021348e-01, 2.4757649834958393e+00, -5.3435309039155832e-01, -2.4349028761454212e+00}; + + std::array im{1.2035170114779326e+00, -3.5932055457567835e-01, 8.3662582734559288e-01, 2.2469103234867700e+00, -9.1023149008937376e-01, 2.4818838321237546e-01, 1.5847819595769823e+00, -1.3967310652706204e+00}; + + std::array reres{-8.1045682169396482e-01, 2.3820060953407891e-02, -9.3267933254955498e-02, -1.5790753753124993e+00, -7.5013324227681488e-01, 1.4423021239534972e-01, -8.5123355830924607e-01, 4.3341743806150423e-01}; + + std::array imres{6.3354487728865994e-01, 8.9036262959787016e-01, 5.0186176942328398e-01, -1.1543670775735402e+00, -5.8526040753707131e-01, 1.1142422613715078e-01, -4.3418909281427642e-01, 1.0069485325105600e+00}; + + for(int i=0; i < 8; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_yn(1, c), res, 1.0e-4) << i << " <- " << c << '\n'; + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_y1 over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + auto c = kyosu::complex(a0, a1); + auto cb= kyosu::conj(c); + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + using kyosu::cyl_bessel_y1; + auto y1c = cyl_bessel_y1(c); + TTS_IEEE_EQUAL(y1c, kyosu::conj(cyl_bessel_y1(cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + auto z = kyosu::complex(T(0), T(0)); + auto minf= kyosu::complex(eve::minf(eve::as())); + TTS_IEEE_EQUAL(cyl_bessel_y1(z), minf); +}; diff --git a/test/unit/function/cyl_bessel_yn.cpp b/test/unit/function/cyl_bessel_yn.cpp new file mode 100644 index 00000000..9a9c1d37 --- /dev/null +++ b/test/unit/function/cyl_bessel_yn.cpp @@ -0,0 +1,63 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_yn over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + if constexpr(sizeof(T) == 8) + { + std::cout.precision(17); + std::array re{2.0309144503737393e+00, -1.5992515729134549e+00, 1.6982551054108397e+00, -1.3327549270623136e+00, 2.6585368919111341e-01, 7.3163244235853131e-01, -1.5491151507087859e+00, 2.4846599242220186e+00}; + + std::array im{-1.4137922349350807e+00, -1.7357125614937789e+00, 4.0762417134967421e-02, -1.4645033242641059e+00, 1.5914290999536518e+00, 2.3057243858262124e+00, -2.2075167685839387e+00, -5.7866716195750623e-01}; + + std::array reres{-4.4062036092080337e-01, -3.3967954461649441e-01, -1.5667619074526007e+00, -4.0477676250027339e-01, 5.9271179286244735e-01, 3.7250752149022126e-01, -4.0842410865380668e-01, -6.8424644495325704e-01}; + + std::array imres{-3.6401836760736062e-01, -4.6948726732250373e-01, 8.0811726440696674e-02, -6.2890286218778024e-01, -8.1624047120865173e-01, -4.2529834100001956e-01,-5.2072912299547713e-01, -2.6673624676594038e-01}; + + for(int i=0; i < 8; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_yn(3, c), res, 1.0e-6) << i << " <- " << c << '\n'; + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_yn over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + auto z = kyosu::complex(a0, a1); + auto re = kyosu::complex(eve::abs(a0)); + auto zer = kyosu::complex(T(0), T(0)); + auto minf= kyosu::complex(eve::minf(eve::as())); + { + auto conjz = kyosu::conj(z); + + for(int i=0; i <10 ; ++i) + { + auto jnz = kyosu::cyl_bessel_yn(i, z); + TTS_IEEE_EQUAL(jnz, kyosu::conj(kyosu::cyl_bessel_yn(i, conjz))); + auto yre = kyosu::cyl_bessel_yn(i, re); + TTS_EXPECT(eve::all(kyosu::is_real(yre))); + auto yzer = kyosu::cyl_bessel_yn(i, zer); + TTS_IEEE_EQUAL(yzer, minf); + } + } +}; diff --git a/test/unit/function/cyl_lbessel_j0.cpp b/test/unit/function/cyl_lbessel_j0.cpp index 8f67b136..e4aeb996 100644 --- a/test/unit/function/cyl_lbessel_j0.cpp +++ b/test/unit/function/cyl_lbessel_j0.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_lbessel_j0 over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -51,7 +51,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_lbessel_j0 over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) diff --git a/test/unit/function/cyl_lbessel_j1.cpp b/test/unit/function/cyl_lbessel_j1.cpp index 2c71a74a..f1c552ed 100644 --- a/test/unit/function/cyl_lbessel_j1.cpp +++ b/test/unit/function/cyl_lbessel_j1.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_lbessel_j1 over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) diff --git a/test/unit/function/cyl_lbessel_jn.cpp b/test/unit/function/cyl_lbessel_jn.cpp index 1696861a..a758e85b 100644 --- a/test/unit/function/cyl_lbessel_jn.cpp +++ b/test/unit/function/cyl_lbessel_jn.cpp @@ -8,7 +8,7 @@ #include #include -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_jn over real" , kyosu::scalar_real_types , tts::generate(tts::randoms(-10,10)) ) @@ -53,7 +53,7 @@ TTS_CASE_WITH ( "Check kyosu::abs over real" }; -TTS_CASE_WITH ( "Check kyosu::abs over real" +TTS_CASE_WITH ( "Check kyosu::cyl_lbessel_jn over real" , kyosu::real_types , tts::generate(tts::randoms(-10,10), tts::randoms(-10,10) From 640d5b5a5cb6ade6851f7fd2f4d27c46022046f7 Mon Sep 17 00:00:00 2001 From: jtlap Date: Wed, 1 Nov 2023 22:16:56 +0100 Subject: [PATCH 09/17] corretion to jn more tests k0 k1 --- include/kyosu/functions.hpp | 2 + include/kyosu/functions/cyl_bessel_k0.hpp | 73 ++++++++++++++++++++ include/kyosu/functions/cyl_bessel_k1.hpp | 73 ++++++++++++++++++++ include/kyosu/types/impl/detail/besselhn.hpp | 2 +- include/kyosu/types/impl/detail/besseljn.hpp | 9 +-- include/kyosu/types/impl/detail/besselkn.hpp | 22 ++++-- include/kyosu/types/impl/detail/bessely0.hpp | 4 ++ include/kyosu/types/impl/detail/bessely1.hpp | 3 + test/unit/function/cyl_bessel_k0.cpp | 55 +++++++++++++++ test/unit/function/cyl_bessel_k1.cpp | 56 +++++++++++++++ test/unit/function/cyl_bessel_y1.cpp | 2 +- 11 files changed, 289 insertions(+), 12 deletions(-) create mode 100644 include/kyosu/functions/cyl_bessel_k0.hpp create mode 100644 include/kyosu/functions/cyl_bessel_k1.hpp create mode 100644 test/unit/function/cyl_bessel_k0.cpp create mode 100644 test/unit/function/cyl_bessel_k1.cpp diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index ca073505..c3d675f2 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -58,6 +58,8 @@ #include #include #include +#include +#include #include #include diff --git a/include/kyosu/functions/cyl_bessel_k0.hpp b/include/kyosu/functions/cyl_bessel_k0.hpp new file mode 100644 index 00000000..b1bdb2fd --- /dev/null +++ b/include/kyosu/functions/cyl_bessel_k0.hpp @@ -0,0 +1,73 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_bessel_k0: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_k0; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_k0_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::cyl_bessel_k0(v); } + + template + KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) + { + return eve::tag_invoke(*this, target); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_bessel_k0 +//! @brief Computes the Bessel function of the second kind, +//! \f$ K_0(x)\f$ extended to the complex plane and cayley_dickson values. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_k0(T z) noexcept; +//! template constexpr auto cyl_bessel_k0(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return the cylindrical \f$K_0(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_k0.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_k0 cyl_bessel_k0 = {}; +} diff --git a/include/kyosu/functions/cyl_bessel_k1.hpp b/include/kyosu/functions/cyl_bessel_k1.hpp new file mode 100644 index 00000000..e7e3609a --- /dev/null +++ b/include/kyosu/functions/cyl_bessel_k1.hpp @@ -0,0 +1,73 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once + +#include +#include + +namespace kyosu::tags +{ + struct callable_cyl_bessel_k1: eve::elementwise + { + using callable_tag_type = callable_cyl_bessel_k1; + + KYOSU_DEFERS_CALLABLE(cyl_bessel_k1_); + + template + static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept { return eve::cyl_bessel_k1(v); } + + template + KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) + { + return eve::tag_invoke(*this, target); + } + + template + eve::unsupported_call operator()(T&&... x) const + requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; + }; +} + +namespace kyosu +{ +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var cyl_bessel_k1 +//! @brief Computes the Bessel function of the second kind, +//! \f$ K_1(x)\f$ extended to the complex plane and cayley_dickson values. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! template constexpr auto cyl_bessel_k1(T z) noexcept; +//! template constexpr auto cyl_bessel_k1(T z) noexcept; +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `z`: Value to process. +//! +//! **Return value** +//! +//! * return the cylindrical \f$K_1(z)\f$. +//! +//! @groupheader{Example} +//! +//! @godbolt{doc/cyl_bessel_k1.cpp} +//! @} +//====================================================================================================================== +inline constexpr tags::callable_cyl_bessel_k1 cyl_bessel_k1 = {}; +} diff --git a/include/kyosu/types/impl/detail/besselhn.hpp b/include/kyosu/types/impl/detail/besselhn.hpp index 9e7b5ea2..95ddd8be 100644 --- a/include/kyosu/types/impl/detail/besselhn.hpp +++ b/include/kyosu/types/impl/detail/besselhn.hpp @@ -40,7 +40,7 @@ namespace kyosu::_ } else { - return cayley_extend_rev(cyl_bessel_h1n, n, z); + return cayley_extend_rev(cyl_bessel_h2n, n, z); } } } diff --git a/include/kyosu/types/impl/detail/besseljn.hpp b/include/kyosu/types/impl/detail/besseljn.hpp index 21a96158..6e6e608b 100644 --- a/include/kyosu/types/impl/detail/besseljn.hpp +++ b/include/kyosu/types/impl/detail/besseljn.hpp @@ -66,7 +66,7 @@ namespace kyosu::_ // Starting point for backward recurrence // for when |Jn(x)|~10e-mg // using the secant method. - auto n0 = inc(eve::nearest( u_t(1.1)*az)); + auto n0 = inc(eve::ceil( u_t(1.1)*az)); auto f0 = minus_log10_cyl_j_at_infinity(n0, az) - mg; auto n1 = n0 + 5; auto f1 = minus_log10_cyl_j_at_infinity(n1, az) - mg; @@ -94,7 +94,7 @@ namespace kyosu::_ auto ejn = minus_log10_cyl_j_at_infinity(n, az); auto t = ejn <= hmp; auto obj = eve::if_else(t, sd, hmp+ejn); - auto n0 = eve::if_else(t, eve::nearest(e_t(1.1)*az), n); + auto n0 = eve::if_else(t, eve::ceil(e_t(1.1)*az), n); auto f0 = minus_log10_cyl_j_at_infinity(n0, az) - obj; auto n1 = n0 + 5; auto f1 = minus_log10_cyl_j_at_infinity(n1, az) - obj; @@ -115,8 +115,9 @@ namespace kyosu::_ }; auto backward = [az, n, ini_for_br_1, ini_for_br_2](auto z){ - auto m = ini_for_br_1(az, e_t(200)); - m = eve::if_else ( m >= n, ini_for_br_2(n, az, e_t(15)), m); + auto m1 = ini_for_br_1(az, e_t(200)); + auto m2 = ini_for_br_2(n, az, e_t(15)); + auto m = eve::if_else( m1 >= n && eve::is_not_nan(m2), m2, m1); auto cf2 = Z(0); auto cf1 = complex(eve::sqrtsmallestposval(eve::as< e_t>())); Z cf(cf2); diff --git a/include/kyosu/types/impl/detail/besselkn.hpp b/include/kyosu/types/impl/detail/besselkn.hpp index 5f843f9d..3e563519 100644 --- a/include/kyosu/types/impl/detail/besselkn.hpp +++ b/include/kyosu/types/impl/detail/besselkn.hpp @@ -6,15 +6,13 @@ */ //====================================================================================================================== #pragma once -#include -#include -#include +#include namespace kyosu::_ { //===------------------------------------------------------------------------------------------- - // cyl_bessel_h1n + // cyl_bessel_kn //===------------------------------------------------------------------------------------------- template auto dispatch(eve::tag_of, int n, Z z) noexcept @@ -32,12 +30,24 @@ namespace kyosu::_ auto r = if_else(eve::is_ltz(argz) , cpi*cyl_bessel_h1n(n, z*epio2) , cmi*cyl_bessel_h2n(n, z*empio2)); - - return r; +// std::cout << "zut " << cyl_bessel_h2n(n, z*empio2) << " === " << z*empio2 << std::endl; + return if_else(is_eqz(z), complex(eve::one(eve::as())), r); } else { return cayley_extend_rev(cyl_bessel_kn, n, z); } } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return cyl_bessel_kn(0, z); + } + + template + auto dispatch(eve::tag_of, Z z) noexcept + { + return cyl_bessel_kn(1, z); + } } diff --git a/include/kyosu/types/impl/detail/bessely0.hpp b/include/kyosu/types/impl/detail/bessely0.hpp index b34d0320..9d9b291b 100644 --- a/include/kyosu/types/impl/detail/bessely0.hpp +++ b/include/kyosu/types/impl/detail/bessely0.hpp @@ -22,12 +22,16 @@ namespace kyosu::_ auto egamma = eve::egamma(eve::as()); auto eps = eve::eps(eve::as()); auto j0z = cyl_bessel_j0(z); +// std::cout << "z " << z << " j0z " << j0z << std::endl; auto bd = bound(z); Z s{}, sk{}; auto sgn = -rec(j0z); +// std::cout << "sgn " << sgn << std::endl; int k = 1; do { +// std::cout << "jn(2k, z) " << cyl_bessel_jn(k+k, z) << std::endl; sk = sgn*cyl_bessel_jn(k+k, z)/k; +// std::cout << "k " << k << " sk " << sk << std::endl; ++k; sgn = -sgn; s+= sk; diff --git a/include/kyosu/types/impl/detail/bessely1.hpp b/include/kyosu/types/impl/detail/bessely1.hpp index d2705a28..86d49380 100644 --- a/include/kyosu/types/impl/detail/bessely1.hpp +++ b/include/kyosu/types/impl/detail/bessely1.hpp @@ -21,7 +21,10 @@ namespace kyosu::_ using u_t = eve::underlying_type_t; auto twoopi = eve::two_o_pi(eve::as()); auto r1 = R(1, z); +// std::cout << "r1 " << r1 << std::endl; auto y0 = cyl_bessel_y0(z); +// std::cout << "y0 " << y0 << std::endl; +// std::cout << "j0 " << cyl_bessel_j0(z) << std::endl; 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/test/unit/function/cyl_bessel_k0.cpp b/test/unit/function/cyl_bessel_k0.cpp new file mode 100644 index 00000000..35859834 --- /dev/null +++ b/test/unit/function/cyl_bessel_k0.cpp @@ -0,0 +1,55 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k0 over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + using e_t = eve::element_type_t; + if constexpr(sizeof(e_t) == 8) + { + std::cout.precision(17); + std::array re{2.3584877349134130e+00, 9.0164385767932564e-01, -1.5393599740361523e+00, 1.5668095462351217e+00, 1.6179702389500128e+00, -5.5681137226619104e-01, -2.2990386896893691e+00, 1.9300383451287466e+00}; + + std::array im{-8.0071372624748560e-01, 1.9823595245805736e+00, -1.3978761922193099e+00, -1.4654991890321578e+00, -5.3416883666219606e-01, -1.9265301646062505e+00, -1.6132600067568537e+00, -2.3575615216282779e+00}; + + std::array reres{4.1782770293350736e-02, -2.7083150384809346e-01, -3.4769747409921696e+00, -4.0882933553504530e-02, 1.4002108258106824e-01, -1.4548439736185870e+00, -7.3832781780344741e+00, -9.4738063303879627e-02}; + + std::array imres{5.8709121187673685e-02, -1.9749705340577356e-01, 2.3664976923056558e+00, 1.6665625358420463e-01, 1.1276125697599239e-01, 5.3154381678401141e-01, 2.2763972412185578e+00, 3.6519862966414102e-02}; + for(int i=0; i < 8; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_k0(c), res, 1.0e-4) << i << " <- " << c << '\n'; + } + } +}; + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k0 over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + auto c = kyosu::complex(a0, a1); + auto cb= kyosu::conj(c); + auto cr= kyosu::complex(a0); + using kyosu::cyl_bessel_k0; + auto k0c = cyl_bessel_k0(c); + TTS_IEEE_EQUAL(k0c, kyosu::conj(cyl_bessel_k0(cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + auto z = kyosu::complex(T(0), T(0)); + auto o = kyosu::complex(T(1), T(0)); + TTS_IEEE_EQUAL(cyl_bessel_k0(z), o); +}; diff --git a/test/unit/function/cyl_bessel_k1.cpp b/test/unit/function/cyl_bessel_k1.cpp new file mode 100644 index 00000000..0196bb02 --- /dev/null +++ b/test/unit/function/cyl_bessel_k1.cpp @@ -0,0 +1,56 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k1 over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + using e_t = eve::element_type_t; + if constexpr(sizeof(e_t) == 8) + { + std::cout.precision(17); + std::array re{-7.3954299670486945e-01, -3.7928780215327196e-02, -1.4509312401953633e+00, 2.2964197111982774e+00, -7.9699549838500194e-01, 5.4099472045802843e-01, -5.8957905551605738e-01, -5.1146527014834298e-01}; + + std::array im{-2.1345822747239773e+00, -3.3484907524396035e-01, 2.2665940322875286e+00, -1.0045032564810230e+00, -2.4392625060316648e+00, 2.6547418070570761e-01, -2.3352892270792931e+00, 1.3879256025736524e+00}; + + std::array reres{-1.7027762686598518e+00, -5.8170412824877404e-01, -2.9141064941382457e+00, 2.7872334701787467e-02, -1.5779496640047279e+00, 1.1145462290717538e+00, -1.3982594001564612e+00, -1.4628928686130445e+00}; + + std::array imres{-1.9261773906655510e-01, 3.2589115393837513e+00, 1.0075935641554303e+00, 8.5373399021791460e-02, -6.4553389686343909e-01, -7.8086954815494913e-01, -3.6962699681190719e-01, -8.7620662002477645e-01}; + for(int i=0; i < 8; ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_k1(c), res, 1.1e-6) << i << " <- " << c << " -- "<< kyosu::arg(c) << '\n'; + } + } +}; + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k1 over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ + auto c = kyosu::complex(a0, a1); + auto cb= kyosu::conj(c); + auto cr= kyosu::complex(a0); + using kyosu::cyl_bessel_k1; + auto k1c = cyl_bessel_k1(c); + TTS_IEEE_EQUAL(k1c, kyosu::conj(cyl_bessel_k1(cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + auto z = kyosu::complex(T(0), T(0)); + auto o = kyosu::complex(T(1), T(0)); + TTS_IEEE_EQUAL(cyl_bessel_k1(z), o); + +}; diff --git a/test/unit/function/cyl_bessel_y1.cpp b/test/unit/function/cyl_bessel_y1.cpp index 8e5d1a52..11a37499 100644 --- a/test/unit/function/cyl_bessel_y1.cpp +++ b/test/unit/function/cyl_bessel_y1.cpp @@ -30,7 +30,7 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_y1 over real" { auto c = kyosu::complex(re[i], im[i]); auto res = kyosu::complex(reres[i], imres[i]); - TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_yn(1, c), res, 1.0e-4) << i << " <- " << c << '\n'; + TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_y1(c), res, 1.0e-4) << i << " <- " << c << '\n'; } } }; From 1b877c5cd08689c72736a8ec1fc1df5c672a3912 Mon Sep 17 00:00:00 2001 From: jtlap Date: Thu, 2 Nov 2023 10:13:52 +0100 Subject: [PATCH 10/17] doc test for j0 --- test/doc/cyl_bessel_j0.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/doc/cyl_bessel_j0.cpp b/test/doc/cyl_bessel_j0.cpp index c6fc42f9..ef0deab7 100644 --- a/test/doc/cyl_bessel_j0.cpp +++ b/test/doc/cyl_bessel_j0.cpp @@ -29,5 +29,13 @@ int main() eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); auto z0 = kyosu::complex(z1, z2); std::cout << z0 << "\n -> " << cyl_bessel_j0(z0) << "\n"; + + eve::wide> l(1.0, 15.0, 40.0, 60.0, eve::inf(eve::as()),eve::minf(eve::as()), eve::nan(eve::as()), 0.0); + auto lims = kyosu::complex(l); + std::cout << lims << " -> \n" << cyl_bessel_j0(lims) << "\n"; + std::cout << l << " -> \n" << cyl_bessel_j0(l) << "\n"; + auto nanan = kyosu::complex(eve::nan(eve::as()), eve::nan(eve::as())); + std::cout << nanan << " -> \n" << cyl_bessel_j0(nanan) << std::endl; + return 0; } From c035f4745eaf2230c8db49057031fff03e725dc7 Mon Sep 17 00:00:00 2001 From: jtlap Date: Thu, 2 Nov 2023 16:27:53 +0100 Subject: [PATCH 11/17] test for kn and accelerating yn --- .../kyosu/types/impl/detail/bessel_utils.hpp | 47 +++++--------- include/kyosu/types/impl/detail/besselin.hpp | 4 +- include/kyosu/types/impl/detail/besseljn.hpp | 6 +- include/kyosu/types/impl/detail/besselkn.hpp | 4 +- include/kyosu/types/impl/detail/bessely0.hpp | 9 +-- include/kyosu/types/impl/detail/besselyn.hpp | 18 ++--- test/doc/rstest.cpp | 37 ++++++++--- test/unit/function/cyl_bessel_in.cpp | 1 + test/unit/function/cyl_bessel_jn.cpp | 1 + test/unit/function/cyl_bessel_kn.cpp | 65 +++++++++++++++++++ test/unit/function/cyl_bessel_yn.cpp | 10 +-- 11 files changed, 138 insertions(+), 64 deletions(-) create mode 100644 test/unit/function/cyl_bessel_kn.cpp diff --git a/include/kyosu/types/impl/detail/bessel_utils.hpp b/include/kyosu/types/impl/detail/bessel_utils.hpp index c746c45f..6798bdef 100644 --- a/include/kyosu/types/impl/detail/bessel_utils.hpp +++ b/include/kyosu/types/impl/detail/bessel_utils.hpp @@ -85,16 +85,8 @@ namespace kyosu::_ result_type f, C, delta; f = traits::b(v); - // if constexpr(b_only) -// { -// f = kyosu::if_else(kyosu::is_eqz(v), tiny, v) ; -// } -// else -// { -// f = kyosu::if_else(kyosu::is_eqz(kumi::get<0>(v)), tiny, kumi::get<0>(v)) ; -// } C = f; - result_type D{}; //(u_t(0)); + result_type D{}; std::uintmax_t counter(max_terms); do{ @@ -161,28 +153,9 @@ namespace kyosu::_ return r; } -// template -// inline auto R(size_t n, size_t p, Z z) noexcept -// // compute the ratio Jn(n, z)/Jp(n-1, z)J -// { -// using u_t = eve::underlying_type_t; -// if (is_eqz(n)) return Z(eve::nan(as())); -// if (n == p) return Z(eve::one(as())); -// Z r{}, s{}; -// auto rz = rec(z); -// auto rn = R(n, z); -// for(int i=n-1; i > p; --i) -// { -// rn = 2*i*rz-rec(rn); -// pri *= rn; -// } -// return lri; -// } - - template inline auto Rs(size_t n, Z z) noexcept - // compute the ratio Jn(i, z)/Jn(i-1, z)J for i = n:-1:1 + // compute the ratio Jn(i, z)/Jn(i+1, z)J for i = n-1:-1:0 { std::vector rs(n); using u_t = eve::underlying_type_t; @@ -198,6 +171,22 @@ namespace kyosu::_ } + template + inline auto Js(size_t n, Z z) noexcept + // compute Jn(i, z) for i = n:-1:0 + { + auto rs = Rs(n+1, z); + std::vector js(n+1); + js[n] = cyl_bessel_jn(n, z); + for(int i=n-1; i >= 0; --i) + { + js[i] = js[i+1]*rs[i]; + } + return js; + } + + + template inline auto arg_adjust(Z z) noexcept // modify ipart of z modulo 2*pi fit in ]-pi +pi] diff --git a/include/kyosu/types/impl/detail/besselin.hpp b/include/kyosu/types/impl/detail/besselin.hpp index a1659e8b..3e33b171 100644 --- a/include/kyosu/types/impl/detail/besselin.hpp +++ b/include/kyosu/types/impl/detail/besselin.hpp @@ -26,8 +26,8 @@ namespace kyosu::_ else if (n%4 == 2) return complex(eve::mone(eve::as())); else return complex(eve::zero(eve::as()), eve::one(eve::as())); }; - - return miton(n)*cyl_bessel_jn(n,complex(-ipart(z), real(z))); + auto an = eve::abs(n); + return miton(an)*cyl_bessel_jn(an,complex(-ipart(z), real(z))); } else { diff --git a/include/kyosu/types/impl/detail/besseljn.hpp b/include/kyosu/types/impl/detail/besseljn.hpp index 6e6e608b..f976dcbf 100644 --- a/include/kyosu/types/impl/detail/besseljn.hpp +++ b/include/kyosu/types/impl/detail/besseljn.hpp @@ -37,8 +37,7 @@ namespace kyosu::_ { using e_t = as_real_t; using u_t = eve::underlying_type_t; - auto n = u_t(nn); - auto an = eve::abs(n); + auto n = u_t(eve::abs(nn)); auto az = kyosu::abs(z); auto forward = [n](auto z){ @@ -154,7 +153,8 @@ namespace kyosu::_ } } auto sgnaltern = [n](auto x){return eve::if_else(eve::is_ltz(x), eve::one, eve::sign_alternate(n));}; - return sgnaltern(srz)*sgnaltern(n)*r; + r = sgnaltern(srz)*sgnaltern(n)*r; + return nn < 0 ? r*eve::sign_alternate(u_t(nn)) : r; } } else diff --git a/include/kyosu/types/impl/detail/besselkn.hpp b/include/kyosu/types/impl/detail/besselkn.hpp index 3e563519..f9bf7b4c 100644 --- a/include/kyosu/types/impl/detail/besselkn.hpp +++ b/include/kyosu/types/impl/detail/besselkn.hpp @@ -21,6 +21,7 @@ namespace kyosu::_ { using u_t = eve::underlying_type_t; auto argz = arg(z); + n = eve::abs(n); auto piotwo = eve::pio_2(eve::as()); auto i = complex(u_t(0), u_t(1)); auto cpi = piotwo*i*exp_ipi(u_t(n)/2); @@ -30,8 +31,7 @@ namespace kyosu::_ auto r = if_else(eve::is_ltz(argz) , cpi*cyl_bessel_h1n(n, z*epio2) , cmi*cyl_bessel_h2n(n, z*empio2)); -// std::cout << "zut " << cyl_bessel_h2n(n, z*empio2) << " === " << z*empio2 << std::endl; - return if_else(is_eqz(z), complex(eve::one(eve::as())), r); + return if_else(is_eqz(z), complex(eve::inf(eve::as())), r); } else { diff --git a/include/kyosu/types/impl/detail/bessely0.hpp b/include/kyosu/types/impl/detail/bessely0.hpp index 9d9b291b..17ef8ad4 100644 --- a/include/kyosu/types/impl/detail/bessely0.hpp +++ b/include/kyosu/types/impl/detail/bessely0.hpp @@ -22,20 +22,17 @@ namespace kyosu::_ auto egamma = eve::egamma(eve::as()); auto eps = eve::eps(eve::as()); auto j0z = cyl_bessel_j0(z); -// std::cout << "z " << z << " j0z " << j0z << std::endl; auto bd = bound(z); + auto js = Js(2*bd, z); Z s{}, sk{}; auto sgn = -rec(j0z); -// std::cout << "sgn " << sgn << std::endl; int k = 1; do { -// std::cout << "jn(2k, z) " << cyl_bessel_jn(k+k, z) << std::endl; - sk = sgn*cyl_bessel_jn(k+k, z)/k; -// std::cout << "k " << k << " sk " << sk << std::endl; + 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)*cyl_bessel_j0(z)); + 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/besselyn.hpp b/include/kyosu/types/impl/detail/besselyn.hpp index de3734d0..7b484b4b 100644 --- a/include/kyosu/types/impl/detail/besselyn.hpp +++ b/include/kyosu/types/impl/detail/besselyn.hpp @@ -17,24 +17,24 @@ namespace kyosu::_ // cyl_bessel_yn //===------------------------------------------------------------------------------------------- template - auto dispatch(eve::tag_of, int n, Z z) noexcept + auto dispatch(eve::tag_of, int nn, Z z) noexcept { using u_t = eve::underlying_type_t; auto y = cyl_bessel_y0(z); - if (n != 0) + if (nn != 0) { - auto jold = cyl_bessel_j0(z); + 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) { - auto jnew = cyl_bessel_jn(i, z); - y = (jnew*y-twoopi*rec(z))/jold; - jold = jnew; + 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; } - auto r = if_else(eve::is_gtz(real(z)) && is_real(z) && is_nan(y), complex(real(y)), y); - return if_else(is_eqz(z), complex(eve::minf(eve::as())), r); + else return y; } } diff --git a/test/doc/rstest.cpp b/test/doc/rstest.cpp index 7ba83085..0e5ebb73 100644 --- a/test/doc/rstest.cpp +++ b/test/doc/rstest.cpp @@ -1,21 +1,40 @@ #include #include #include -#include +//#include int main() { std::cout.precision(16); - using w_t = eve::wide>; -// auto z = kyosu::complex(1.0, 3.5); - auto z = kyosu::complex(w_t(20.0, 1.0), w_t(0.0, 5.0)); -// auto z = kyosu::complex(1.0, 3.5); - auto rs = kyosu::_::Rs(10, z); - for(int n=1; n < 10; ++n) + size_t n = 10; + auto z = kyosu::complex(1.0, 2.0); + using c_t = decltype(z); + auto rs = kyosu::_::Rs(n, z); + std::vector < c_t> rrs(n); + for(size_t i=0; i < n; ++i) { - auto rn = kyosu::_::R(n, z); - std::cout << "n " << n << " R(n, z) " << rn << " == " << rs[n-1] << std::endl; + rrs[i] = kyosu::_::R(i+1, z); + std::cout << "rs " << rs[i] << std::endl; + std::cout << "rrs " << rrs[i] << std::endl; + std::cout << "j/jm1 " << kyosu::cyl_bessel_jn(i, z)/kyosu::cyl_bessel_jn(i+1, z)<< std::endl; } + auto js = kyosu::_::Js(n, z); + for(size_t i=1; i < n; ++i) + { + std::cout << i << " -> js " << js[i] << std::endl; + std::cout << i << " -> jn " << kyosu::cyl_bessel_jn(i, z) << std::endl; + } +// std::vector js(n+1); +// for(size_t i=0; i < n+1; ++i) +// { +// js[i] = kyosu::cyl_bessel_jn(i, z); +// } +// for(size_t i=1; i < n; ++i) +// { +// std::cout << i << " -> rs " << rs[i] << std::endl; +// std::cout << i << " -> jsp1/rs " << js[i]/rs[i-1] << std::endl; +// std::cout << i << " -> js " << js[i]/js[i-1] << std::endl; +// } return 0; } diff --git a/test/unit/function/cyl_bessel_in.cpp b/test/unit/function/cyl_bessel_in.cpp index 2fe133b2..42aa423c 100644 --- a/test/unit/function/cyl_bessel_in.cpp +++ b/test/unit/function/cyl_bessel_in.cpp @@ -79,5 +79,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_in over real" TTS_EXPECT(eve::all(kyosu::is_real(cr))); TTS_EXPECT(eve::all(kyosu::is_pure(ci))); TTS_IEEE_EQUAL(cyl_bessel_in(i, zer), i ? zer : one); + TTS_IEEE_EQUAL(inc, cyl_bessel_in(-i, c)) << i << '\n'; } }; diff --git a/test/unit/function/cyl_bessel_jn.cpp b/test/unit/function/cyl_bessel_jn.cpp index 946fed45..a0c02310 100644 --- a/test/unit/function/cyl_bessel_jn.cpp +++ b/test/unit/function/cyl_bessel_jn.cpp @@ -80,5 +80,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_jn over real" TTS_EXPECT(eve::all(kyosu::is_real(cr))); TTS_EXPECT(eve::all(kyosu::is_pure(ci))); TTS_IEEE_EQUAL(cyl_bessel_jn(i, zer), i ? zer : one); + TTS_IEEE_EQUAL(jnc, eve::sign_alternate(u_t(i))*cyl_bessel_jn(-i, c)) << i << '\n'; } }; diff --git a/test/unit/function/cyl_bessel_kn.cpp b/test/unit/function/cyl_bessel_kn.cpp new file mode 100644 index 00000000..21d7c910 --- /dev/null +++ b/test/unit/function/cyl_bessel_kn.cpp @@ -0,0 +1,65 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_in over real" + , kyosu::scalar_real_types + , tts::generate(tts::randoms(-10,10)) + ) +(T ) +{ + if constexpr(sizeof(T) == 8) + { + std::cout.precision(17); + std::array re{1.4376454063874640e+00, -1.0673204704071071e-01, -1.0828666774506268e+00, -7.9765980000519210e-01, 4.2409189866361996e-01,7.6608496996954323e-01, -2.2849830522845975e+00, 4.3813679108628711e-01}; + + std::array im{-2.5175040147125849e-01, 6.2364116715495577e-01, 1.1653810042557073e+00, 1.3466103438716175e+00, -1.8438471344336298e+00, -2.0952595273684760e+00, -7.7012218588605907e-01, 1.3705948670934771e+00}; + + std::array reres{1.6867419926930389e+00, 1.5635267649025545e+01, 1.7030997128713221e+00, 2.2880149616855170e+00, -6.6737532774450159e-01, -5.4294479427715847e-01, -1.2185670369281292e+00, -2.2366901615818047e+00}; + + std::array imres{1.1484897053630303e+00, 2.9235356194719746e+01, -4.8470591465484553e-01, 8.1442309859938877e-01, -1.6108318991995840e+00, -9.0484737739480448e-01, 7.8990385458658041e-01, 2.4137425589017183e+00}; + + for(size_t i=0; i < re.size(); ++i) + { + auto c = kyosu::complex(re[i], im[i]); + auto res = kyosu::complex(reres[i], imres[i]); + TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_kn(3, c), res, 1.0e-6) << i << " <- " << c << '\n'; + } + } +}; + + +TTS_CASE_WITH ( "Check kyosu::cyl_bessel_kn over real" + , kyosu::real_types + , tts::generate(tts::randoms(-10,10), + tts::randoms(-10,10) + ) + ) +(T a0, T a1) +{ +// using u_t = eve::underlying_type_t; + auto c = kyosu::complex(a0, a1); + auto cb= kyosu::conj(c); + auto cr= kyosu::complex(a0); + auto ci= kyosu::complex(T(0), a1); + auto zer = kyosu::complex(T(0), T(0)); + auto inf = kyosu::complex(eve::inf(eve::as())); + + using kyosu::cyl_bessel_kn; + + for(int i=0; i < 10; ++i) + { + auto knc = cyl_bessel_kn(i, c); + TTS_IEEE_EQUAL(knc, kyosu::conj(cyl_bessel_kn(i, cb))); + TTS_EXPECT(eve::all(kyosu::is_real(cr))); + TTS_EXPECT(eve::all(kyosu::is_pure(ci))); + TTS_IEEE_EQUAL(cyl_bessel_kn(i, zer), inf) << i << '\n'; + TTS_IEEE_EQUAL(knc, cyl_bessel_kn(-i, c)); + } +}; diff --git a/test/unit/function/cyl_bessel_yn.cpp b/test/unit/function/cyl_bessel_yn.cpp index 9a9c1d37..c9ad9e41 100644 --- a/test/unit/function/cyl_bessel_yn.cpp +++ b/test/unit/function/cyl_bessel_yn.cpp @@ -43,6 +43,7 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_yn over real" ) (T a0, T a1) { + using u_t = eve::underlying_type_t; auto z = kyosu::complex(a0, a1); auto re = kyosu::complex(eve::abs(a0)); auto zer = kyosu::complex(T(0), T(0)); @@ -52,12 +53,13 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_yn over real" for(int i=0; i <10 ; ++i) { - auto jnz = kyosu::cyl_bessel_yn(i, z); - TTS_IEEE_EQUAL(jnz, kyosu::conj(kyosu::cyl_bessel_yn(i, conjz))); - auto yre = kyosu::cyl_bessel_yn(i, re); - TTS_EXPECT(eve::all(kyosu::is_real(yre))); + auto ynz = kyosu::cyl_bessel_yn(i, z); + TTS_IEEE_EQUAL(ynz, kyosu::conj(kyosu::cyl_bessel_yn(i, conjz))); + auto yre = kyosu::cyl_bessel_yn(i, kyosu::abs(re)); + TTS_EXPECT(eve::all(kyosu::is_real(yre)))<< i << " -> " << yre << '\n'; auto yzer = kyosu::cyl_bessel_yn(i, zer); TTS_IEEE_EQUAL(yzer, minf); + TTS_IEEE_EQUAL(ynz, eve::sign_alternate(u_t(i))*kyosu::cyl_bessel_yn(-i, z)) << i << '\n'; } } }; From 2411a34d0da2029fe52dc0619e8a6e1fa53f2281 Mon Sep 17 00:00:00 2001 From: jtlap Date: Thu, 2 Nov 2023 16:34:27 +0100 Subject: [PATCH 12/17] suppressing lbessel temporarily --- include/kyosu/functions.hpp | 7 -- include/kyosu/functions/cyl_lbessel_i0.hpp | 80 --------------------- include/kyosu/functions/cyl_lbessel_i1.hpp | 80 --------------------- include/kyosu/functions/cyl_lbessel_in.hpp | 79 --------------------- include/kyosu/functions/cyl_lbessel_j0.hpp | 80 --------------------- include/kyosu/functions/cyl_lbessel_j1.hpp | 80 --------------------- include/kyosu/functions/cyl_lbessel_jn.hpp | 81 ---------------------- include/kyosu/types/impl/detail/li1.hpp | 29 -------- include/kyosu/types/impl/detail/lin.hpp | 38 ---------- include/kyosu/types/impl/detail/lj0.hpp | 29 -------- include/kyosu/types/impl/detail/lj1.hpp | 29 -------- include/kyosu/types/impl/detail/ljn.hpp | 47 ------------- 12 files changed, 659 deletions(-) delete mode 100644 include/kyosu/functions/cyl_lbessel_i0.hpp delete mode 100644 include/kyosu/functions/cyl_lbessel_i1.hpp delete mode 100644 include/kyosu/functions/cyl_lbessel_in.hpp delete mode 100644 include/kyosu/functions/cyl_lbessel_j0.hpp delete mode 100644 include/kyosu/functions/cyl_lbessel_j1.hpp delete mode 100644 include/kyosu/functions/cyl_lbessel_jn.hpp delete mode 100644 include/kyosu/types/impl/detail/li1.hpp delete mode 100644 include/kyosu/types/impl/detail/lin.hpp delete mode 100644 include/kyosu/types/impl/detail/lj0.hpp delete mode 100644 include/kyosu/types/impl/detail/lj1.hpp delete mode 100644 include/kyosu/types/impl/detail/ljn.hpp diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index c3d675f2..92af6a8f 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -61,13 +61,6 @@ #include #include #include - -#include -#include -#include -#include -#include -#include #include #include #include diff --git a/include/kyosu/functions/cyl_lbessel_i0.hpp b/include/kyosu/functions/cyl_lbessel_i0.hpp deleted file mode 100644 index 74fcc6ff..00000000 --- a/include/kyosu/functions/cyl_lbessel_i0.hpp +++ /dev/null @@ -1,80 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright: KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once - -#include -#include - -namespace kyosu::tags -{ - struct callable_cyl_lbessel_i0: eve::elementwise - { - using callable_tag_type = callable_cyl_lbessel_i0; - - KYOSU_DEFERS_CALLABLE(cyl_lbessel_i0_); - - template - static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept - { - callable_cyl_lbessel_i0 fn{}; - return fn(complex(v)); - } - - template - KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) - { - return eve::tag_invoke(*this, target); - } - - template - eve::unsupported_call operator()(T&&... x) const - requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; - }; -} - -namespace kyosu -{ -//====================================================================================================================== -//! @addtogroup functions -//! @{ -//! @var cyl_lbessel_i0 -//! @brief Computes the logarithm of the Bessel function of the first kind, -//! \f$ J_1(x)=\frac1{\pi }\int _{0}^{\pi}\cos(\tau-x\sin \tau )\,\mathrm {d} \tau \f$ -//! extended to the complex plane and cayley_dickson values. -//! -//! \f$ J_1(\f$ is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. -//! -//! @code -//! #include -//! @endcode -//! -//! @groupheader{Callable Signatures} -//! -//! @code -//! namespace kyosu -//! { -//! template constexpr auto cyl_lbessel_i0(T z) noexcept; -//! template constexpr T cyl_lbessel_i0(T z) noexcept; -//! } -//! @endcode -//! -//! **Parameters** -//! -//! * `z`: Value to process. -//! -//! **Return value** -//! -//! * return the cylindrical \f$J_0(z)\f$. If T is a floating point value (real), a complex value is returned. -//! -//! @groupheader{Example} -//! -//! @godbolt{doc/cyl_lbessel_i0.cpp} -//! @} -//====================================================================================================================== -inline constexpr tags::callable_cyl_lbessel_i0 cyl_lbessel_i0 = {}; -} diff --git a/include/kyosu/functions/cyl_lbessel_i1.hpp b/include/kyosu/functions/cyl_lbessel_i1.hpp deleted file mode 100644 index 636de4e9..00000000 --- a/include/kyosu/functions/cyl_lbessel_i1.hpp +++ /dev/null @@ -1,80 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright: KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once - -#include -#include - -namespace kyosu::tags -{ - struct callable_cyl_lbessel_i1: eve::elementwise - { - using callable_tag_type = callable_cyl_lbessel_i1; - - KYOSU_DEFERS_CALLABLE(cyl_lbessel_i1_); - - template - static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept - { - callable_cyl_lbessel_i1 fn{}; - return fn(complex(v)); - } - - template - KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) - { - return eve::tag_invoke(*this, target); - } - - template - eve::unsupported_call operator()(T&&... x) const - requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; - }; -} - -namespace kyosu -{ -//====================================================================================================================== -//! @addtogroup functions -//! @{ -//! @var cyl_lbessel_i1 -//! @brief Computes the logarithm of the Bessel function of the first kind, -//! \f$ J_1(x)=\frac1{\pi }\int _{0}^{\pi}\cos(\tau-x\sin \tau )\,\mathrm {d} \tau \f$ -//! extended to the complex plane and cayley_dickson values. -//! -//! \f$ J_1(\f$ is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. -//! -//! @code -//! #include -//! @endcode -//! -//! @groupheader{Callable Signatures} -//! -//! @code -//! namespace kyosu -//! { -//! template constexpr auto cyl_lbessel_i1(T z) noexcept; -//! template constexpr T cyl_lbessel_i1(T z) noexcept; -//! } -//! @endcode -//! -//! **Parameters** -//! -//! * `z`: Value to process. -//! -//! **Return value** -//! -//! * return the cylindrical \f$J_0(z)\f$. If T is a floating point value (real), a complex value is returned. -//! -//! @groupheader{Example} -//! -//! @godbolt{doc/cyl_lbessel_i1.cpp} -//! @} -//====================================================================================================================== -inline constexpr tags::callable_cyl_lbessel_i1 cyl_lbessel_i1 = {}; -} diff --git a/include/kyosu/functions/cyl_lbessel_in.hpp b/include/kyosu/functions/cyl_lbessel_in.hpp deleted file mode 100644 index 748d6741..00000000 --- a/include/kyosu/functions/cyl_lbessel_in.hpp +++ /dev/null @@ -1,79 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright: KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once - -#include -#include - -namespace kyosu::tags -{ - struct callable_cyl_lbessel_in: eve::elementwise - { - using callable_tag_type = callable_cyl_lbessel_in; - - KYOSU_DEFERS_CALLABLE(cyl_lbessel_in_); - - template - static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& v) noexcept { - callable_cyl_lbessel_in fn{}; - return 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 -//! @{ -//! @brief Computes the logarithm of the Bessel functions of the first kind, -//! \f$ I_{n}(z)\f$. -//! -//! -//! @code -//! #include -//! @endcode -//! -//! @groupheader{Callable Signatures} -//! -//! @code -//! namespace kyosu -//! { -//! template constexpr auto cyl_lbessel_in(int n, T z) noexcept; -//! template constexpr T cyl_lbessel_in(int n, T z) noexcept; -//! } -//! @endcode -//! -//! **Parameters** -//! -//! * `z`: Value to process. -//! -//! **Return value** -//! -//! * return the logarithm of \f$I_n(z)\f$. -//! -//! @warning Up to now the cayley_dickson versions have only been implemented for scalar int values of n. -//! -//! @groupheader{Example} -//! -//! @godbolt{doc/cyl_lbessel_in.cpp} -//! @} -//====================================================================================================================== -inline constexpr tags::callable_cyl_lbessel_in cyl_lbessel_in = {}; -} diff --git a/include/kyosu/functions/cyl_lbessel_j0.hpp b/include/kyosu/functions/cyl_lbessel_j0.hpp deleted file mode 100644 index 0e1727de..00000000 --- a/include/kyosu/functions/cyl_lbessel_j0.hpp +++ /dev/null @@ -1,80 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright: KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once - -#include - - -namespace kyosu::tags -{ - struct callable_cyl_lbessel_j0: eve::elementwise - { - using callable_tag_type = callable_cyl_lbessel_j0; - - KYOSU_DEFERS_CALLABLE(cyl_lbessel_j0_); - - template - static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept - { - auto fn = callable_cyl_lbessel_j0{}; - return fn(complex(v)); - } - - template - KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) - { - return eve::tag_invoke(*this, target); - } - - template - eve::unsupported_call operator()(T&&... x) const - requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; - }; -} - -namespace kyosu -{ -//====================================================================================================================== -//! @addtogroup functions -//! @{ -//! @var cyl_lbessel_j0 -//! @brief Computes the logarithm of the bessel function of the first kind, -//! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) -//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. -//! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. -//! -//! @code -//! #include -//! @endcode -//! -//! @groupheader{Callable Signatures} -//! -//! @code -//! namespace kyosu -//! { -//! template constexpr auto cyl_lbessel_j0(T z) noexcept; -//! template constexpr T cyl_lbessel_j0(T z) noexcept; -//! } -//! @endcode -//! -//! **Parameters** -//! -//! * `z`: Value to process. -//! -//! **Return value** -//! -//! * return \f$\log(J_0(z))\f$. If T is a floating point value (real), a complex value is returned. -//! -//! @groupheader{Example} -//! -//! @godbolt{doc/cyl_lbessel_j0.cpp} -//! @} -//====================================================================================================================== -inline constexpr tags::callable_cyl_lbessel_j0 cyl_lbessel_j0 = {}; -} diff --git a/include/kyosu/functions/cyl_lbessel_j1.hpp b/include/kyosu/functions/cyl_lbessel_j1.hpp deleted file mode 100644 index 88cc8569..00000000 --- a/include/kyosu/functions/cyl_lbessel_j1.hpp +++ /dev/null @@ -1,80 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright: KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once - -#include -#include - -namespace kyosu::tags -{ - struct callable_cyl_lbessel_j1: eve::elementwise - { - using callable_tag_type = callable_cyl_lbessel_j1; - - KYOSU_DEFERS_CALLABLE(cyl_lbessel_j1_); - - template - static KYOSU_FORCEINLINE auto deferred_call(auto, T const& v) noexcept - { - callable_cyl_lbessel_j1 fn{}; - return fn(complex(v)); - } - - template - KYOSU_FORCEINLINE auto operator()(T const& target) const noexcept -> decltype(eve::tag_invoke(*this, target)) - { - return eve::tag_invoke(*this, target); - } - - template - eve::unsupported_call operator()(T&&... x) const - requires(!requires { eve::tag_invoke(*this, KYOSU_FWD(x)...); }) = delete; - }; -} - -namespace kyosu -{ -//====================================================================================================================== -//! @addtogroup functions -//! @{ -//! @var cyl_lbessel_j1 -//! @brief Computes the logarithm of the Bessel function of the first kind, -//! \f$ J_1(x)=\frac1{\pi }\int _{0}^{\pi}\cos(\tau-x\sin \tau )\,\mathrm {d} \tau \f$ -//! extended to the complex plane and cayley_dickson values. -//! -//! \f$ J_1(\f$ is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. -//! -//! @code -//! #include -//! @endcode -//! -//! @groupheader{Callable Signatures} -//! -//! @code -//! namespace kyosu -//! { -//! template constexpr auto cyl_lbessel_j1(T z) noexcept; -//! template constexpr T cyl_lbessel_j1(T z) noexcept; -//! } -//! @endcode -//! -//! **Parameters** -//! -//! * `z`: Value to process. -//! -//! **Return value** -//! -//! * return the cylindrical \f$J_0(z)\f$. If T is a floating point value (real), a complex value is returned. -//! -//! @groupheader{Example} -//! -//! @godbolt{doc/cyl_lbessel_j1.cpp} -//! @} -//====================================================================================================================== -inline constexpr tags::callable_cyl_lbessel_j1 cyl_lbessel_j1 = {}; -} diff --git a/include/kyosu/functions/cyl_lbessel_jn.hpp b/include/kyosu/functions/cyl_lbessel_jn.hpp deleted file mode 100644 index 5cd51c65..00000000 --- a/include/kyosu/functions/cyl_lbessel_jn.hpp +++ /dev/null @@ -1,81 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright: KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once - -#include -#include - -namespace kyosu::tags -{ - struct callable_cyl_lbessel_jn: eve::elementwise - { - using callable_tag_type = callable_cyl_lbessel_jn; - - KYOSU_DEFERS_CALLABLE(cyl_lbessel_jn_); - - template - static KYOSU_FORCEINLINE auto deferred_call(auto, N n, T const& v) noexcept - { - callable_cyl_lbessel_jn fn{}; - return 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 -//! @{ -//! @brief Computes the logarithm of the Bessel functions of the first kind, -//! \f$ J_{n}(x)=\sum_{p=0}^{\infty}{\frac{(-1)^p}{p!\,\Gamma (p+n +1)}} -//! {\left({x \over 2}\right)}^{2p+n }\f$. -//! -//! \f$ J_{n}\f$ is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which -//! \f$ y(0) = 0\f$ if \f$n \ne 0\f$ else \f$1\f$. -//! -//! @code -//! #include -//! @endcode -//! -//! @groupheader{Callable Signatures} -//! -//! @code -//! namespace kyosu -//! { -//! template constexpr auto cyl_lbessel_jn(int n, T z) noexcept; -//! template constexpr T cyl_lbessel_jn(int n, T z) noexcept; -//! } -//! @endcode -//! -//! **Parameters** -//! -//! * `z`: Value to process. -//! -//! **Return value** -//! -//! * return the cylindrical \f$J_n(z)\f$. If T is a floating point value (real), a complex value is returned. -//! @groupheader{Example} -//! -//! @godbolt{doc/cyl_lbessel_jn.cpp} -//! @} -//====================================================================================================================== -inline constexpr tags::callable_cyl_lbessel_jn cyl_lbessel_jn = {}; -} diff --git a/include/kyosu/types/impl/detail/li1.hpp b/include/kyosu/types/impl/detail/li1.hpp deleted file mode 100644 index a8d586b6..00000000 --- a/include/kyosu/types/impl/detail/li1.hpp +++ /dev/null @@ -1,29 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include - -namespace kyosu::_ -{ - - //===------------------------------------------------------------------------------------------- - // cyl_lbessel_j1 - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, Z z) noexcept - { - if constexpr(concepts::complex ) - { - return cyl_lbessel_in(1, z); - } - else - { - return cayley_extend(cyl_bessel_i1, z); - } - } -} diff --git a/include/kyosu/types/impl/detail/lin.hpp b/include/kyosu/types/impl/detail/lin.hpp deleted file mode 100644 index 320fd544..00000000 --- a/include/kyosu/types/impl/detail/lin.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_in - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, N n, Z z) noexcept - { - if constexpr(concepts::complex ) - { - using e_t = as_real_t; - auto logmiton = [](N n){ - auto nmod4 = n%4; - if (nmod4 == 0) return complex(eve::zero(eve::as())); - else if (nmod4 == 1) return complex(eve::zero(eve::as()), -eve::pio_2(eve::as())); - else if (nmod4 == 2) return complex(eve::zero(eve::as()), eve::pi(eve::as())); - else return complex(eve::zero(eve::as()), eve::pio_2(eve::as())); - }; - - return arg_adjust(logmiton(n)+cyl_lbessel_jn(n,complex(-ipart(z), real(z)))); - } - else - { - return cayley_extend_rev(cyl_bessel_in, n, z); - } - } -} diff --git a/include/kyosu/types/impl/detail/lj0.hpp b/include/kyosu/types/impl/detail/lj0.hpp deleted file mode 100644 index bedb8c5e..00000000 --- a/include/kyosu/types/impl/detail/lj0.hpp +++ /dev/null @@ -1,29 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -//#include - -namespace kyosu::_ -{ - - //===------------------------------------------------------------------------------------------- - // cyl_lbessel_j0 - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, Z z) noexcept - { - if constexpr(concepts::complex ) - { - return log(cyl_bessel_j0(z)); - } - else - { - return cayley_extend(cyl_bessel_j0, z); - } - } -} diff --git a/include/kyosu/types/impl/detail/lj1.hpp b/include/kyosu/types/impl/detail/lj1.hpp deleted file mode 100644 index aa0db7b3..00000000 --- a/include/kyosu/types/impl/detail/lj1.hpp +++ /dev/null @@ -1,29 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include - -namespace kyosu::_ -{ - - //===------------------------------------------------------------------------------------------- - // cyl_lbessel_j1 - //===------------------------------------------------------------------------------------------- - template - auto dispatch(eve::tag_of, Z z) noexcept - { - if constexpr(concepts::complex ) - { - return cyl_lbessel_jn(1, z); - } - else - { - return cayley_extend(cyl_bessel_j1, z); - } - } -} diff --git a/include/kyosu/types/impl/detail/ljn.hpp b/include/kyosu/types/impl/detail/ljn.hpp deleted file mode 100644 index 2b8f560a..00000000 --- a/include/kyosu/types/impl/detail/ljn.hpp +++ /dev/null @@ -1,47 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#pragma once -#include -#include - -namespace kyosu::_ -{ - - //===------------------------------------------------------------------------------------------- - // cyl_lbessel_jn - //===------------------------------------------------------------------------------------------- - - template - auto dispatch(eve::tag_of, N n, Z z) noexcept - { - if constexpr(concepts::complex ) - { - Z r{}, s{}; - if (n > 0) - { - auto rz = rec(z); - auto rn = R(n, z); - auto lri = log(rn); - for(int i=n-1; i > 0; --i) - { - rn = 2*i*rz-rec(rn); - lri += log(rn); - } - s = log(kyosu::cyl_bessel_j0(z))-lri; - } - else s = kyosu::log(cyl_bessel_j0(z)); - using e_t = eve::element_type_t; - return if_else(is_nez(z), arg_adjust(s), kyosu::complex(eve::minf(eve::as()))); - } - else - { - return cayley_extend_rev(cyl_lbessel_jn, n, z); - } - } - -} From e87f82ac595031999d2e5e86ffbea11ebb0bb944 Mon Sep 17 00:00:00 2001 From: jtlap Date: Thu, 2 Nov 2023 20:58:22 +0100 Subject: [PATCH 13/17] log temporary suppression --- include/kyosu/types/impl/bessel.hpp | 7 --- test/doc/cyl_lbessel_i0.cpp | 18 ------ test/doc/cyl_lbessel_i1.cpp | 18 ------ test/doc/cyl_lbessel_in.cpp | 27 -------- test/doc/cyl_lbessel_j0.cpp | 33 ---------- test/doc/cyl_lbessel_j1.cpp | 19 ------ test/doc/cyl_lbessel_jn.cpp | 22 ------- test/unit/function/cyl_bessel_k0.cpp | 8 ++- test/unit/function/cyl_bessel_k1.cpp | 6 +- test/unit/function/cyl_lbessel_j0.cpp | 85 ------------------------- test/unit/function/cyl_lbessel_j1.cpp | 83 ------------------------ test/unit/function/cyl_lbessel_jn.cpp | 90 --------------------------- 12 files changed, 9 insertions(+), 407 deletions(-) delete mode 100644 test/doc/cyl_lbessel_i0.cpp delete mode 100644 test/doc/cyl_lbessel_i1.cpp delete mode 100644 test/doc/cyl_lbessel_in.cpp delete mode 100644 test/doc/cyl_lbessel_j0.cpp delete mode 100644 test/doc/cyl_lbessel_j1.cpp delete mode 100644 test/doc/cyl_lbessel_jn.cpp delete mode 100644 test/unit/function/cyl_lbessel_j0.cpp delete mode 100644 test/unit/function/cyl_lbessel_j1.cpp delete mode 100644 test/unit/function/cyl_lbessel_jn.cpp diff --git a/include/kyosu/types/impl/bessel.hpp b/include/kyosu/types/impl/bessel.hpp index 7d73a0fc..01cb5975 100644 --- a/include/kyosu/types/impl/bessel.hpp +++ b/include/kyosu/types/impl/bessel.hpp @@ -24,10 +24,3 @@ #include #include #include - -#include -#include -#include -#include -#include -#include diff --git a/test/doc/cyl_lbessel_i0.cpp b/test/doc/cyl_lbessel_i0.cpp deleted file mode 100644 index 44d6c29b..00000000 --- a/test/doc/cyl_lbessel_i0.cpp +++ /dev/null @@ -1,18 +0,0 @@ -#include -#include -#include - -int main() -{ - std::cout.precision(16); - using w_t = eve::wide>; - auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); -// auto z = kyosu::complex(30.0, 0.0); -// auto zi= kyosu::complex(00.0, 30.0); - // auto z = kyosu::complex(1.0); - auto in = kyosu::cyl_lbessel_i0(z); - std::cout << "0 " << in << std::endl; - auto login = kyosu::log(kyosu::cyl_bessel_i0(z)); - std::cout << "1 " << login << std::endl; - return 0; -} diff --git a/test/doc/cyl_lbessel_i1.cpp b/test/doc/cyl_lbessel_i1.cpp deleted file mode 100644 index 180ff649..00000000 --- a/test/doc/cyl_lbessel_i1.cpp +++ /dev/null @@ -1,18 +0,0 @@ -#include -#include -#include - -int main() -{ - std::cout.precision(16); - using w_t = eve::wide>; - auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); -// auto z = kyosu::complex(30.0, 0.0); -// auto zi= kyosu::complex(00.0, 30.0); - // auto z = kyosu::complex(1.0); - auto in = kyosu::cyl_lbessel_i1(z); - std::cout << "0 " << in << std::endl; - auto login = kyosu::log(kyosu::cyl_bessel_i1(z)); - std::cout << "1 " << login << std::endl; - return 0; -} diff --git a/test/doc/cyl_lbessel_in.cpp b/test/doc/cyl_lbessel_in.cpp deleted file mode 100644 index c414fb83..00000000 --- a/test/doc/cyl_lbessel_in.cpp +++ /dev/null @@ -1,27 +0,0 @@ -#include -#include -#include - -int main() -{ - std::cout.precision(16); - // using w_t = eve::wide>; - //auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); - auto z = kyosu::complex(30.0, 0.0); - auto zi= kyosu::complex(00.0, 30.0); -// auto z = kyosu::complex(1.0); - for(int n=2; n <= 2; ++n) - { - auto in = kyosu::cyl_lbessel_in(n, z); - std::cout << "n " << n << " z " << z << std::endl; - std::cout << "0 " << in << std::endl; - std::cout << "1 " << kyosu::cyl_bessel_in(n, z) << std::endl; - std::cout << "2 " << kyosu::cyl_bessel_jn(n, zi) << std::endl; - std::cout << "3 " << kyosu::cyl_bessel_in(n, 20.0) << std::endl; - auto login = kyosu::log(kyosu::cyl_bessel_in(n, z)); - std::cout << "4 " << login << std::endl; - std::cout << "5 " << kyosu::exp(in) << std::endl; - std::cout << "6 " << kyosu::exp(login) << std::endl; - } - return 0; -} diff --git a/test/doc/cyl_lbessel_j0.cpp b/test/doc/cyl_lbessel_j0.cpp deleted file mode 100644 index bc1b376d..00000000 --- a/test/doc/cyl_lbessel_j0.cpp +++ /dev/null @@ -1,33 +0,0 @@ -#include -#include -#include - -int main() -{ - using kyosu::cyl_lbessel_j0; - std::cout.precision(16); - - std::cout << "1.0 " << " -> " << cyl_lbessel_j0(1.0) << "\n"; - std::cout << "1+0i " << " -> " << cyl_lbessel_j0(kyosu::complex(1.0, 0.0)) << "\n"; - std::cout << "15.0 " << " -> " << cyl_lbessel_j0(15.0) << "\n"; - std::cout << "15+0i " << " -> " << cyl_lbessel_j0(kyosu::complex(15.0, 0.0)) << "\n"; - std::cout << "40.0 " << " -> " << cyl_lbessel_j0(40.0) << "\n"; - std::cout << "40+0i " << " -> " << cyl_lbessel_j0(kyosu::complex(40.0, 0.0)) << "\n"; - std::cout << "60.0 " << " -> " << cyl_lbessel_j0(60.0) << "\n"; - std::cout << "60+0i " << " -> " << cyl_lbessel_j0(kyosu::complex(60.0, 0.0)) << "\n"; - - eve::wide> z(1.0, 15.0, 40.0, 60.0); - auto zz = kyosu::complex(z); - std::cout << z << " -> \n" << cyl_lbessel_j0(z) << "\n"; - std::cout << zz << " -> \n" << cyl_lbessel_j0(zz) << "\n"; - - std::cout << "1.0 " << " -> " << cyl_lbessel_j0(1.0) << "\n"; - std::cout << "1.0+36.0i " << " -> " << cyl_lbessel_j0(kyosu::complex(1.0, 36.0)) << "\n"; - std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_lbessel_j0(kyosu::quaternion(1.0, 36.0, 2.0, 1.5)) << "\n"; - - - eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); - auto z0 = kyosu::complex(z1, z2); - std::cout << z0 << "\n -> " << cyl_lbessel_j0(z0) << "\n"; - return 0; -} diff --git a/test/doc/cyl_lbessel_j1.cpp b/test/doc/cyl_lbessel_j1.cpp deleted file mode 100644 index 5ec32440..00000000 --- a/test/doc/cyl_lbessel_j1.cpp +++ /dev/null @@ -1,19 +0,0 @@ -#include -#include -#include - -int main() -{ - std::cout.precision(16); - using w_t = eve::wide>; - auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); -// auto z = kyosu::complex(30.0, 0.0); -// auto zi= kyosu::complex(00.0, 30.0); - // auto z = kyosu::complex(1.0); - auto in = kyosu::cyl_lbessel_i1(z); - std::cout << "0 " << in << std::endl; - std::cout << "1 " << kyosu::cyl_bessel_i1(z) << std::endl; - auto login = kyosu::log(kyosu::cyl_bessel_i1(z)); - std::cout << "2 " << login << std::endl; - return 0; -} diff --git a/test/doc/cyl_lbessel_jn.cpp b/test/doc/cyl_lbessel_jn.cpp deleted file mode 100644 index 0205c20b..00000000 --- a/test/doc/cyl_lbessel_jn.cpp +++ /dev/null @@ -1,22 +0,0 @@ -#include -#include -#include -//#include - - -int main() -{ - std::cout.precision(16); - using w_t = eve::wide>; - auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); -// auto z = kyosu::complex(1.0); - for(int n=0; n <= 10; ++n) - { - auto jn = kyosu::cyl_lbessel_jn(n, z); - std::cout << "n " << n << " z " << z << std::endl; - std::cout << jn << std::endl; - std::cout << kyosu::log(kyosu::cyl_bessel_jn(n, z)) << std::endl; -// std::cout << kyosu::_::lbesseljn(n, z) << std::endl; - } - return 0; -} diff --git a/test/unit/function/cyl_bessel_k0.cpp b/test/unit/function/cyl_bessel_k0.cpp index 35859834..acbd42a9 100644 --- a/test/unit/function/cyl_bessel_k0.cpp +++ b/test/unit/function/cyl_bessel_k0.cpp @@ -47,9 +47,11 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k0 over real" auto cr= kyosu::complex(a0); using kyosu::cyl_bessel_k0; auto k0c = cyl_bessel_k0(c); + auto inf = kyosu::complex(eve::inf(eve::as())); + auto z = kyosu::complex(T(0), T(0)); + TTS_IEEE_EQUAL(k0c, kyosu::conj(cyl_bessel_k0(cb))); TTS_EXPECT(eve::all(kyosu::is_real(cr))); - auto z = kyosu::complex(T(0), T(0)); - auto o = kyosu::complex(T(1), T(0)); - TTS_IEEE_EQUAL(cyl_bessel_k0(z), o); + TTS_IEEE_EQUAL(cyl_bessel_k0(z), inf); + TTS_IEEE_EQUAL(k0c, cyl_bessel_k0(c)); }; diff --git a/test/unit/function/cyl_bessel_k1.cpp b/test/unit/function/cyl_bessel_k1.cpp index 0196bb02..fa2e5574 100644 --- a/test/unit/function/cyl_bessel_k1.cpp +++ b/test/unit/function/cyl_bessel_k1.cpp @@ -47,10 +47,12 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k1 over real" auto cr= kyosu::complex(a0); using kyosu::cyl_bessel_k1; auto k1c = cyl_bessel_k1(c); + auto inf = kyosu::complex(eve::inf(eve::as())); + TTS_IEEE_EQUAL(k1c, kyosu::conj(cyl_bessel_k1(cb))); TTS_EXPECT(eve::all(kyosu::is_real(cr))); auto z = kyosu::complex(T(0), T(0)); - auto o = kyosu::complex(T(1), T(0)); - TTS_IEEE_EQUAL(cyl_bessel_k1(z), o); + TTS_IEEE_EQUAL(cyl_bessel_k1(z), inf); + TTS_IEEE_EQUAL(k1c, cyl_bessel_k1(c)); }; diff --git a/test/unit/function/cyl_lbessel_j0.cpp b/test/unit/function/cyl_lbessel_j0.cpp deleted file mode 100644 index e4aeb996..00000000 --- a/test/unit/function/cyl_lbessel_j0.cpp +++ /dev/null @@ -1,85 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#include -#include - -TTS_CASE_WITH ( "Check kyosu::cyl_lbessel_j0 over real" - , kyosu::scalar_real_types - , tts::generate(tts::randoms(-10,10)) - ) -(T ) -{ - using e_t = eve::element_type_t; - if constexpr(sizeof(e_t) == 8) - { - std::cout.precision(17); - std::array re{ e_t(-9.8403264736233908e+01), e_t(8.6287324687322851e+01), e_t(2.8677183208357349e+01), e_t(7.5642931882481321e+01), e_t(9.6351793077594323e+01) - , e_t(-7.8716695947540288e+01), e_t(-5.0019457963533974e+01), e_t(5.9030075032843811e+01), e_t(-7.5596124248257235e+00), e_t(8.5810875418338028e+01)}; - - std::array im{ e_t(-1.9271037994246587e+01), e_t(7.9043790847200569e+01), e_t(-3.0599555605285602e+01), e_t(-7.1720588923445590e+01), e_t(-6.1110206702938363e+01) - , e_t(9.0167523214068439e+01), e_t(-7.0131934820406244e+01), e_t(-9.3405497312943496e+01), e_t(8.3431692318213322e+01), e_t(-4.1557438097347173e+01)}; - - std::array reres{e_t(-8.8414975251226071e+06), e_t(-3.9177291392384149e+32), e_t(-1.2012865579185671e+12), e_t(5.4246749994831845e+29), e_t(-3.8201314596721725e+23) - , e_t(-5.1783207689008547e+37), e_t(1.0487044863432896e+29), e_t(-8.2091088236413212e+38), e_t(2.4924320903438798e+34), e_t(-4.1577568667728512e+16)}; - std::array imres{e_t(2.9689225084136692e+06), e_t(6.8130876628410661e+32), e_t(-2.9179495183191563e+10), e_t(-8.8752530111744507e+28), e_t(1.2946830695756804e+25) - , e_t(9.5653949361454081e+36), e_t(6.5259614844547437e+28), e_t(1.1301644890532119e+39), e_t(7.0530420820017420e+34), e_t(-1.8907103618909208e+16)}; - for(int i=0; i < 10; ++i) - { - auto c = kyosu::complex(re[i], im[i]); - auto res = kyosu::log(kyosu::complex(reres[i], imres[i])); - TTS_RELATIVE_EQUAL(kyosu::cyl_lbessel_j0(c), res, 1.0e-4) << i << " <- " << c << '\n'; - } - } - std::cout.precision(17); - std::array re{2.175102683907268e+00, 3.007738475083643e-01, 3.571791197760752e+00, 1.916834892035037e+00, 2.309329011574615e+00, 1.116302757926106e+00, 3.454188251414261e+00, 2.504516854443353e+00}; - std::array im{ 3.008362577802520e+00, 3.733804275914350e+00, 2.988312632706299e+00, 4.411110052057489e+00, 3.055249958745497e+00, 2.508994160634028e+00, 3.100976745146938e-01, 3.854880352108083e+00}; - - std::array reres{-1.1271605877499498e+00, 8.6819450415271469e+00, -3.7342482747514563e+00, -1.9594143015933545e+00, -1.6842573645694947e+00, 1.9975730077018730e+00, -3.9366090510718693e-01, -5.2795374955460659e+00}; - std::array imres{-4.1516002202378823e+00, -2.2772457219968798e+00, -1.0734738341492034e-01, -1.5246293957752114e+01, -4.0999678002136202e+00, -2.3878228304877860e+00, -4.8664505025944993e-02, -7.2799291393660077e+00}; - for(int i=0; i < 8; ++i) - { - auto c = kyosu::complex(re[i], im[i]); - auto res = kyosu::log(kyosu::complex(reres[i], imres[i])); - double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; - TTS_RELATIVE_EQUAL(kyosu::cyl_lbessel_j0(c), res, tol) << i << " <- " << c << '\n'; - } -}; - - -TTS_CASE_WITH ( "Check kyosu::cyl_lbessel_j0 over real" - , kyosu::real_types - , tts::generate(tts::randoms(-10,10), - tts::randoms(-10,10) - ) - ) -(T a0, T a1) -{ - if constexpr(sizeof(T) == 8) - { - using e_t = eve::element_type_t; - double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; - auto c = kyosu::complex(a0, a1); - auto cb= kyosu::conj(c); - auto cm= -c; - using kyosu::cyl_lbessel_j0; - using kyosu::cyl_bessel_j0; - auto lj0c = cyl_lbessel_j0(c); - auto logj0c = kyosu::log(cyl_bessel_j0(c)); - TTS_RELATIVE_EQUAL(lj0c, logj0c, tol) << "c " << c << "\n"; - - auto lj0cb = cyl_lbessel_j0(cb); - auto logj0cb = kyosu::log(cyl_bessel_j0(cb)); - TTS_RELATIVE_EQUAL(lj0cb, logj0cb, tol) << "c " << c << "\n"; - - auto lj0cm = cyl_lbessel_j0(cm); - auto logj0cm = kyosu::log(cyl_bessel_j0(cm)); - TTS_RELATIVE_EQUAL(lj0cm, logj0cm, tol) << "c " << c << "\n"; - auto z = kyosu::complex(T(0), T(0)); - TTS_IEEE_EQUAL(cyl_lbessel_j0(z), kyosu::complex(eve::zero(eve::as()))); - } -}; diff --git a/test/unit/function/cyl_lbessel_j1.cpp b/test/unit/function/cyl_lbessel_j1.cpp deleted file mode 100644 index f1c552ed..00000000 --- a/test/unit/function/cyl_lbessel_j1.cpp +++ /dev/null @@ -1,83 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#include -#include - -TTS_CASE_WITH ( "Check kyosu::cyl_lbessel_j1 over real" - , kyosu::scalar_real_types - , tts::generate(tts::randoms(-10,10)) - ) -(T ) -{ - using e_t = eve::element_type_t; - if constexpr(sizeof(e_t) == 8) - { - std::cout.precision(17); - std::array re{ -9.8403264736233908e+01, 8.6287324687322851e+01, 2.8677183208357349e+01, 7.5642931882481321e+01, 9.6351793077594323e+01 - , -7.8716695947540288e+01, -5.0019457963533974e+01, 5.9030075032843811e+01, -7.5596124248257235e+00, 8.5810875418338028e+01}; - std::array im{-1.9271037994246587e+01, 7.9043790847200569e+01, -3.0599555605285602e+01, -7.1720588923445590e+01, -6.1110206702938363e+01 - , 9.0167523214068439e+01, -7.0131934820406244e+01, -9.3405497312943496e+01, 8.3431692318213322e+01, -4.1557438097347173e+01}; - - std::array reres{3.0094177278540167e+06, -6.8058089697014938e+32, -3.8806283458664719e+10, -8.6565181213118240e+28, 1.2915079982169630e+25, - -9.3925711342553668e+36, 6.4595726071319036e+28, 1.1238442153777222e+39, -7.0123433800633681e+34, -1.9060735303763808e+16}; - std::array imres{8.8185850069141667e+06, -3.8848970122848989e+32, 1.1905919032156870e+12, -5.4098751798106477e+29, 4.2914274806719309e+23, - -5.1646555066948327e+37, -1.0459530606430648e+29, 8.2050979033727447e+38, 2.4737508957685513e+34, 4.1393451863209432e+16}; - for(int i=0; i < 10; ++i) - { - auto c = kyosu::complex(re[i], im[i]); - auto res = kyosu::log(kyosu::complex(reres[i], imres[i])); - TTS_RELATIVE_EQUAL(kyosu::cyl_lbessel_j1(c), res, 1.0e-7) << i << " <- " << c << '\n'; - } - } - std::cout.precision(17); - std::array re{2.175102683907268e+00, 3.007738475083643e-01, 3.571791197760752e+00, 1.916834892035037e+00, 2.309329011574615e+00, 1.116302757926106e+00, 3.454188251414261e+00, 2.504516854443353e+00}; - std::array im{ 3.008362577802520e+00, 3.733804275914350e+00, 2.988312632706299e+00, 4.411110052057489e+00, 3.055249958745497e+00, 2.508994160634028e+00, 3.100976745146938e-01, 3.854880352108083e+00}; - - std::array reres{3.6097545599299043e+00, 2.0607036224436954e+00, -2.1724057594208968e-01, 1.3644664632775582e+01, 3.5380568884368788e+00, 2.1507545829381751e+00, 1.5764037998430128e-01, 6.2619292962097814e+00}; - std::array imres{-1.3891357467781533e+00, 7.3810832876707710e+00, -3.4818290031753225e+00, -2.4845686655182946e+00, -1.8826359832772270e+00, 1.3851339474357725e+00, -1.3157891968242472e-01, -5.2738074786282327e+00}; - for(int i=0; i < 8; ++i) - { - auto c = kyosu::complex(re[i], im[i]); - auto res = kyosu::log(kyosu::complex(reres[i], imres[i])); - double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; - TTS_RELATIVE_EQUAL(kyosu::cyl_lbessel_j1(c), res, tol) << i << " <- " << c << '\n'; - } -}; - -// TTS_CASE_WITH ( "Check kyosu::cyl_lbessel_j1 over real" -// , kyosu::real_types -// , tts::generate(tts::randoms(-10,10), -// tts::randoms(-10,10) -// ) -// ) -// (T a0, T a1) -// { -// if constexpr(sizeof(T) == 8) -// { -// using e_t = eve::element_type_t; -// double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; -// auto c = kyosu::complex(a0, a1); -// auto cb= kyosu::conj(c); -// auto cm= -c; -// using kyosu::cyl_lbessel_j1; -// using kyosu::cyl_bessel_j1; -// auto lj1c = cyl_lbessel_j1(c); -// auto logj1c = kyosu::log(cyl_bessel_j1(c)); -// TTS_RELATIVE_EQUAL(lj1c, logj1c, tol) << "c " << c << "\n"; - -// auto lj1cb = cyl_lbessel_j1(cb); -// auto logj1cb = kyosu::log(cyl_bessel_j1(cb)); -// TTS_RELATIVE_EQUAL(lj1cb, logj1cb, tol) << "c " << c << "\n"; - -// auto lj1cm = cyl_lbessel_j1(cm); -// auto logj1cm = kyosu::log(cyl_bessel_j1(cm)); -// TTS_RELATIVE_EQUAL(lj1cm, logj1cm, tol) << "c " << c << "\n"; -// auto z = kyosu::complex(T(0), T(0)); -// TTS_IEEE_EQUAL(cyl_lbessel_j1(z), kyosu::complex(eve::minf(eve::as()))); -// } -// }; diff --git a/test/unit/function/cyl_lbessel_jn.cpp b/test/unit/function/cyl_lbessel_jn.cpp deleted file mode 100644 index a758e85b..00000000 --- a/test/unit/function/cyl_lbessel_jn.cpp +++ /dev/null @@ -1,90 +0,0 @@ -//====================================================================================================================== -/* - Kyosu - Complex Without Complexes - Copyright : KYOSU Contributors & Maintainers - SPDX-License-Identifier: BSL-1.0 -*/ -//====================================================================================================================== -#include -#include - -TTS_CASE_WITH ( "Check kyosu::cyl_bessel_jn over real" - , kyosu::scalar_real_types - , tts::generate(tts::randoms(-10,10)) - ) -(T ) -{ - if constexpr(sizeof(T) == 8) - { - std::cout.precision(17); - std::array re{ -9.8403264736233908e+01, 8.6287324687322851e+01, 2.8677183208357349e+01, 7.5642931882481321e+01, 9.6351793077594323e+01 - , -7.8716695947540288e+01, -5.0019457963533974e+01, 5.9030075032843811e+01, -7.5596124248257235e+00, 8.5810875418338028e+01}; - std::array im{-1.9271037994246587e+01, 7.9043790847200569e+01, -3.0599555605285602e+01, -7.1720588923445590e+01, -6.1110206702938363e+01 - , 9.0167523214068439e+01, -7.0131934820406244e+01, -9.3405497312943496e+01, 8.3431692318213322e+01, -4.1557438097347173e+01}; - - std::array reres{-3.3279158735087640e+06, 6.7446373092398792e+32, 1.0973422267168294e+11, 6.9510412682541761e+28, -1.2657523599731811e+25, - 8.0429582124849094e+36, -5.9430267579521807e+28, -1.0741846420111567e+39, 6.6950456879149885e+34, 2.0245919743884628e+16}; - std::array imres{-8.6289796115820184e+06, 3.6264735895350172e+32, -1.1056118662675706e+12, 5.2905159458975805e+29, -7.9797425680943093e+23, - 5.0551834834322169e+37, 1.0236118318065309e+29, -8.1681100953492296e+38, -2.3290495380992445e+34, -3.9909835121233848e+16}; - for(int i=0; i < 10; ++i) - { - auto c = kyosu::complex(re[i], im[i]); - auto res = kyosu::complex(reres[i], imres[i]); - TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_jn(3, c), res, 1.0e-7) << i << " <- " << c << '\n'; - } - } - - std::array re{-2.2013892240657134e+00, -2.8382367527759522e-01, -7.0737598522551182e-01, 8.6178612999047166e-02, 1.8217799341954277e-01, -7.6718326053284047e-01, -5.9167207480492257e-01, 7.5911892713725948e-01}; - - std::array im{8.7228967955690240e-01, -2.1140894389723703e+00, 2.2905192466959483e+00, -1.5275122978034243e+00, 5.5435596674685006e-01, -1.1418929858153533e+00, 1.1666506779144503e+00, -1.9198575359478032e+00}; - - std::array reres{-1.3672173596226031e-01, 1.2004993008474076e-01, 3.4070032613198098e-01, -1.5845109266378627e-02, -3.4617874091075179e-03, 5.4024365099246512e-02, 4.9541047773116346e-02, -2.1467277569289250e-01}; - - std::array imres{1.6497274394928027e-01, 2.3540711137379494e-01, -1.7769121638069371e-01, 8.4636178354763933e-02, -2.3973187185149024e-03, -1.7452629167850800e-02, -3.8485807736362864e-03, 5.8435910432745555e-02}; - using e_t = eve::element_type_t; - - for(size_t i=0; i < re.size(); ++i) - { - auto c = kyosu::complex(re[i], im[i]); - auto res = kyosu::log(kyosu::complex(reres[i], imres[i])); - double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; - TTS_RELATIVE_EQUAL(kyosu::cyl_lbessel_jn(3, c), res, tol) << i << " <- " << c << '\n'; - } -}; - - -TTS_CASE_WITH ( "Check kyosu::cyl_lbessel_jn over real" - , kyosu::real_types - , tts::generate(tts::randoms(-10,10), - tts::randoms(-10,10) - ) - ) -(T a0, T a1) -{ - if constexpr(sizeof(T) == 8) - { - using e_t = eve::element_type_t; - double tol = (sizeof(e_t) == 4) ? 1.e-3 :1.e-7; - auto c = kyosu::complex(a0, a1); - auto cb= kyosu::conj(c); - auto cm= -c; - auto cr = kyosu::complex(eve::abs(a0)); - using kyosu::cyl_lbessel_jn; - using kyosu::cyl_bessel_jn; - auto ljnc = cyl_lbessel_jn(3, c); - auto logjnc = kyosu::log(cyl_bessel_jn(3, c)); - TTS_RELATIVE_EQUAL(ljnc, logjnc, tol) << "c " << c << "\n"; - - auto ljncb = cyl_lbessel_jn(3, cb); - auto logjncb = kyosu::log(cyl_bessel_jn(3, cb)); - TTS_RELATIVE_EQUAL(ljncb, logjncb, tol) << "c " << c << "\n"; - - auto ljncm = cyl_lbessel_jn(3, cm); - auto logjncm = kyosu::log(cyl_bessel_jn(3, cm)); - TTS_RELATIVE_EQUAL(ljncm, logjncm, tol) << "c " << c << "\n"; - TTS_EXPECT(eve::all(kyosu::is_real(kyosu::cyl_lbessel_jn(3, cr)))); -// TTS_EXPECT(eve::all(kyosu::is_pure(ci))); - auto z = kyosu::complex(T(0), T(0)); - TTS_IEEE_EQUAL(cyl_lbessel_jn(3, z), kyosu::complex(eve::minf(eve::as()))); - } -}; From 2a6470ff342fd3425760bc6c5cf33aeddac693e6 Mon Sep 17 00:00:00 2001 From: jtlap Date: Fri, 3 Nov 2023 12:33:04 +0100 Subject: [PATCH 14/17] better doc --- include/kyosu/functions/cyl_bessel_h1n.hpp | 7 ++++--- include/kyosu/functions/cyl_bessel_h2n.hpp | 3 ++- include/kyosu/functions/cyl_bessel_i0.hpp | 9 ++++----- include/kyosu/functions/cyl_bessel_i1.hpp | 9 ++++----- include/kyosu/functions/cyl_bessel_in.hpp | 16 +++++++--------- include/kyosu/functions/cyl_bessel_j0.hpp | 6 +++--- include/kyosu/functions/cyl_bessel_j1.hpp | 5 ++--- include/kyosu/functions/cyl_bessel_jn.hpp | 14 +++++++------- include/kyosu/functions/cyl_bessel_k0.hpp | 7 ++++--- include/kyosu/functions/cyl_bessel_k1.hpp | 4 ++-- include/kyosu/functions/cyl_bessel_kn.hpp | 18 +++++++----------- include/kyosu/functions/cyl_bessel_y0.hpp | 8 +++----- include/kyosu/functions/cyl_bessel_y1.hpp | 8 +++----- include/kyosu/functions/cyl_bessel_yn.hpp | 14 +++++--------- test/doc/cyl_bessel_h1n.cpp | 9 +++------ test/doc/cyl_bessel_h2n.cpp | 11 ++++------- test/doc/cyl_bessel_kn.cpp | 11 ++++------- test/doc/cyl_bessel_yn.cpp | 22 +++++++++++----------- 18 files changed, 79 insertions(+), 102 deletions(-) diff --git a/include/kyosu/functions/cyl_bessel_h1n.hpp b/include/kyosu/functions/cyl_bessel_h1n.hpp index 963f5950..70dea117 100644 --- a/include/kyosu/functions/cyl_bessel_h1n.hpp +++ b/include/kyosu/functions/cyl_bessel_h1n.hpp @@ -43,8 +43,9 @@ namespace kyosu //====================================================================================================================== //! @addtogroup functions //! @{ -//! @brief Computes the Bessel/Hankel functions of the third kind , -//! \f$ H_n^{(1)}\f$. +//! @var cyl_bessel_h1n +//! @brief Computes the Bessel/Hankel functions of the third kind, +//! \f$ H_n^{(1)}(z) = J_n(z)+iY_n(z)\f$. //! //! @code //! #include @@ -66,7 +67,7 @@ namespace kyosu //! //! **Return value** //! -//! * return \f$H_n^{(1)}(z)\f$. +//! * returns \f$H_n^{(1)}(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_h2n.hpp b/include/kyosu/functions/cyl_bessel_h2n.hpp index 64680bcd..d0f04f0d 100644 --- a/include/kyosu/functions/cyl_bessel_h2n.hpp +++ b/include/kyosu/functions/cyl_bessel_h2n.hpp @@ -43,8 +43,9 @@ namespace kyosu //====================================================================================================================== //! @addtogroup functions //! @{ +//! @var cyl_bessel_h2n //! @brief Computes the Bessel/Hankel functions of the third kind , -//! \f$ H_n^{(2)}\f$. +//! \f$ H_n^{(2)} = J_n(z)-iY_n(z)\f$. //! //! @code //! #include diff --git a/include/kyosu/functions/cyl_bessel_i0.hpp b/include/kyosu/functions/cyl_bessel_i0.hpp index cab13917..a5f1a21c 100644 --- a/include/kyosu/functions/cyl_bessel_i0.hpp +++ b/include/kyosu/functions/cyl_bessel_i0.hpp @@ -39,11 +39,10 @@ namespace kyosu //! @addtogroup functions //! @{ //! @var cyl_bessel_i0 -//! @brief Computes the Bessel function of the first kind, -//! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) -//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. +//! @brief Computes the modified Bessel function of the first kind \f$I_{0}(x)=J_{0}(ix)\f$ +//! extended to the complex plane and cayley_dickson algebras. //! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! It is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. //! //! @code //! #include @@ -65,7 +64,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_0(z)\f$. +//! * returns \f$I_0(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_i1.hpp b/include/kyosu/functions/cyl_bessel_i1.hpp index dbc8870e..ecbfe37c 100644 --- a/include/kyosu/functions/cyl_bessel_i1.hpp +++ b/include/kyosu/functions/cyl_bessel_i1.hpp @@ -39,11 +39,10 @@ namespace kyosu //! @addtogroup functions //! @{ //! @var cyl_bessel_i1 -//! @brief Computes the Bessel function of the first kind, -//! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) -//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. +//! @brief Computes the modified Bessel function of the first kind, +//! \f$ I_1(x)= _iJ_1(ix) \f$ extended to the complex plane and cayley_dickson algebras. //! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! It is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 0\f$. //! //! @code //! #include @@ -65,7 +64,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_0(z)\f$. +//! * returns \f$I_1(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_in.hpp b/include/kyosu/functions/cyl_bessel_in.hpp index 0002b675..7bef2fcd 100644 --- a/include/kyosu/functions/cyl_bessel_in.hpp +++ b/include/kyosu/functions/cyl_bessel_in.hpp @@ -39,11 +39,11 @@ namespace kyosu //====================================================================================================================== //! @addtogroup functions //! @{ -//! @brief Computes the Bessel functions of the first kind, -//! \f$ J_{n}(x)=\sum_{p=0}^{\infty}{\frac{(-1)^p}{p!\,\Gamma (p+n +1)}} -//! {\left({x \over 2}\right)}^{2p+n }\f$. +//! @var cyl_bessel_in +//! @brief Computes the modified Bessel functions of the first kind \f$I_{n}(x)=i^{-n}J_{n }(ix)\f$, +//! extended to the complex plane and cayley_dickson algebras. //! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which +//! It is the solution of \f$ x^{2}y''+xy'+(x^2+n^2)y=0\f$ for which //! \f$ y(0) = 0\f$ if \f$n \ne 0\f$ else \f$1\f$. //! //! @code @@ -55,8 +55,8 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto cyl_bessel_in(int n, T z) noexcept; -//! template constexpr T cyl_bessel_in(N n, T z) noexcept; +//! template constexpr auto cyl_bessel_in(int n, T z) noexcept; +//! template constexpr T cyl_bessel_in(int n, T z) noexcept; //! } //! @endcode //! @@ -66,9 +66,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_n(z)\f$. -//! -//! @warning Up to now the cayley_dickson versions have only been imlemented dor scalar int values of n. +//! * returns \f$J_n(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_j0.hpp b/include/kyosu/functions/cyl_bessel_j0.hpp index abb35ae6..c15741cd 100644 --- a/include/kyosu/functions/cyl_bessel_j0.hpp +++ b/include/kyosu/functions/cyl_bessel_j0.hpp @@ -41,9 +41,9 @@ namespace kyosu //! @var cyl_bessel_j0 //! @brief Computes the Bessel function of the first kind, //! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) -//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. +//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson algebras. //! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! It is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. //! //! @code //! #include @@ -65,7 +65,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_0(z)\f$. +//! * returns \f$J_0(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_j1.hpp b/include/kyosu/functions/cyl_bessel_j1.hpp index b6442abb..36613d2f 100644 --- a/include/kyosu/functions/cyl_bessel_j1.hpp +++ b/include/kyosu/functions/cyl_bessel_j1.hpp @@ -42,8 +42,7 @@ namespace kyosu //! @brief Computes the Bessel function of the first kind, //! \f$ J_1(x)=\frac1{\pi }\int _{0}^{\pi}\cos(\tau-x\sin \tau )\,\mathrm {d} \tau \f$ //! extended to the complex plane and cayley_dickson values. -//! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! It is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. //! //! @code //! #include @@ -65,7 +64,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_0(z)\f$. +//! * returns \f$J_0(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_jn.hpp b/include/kyosu/functions/cyl_bessel_jn.hpp index f225676d..09e3ffc0 100644 --- a/include/kyosu/functions/cyl_bessel_jn.hpp +++ b/include/kyosu/functions/cyl_bessel_jn.hpp @@ -39,11 +39,13 @@ namespace kyosu //====================================================================================================================== //! @addtogroup functions //! @{ +//! @var cyl_bessel_jn //! @brief Computes the Bessel functions of the first kind, //! \f$ J_{n}(x)=\sum_{p=0}^{\infty}{\frac{(-1)^p}{p!\,\Gamma (p+n +1)}} -//! {\left({x \over 2}\right)}^{2p+n }\f$. +//! {\left({x \over 2}\right)}^{2p+n }\f$ +//! extended to the complex plane and cayley_dickson values. //! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which +//! It is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which //! \f$ y(0) = 0\f$ if \f$n \ne 0\f$ else \f$1\f$. //! //! @code @@ -55,8 +57,8 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto cyl_bessel_jn(int n, T z) noexcept; -//! template constexpr T cyl_bessel_jn(N n, T z) noexcept; +//! template constexpr auto cyl_bessel_jn(int n, T z) noexcept; +//! template constexpr T cyl_bessel_jn(int n, T z) noexcept; //! } //! @endcode //! @@ -66,9 +68,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_n(z)\f$. -//! -//! @warning Up to now the cayley_dickson versions have only been imlemented dor scalar int values of n. +//! * returns \f$J_n(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_k0.hpp b/include/kyosu/functions/cyl_bessel_k0.hpp index b1bdb2fd..a7f38da3 100644 --- a/include/kyosu/functions/cyl_bessel_k0.hpp +++ b/include/kyosu/functions/cyl_bessel_k0.hpp @@ -39,8 +39,9 @@ namespace kyosu //! @addtogroup functions //! @{ //! @var cyl_bessel_k0 -//! @brief Computes the Bessel function of the second kind, -//! \f$ K_0(x)\f$ extended to the complex plane and cayley_dickson values. +//! @brief Computes the modified Bessel function of the second kind, +//! \f$ K_0(x)=\lim_{\alpha\to 0}{\frac {\pi }{2}}{\frac {I_{-\alpha }(x)-I_{\alpha }(x)}{\sin \alpha \pi }}\f$. +//! extended to the complex plane and cayley_dickson values. //! //! @code //! #include @@ -62,7 +63,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$K_0(z)\f$. +//! * returns \f$K_0(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_k1.hpp b/include/kyosu/functions/cyl_bessel_k1.hpp index e7e3609a..d5c57d8e 100644 --- a/include/kyosu/functions/cyl_bessel_k1.hpp +++ b/include/kyosu/functions/cyl_bessel_k1.hpp @@ -40,7 +40,7 @@ namespace kyosu //! @{ //! @var cyl_bessel_k1 //! @brief Computes the Bessel function of the second kind, -//! \f$ K_1(x)\f$ extended to the complex plane and cayley_dickson values. +//! \f$ K_1(x)\lim_{\alpha\to 1}{\frac {\pi }{2}}{\frac {I_{-\alpha }(x)-I_{\alpha }(x)}{\sin \alpha \pi }}\f$ extended to the complex plane and cayley_dickson values. //! //! @code //! #include @@ -62,7 +62,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$K_1(z)\f$. +//! * returns \f$K_1(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_kn.hpp b/include/kyosu/functions/cyl_bessel_kn.hpp index caa0c924..3283b61c 100644 --- a/include/kyosu/functions/cyl_bessel_kn.hpp +++ b/include/kyosu/functions/cyl_bessel_kn.hpp @@ -39,12 +39,10 @@ namespace kyosu //====================================================================================================================== //! @addtogroup functions //! @{ -//! @brief Computes the Bessel functions of the first kind, -//! \f$ J_{n}(x)=\sum_{p=0}^{\infty}{\frac{(-1)^p}{p!\,\Gamma (p+n +1)}} -//! {\left({x \over 2}\right)}^{2p+n }\f$. -//! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which -//! \f$ y(0) = 0\f$ if \f$n \ne 0\f$ else \f$1\f$. +//! @var cyl_bessel_kn +//! @brief Computes the modified Bessel functions of the second kind, +//! \f$ K_{n}(x)=\lim_{\alpha\to n}{\frac {\pi }{2}}{\frac {I_{-\alpha }(x)-I_{\alpha }(x)}{\sin \alpha \pi }}\f$. +//! extended to the complex plane and cayley_dickson algebras. //! //! @code //! #include @@ -55,8 +53,8 @@ namespace kyosu //! @code //! namespace kyosu //! { -//! template constexpr auto cyl_bessel_kn(int n, T z) noexcept; -//! template constexpr T cyl_bessel_kn(N n, T z) noexcept; +//! template constexpr auto cyl_bessel_kn(int n, T z) noexcept; +//! template constexpr T cyl_bessel_kn(int n, T z) noexcept; //! } //! @endcode //! @@ -66,9 +64,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_n(z)\f$. -//! -//! @warning Up to now the cayley_dickson versions have only been imlemented dor scalar int values of n. +//! * returns \f$K_n(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_y0.hpp b/include/kyosu/functions/cyl_bessel_y0.hpp index d59f5481..1511d819 100644 --- a/include/kyosu/functions/cyl_bessel_y0.hpp +++ b/include/kyosu/functions/cyl_bessel_y0.hpp @@ -43,10 +43,8 @@ namespace kyosu //! @{ //! @var cyl_bessel_y0 //! @brief Computes the Bessel function of the second kind, -//! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) -//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. -//! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! \f$ Y_0(x)=\lim_{\alpha\to 0}{{\frac {J_{\alpha }(x)\cos(\alpha\pi)-J_{-\alpha }(x)}{\sin(\alpha\pi)}}}\f$, +//! extended to the complex plane and cayley_dickson algebras. //! //! @code //! #include @@ -68,7 +66,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_0(z)\f$. +//! * returns \f$Y_0(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_y1.hpp b/include/kyosu/functions/cyl_bessel_y1.hpp index 3af150ab..dc179389 100644 --- a/include/kyosu/functions/cyl_bessel_y1.hpp +++ b/include/kyosu/functions/cyl_bessel_y1.hpp @@ -43,10 +43,8 @@ namespace kyosu //! @{ //! @var cyl_bessel_y1 //! @brief Computes the Bessel function of the second kind, -//! \f$ J_0(x)=\frac1{\pi }\int _{0}^{\pi}\cos(x\sin \tau) -//! \,\mathrm {d} \tau \f$ extended to the complex plane and cayley_dickson values. -//! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+x^2y=0\f$ for which \f$ y(0) = 1\f$. +//! \f$ Y_1(x)=\lim_{\alpha\to 1}{{\frac {J_{\alpha }(x)\cos(\alpha\pi)-J_{-\alpha }(x)}{\sin(\alpha\pi)}}}\f$, +//! extended to the complex plane and cayley_dickson algebras. //! //! @code //! #include @@ -68,7 +66,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_0(z)\f$. +//! * returns \f$Y_1(z)\f$. //! //! @groupheader{Example} //! diff --git a/include/kyosu/functions/cyl_bessel_yn.hpp b/include/kyosu/functions/cyl_bessel_yn.hpp index 703ca4c5..075be242 100644 --- a/include/kyosu/functions/cyl_bessel_yn.hpp +++ b/include/kyosu/functions/cyl_bessel_yn.hpp @@ -39,12 +39,10 @@ namespace kyosu //====================================================================================================================== //! @addtogroup functions //! @{ -//! @brief Computes the Bessel functions of the first kind, -//! \f$ J_{n}(x)=\sum_{p=0}^{\infty}{\frac{(-1)^p}{p!\,\Gamma (p+n +1)}} -//! {\left({x \over 2}\right)}^{2p+n }\f$. -//! -//! In the real field, it is the solution of \f$ x^{2}y''+xy'+(x^2-n^2)y=0\f$ for which -//! \f$ y(0) = 0\f$ if \f$n \ne 0\f$ else \f$1\f$. +//! @var cyl_bessel_yn +//! @brief Computes the modified Bessel functions of the second kind, +//! \f$ Y_n(x)=\lim_{\alpha\to n}{{\frac {J_{\alpha }(x)\cos(\alpha\pi)-J_{-\alpha }(x)}{\sin(\alpha\pi)}}}\f$, +//! extended to the complex plane and cayley_dickson algebras. //! //! @code //! #include @@ -66,9 +64,7 @@ namespace kyosu //! //! **Return value** //! -//! * return the cylindrical \f$J_n(z)\f$. -//! -//! @warning Up to now the cayley_dickson versions have only been imlemented dor scalar int values of n. +//! * returns \f$Y_n(z)\f$. //! //! @groupheader{Example} //! diff --git a/test/doc/cyl_bessel_h1n.cpp b/test/doc/cyl_bessel_h1n.cpp index 1abb3ddd..03a21143 100644 --- a/test/doc/cyl_bessel_h1n.cpp +++ b/test/doc/cyl_bessel_h1n.cpp @@ -5,16 +5,13 @@ int main() { std::cout.precision(16); - // using w_t = eve::wide>; - //auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); - auto z = kyosu::complex(30.0, 0.0); -// auto z = kyosu::complex(1.0); + using w_t = eve::wide>; + auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); for(int n=2; n <= 2; ++n) { auto h1n = kyosu::cyl_bessel_h1n(n, z); std::cout << "n " << n << " z " << z << std::endl; - std::cout << "1 " << h1n << std::endl; - std::cout << "2 " << kyosu::cyl_bessel_h1n(n, 30.0) << std::endl; + std::cout << "yl_bessel_h1n(n, z) " << h1n << std::endl; } return 0; } diff --git a/test/doc/cyl_bessel_h2n.cpp b/test/doc/cyl_bessel_h2n.cpp index 62c4f254..b6bf7701 100644 --- a/test/doc/cyl_bessel_h2n.cpp +++ b/test/doc/cyl_bessel_h2n.cpp @@ -5,16 +5,13 @@ int main() { std::cout.precision(16); - // using w_t = eve::wide>; - //auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); - auto z = kyosu::complex(30.0, 0.0); -// auto z = kyosu::complex(1.0); + using w_t = eve::wide>; + auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); for(int n=2; n <= 2; ++n) { - auto h2n = kyosu::cyl_bessel_h2n(n, z); + auto h1n = kyosu::cyl_bessel_h2n(n, z); std::cout << "n " << n << " z " << z << std::endl; - std::cout << "0 " << h2n << std::endl; - std::cout << "2 " << kyosu::cyl_bessel_h2n(n, 30.0) << std::endl; + std::cout << "yl_bessel_h2n(n, z) " << h1n << std::endl; } return 0; } diff --git a/test/doc/cyl_bessel_kn.cpp b/test/doc/cyl_bessel_kn.cpp index 2a6a5193..f255f069 100644 --- a/test/doc/cyl_bessel_kn.cpp +++ b/test/doc/cyl_bessel_kn.cpp @@ -5,16 +5,13 @@ int main() { std::cout.precision(16); - // using w_t = eve::wide>; - //auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); - auto z = kyosu::complex(0.0, 3.0); -// auto z = kyosu::complex(1.0); + using w_t = eve::wide>; + auto z = kyosu::complex(w_t(20.0, 1.5), w_t(0.0, 1.5)); for(int n=2; n <= 10; ++n) { auto kn = kyosu::cyl_bessel_kn(n, z); - std::cout << "n " << n << " z " << z << std::endl; - std::cout << "1 " << kn << std::endl; - std::cout << "2 " << kyosu::cyl_bessel_kn(n, 4.0) << std::endl; + std::cout << "n " << n << ", z " << z << std::endl; + std::cout << "cyl_bessel_kn(n, z) " << kn << std::endl; } return 0; } diff --git a/test/doc/cyl_bessel_yn.cpp b/test/doc/cyl_bessel_yn.cpp index 13c147a4..bc8bf246 100644 --- a/test/doc/cyl_bessel_yn.cpp +++ b/test/doc/cyl_bessel_yn.cpp @@ -17,18 +17,18 @@ int main() 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"; + 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"; + std::cout << "1.0 " << " -> " << cyl_bessel_yn(n,1.0) << "\n"; + std::cout << "1.0+36.0i " << " -> " << cyl_bessel_yn(n,kyosu::complex(1.0, 36.0)) << "\n"; + std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_bessel_yn(n,kyosu::quaternion(1.0, 36.0, 2.0, 1.5)) << "\n"; -// eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); -// auto z0 = kyosu::complex(z1, z2); -// std::cout << z0 << " \n-> " << cyl_bessel_yn(n,z0) << "\n"; -// return 0; + eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); + auto z0 = kyosu::complex(z1, z2); + std::cout << z0 << " \n-> " << cyl_bessel_yn(n,z0) << "\n"; + return 0; } From 5975429faeb3c9d63c5340688570bb16f299411c Mon Sep 17 00:00:00 2001 From: jtlap Date: Fri, 3 Nov 2023 12:38:08 +0100 Subject: [PATCH 15/17] suppressing cout.precision --- test/unit/function/cyl_bessel_i0.cpp | 2 -- test/unit/function/cyl_bessel_i1.cpp | 2 -- test/unit/function/cyl_bessel_in.cpp | 1 - test/unit/function/cyl_bessel_j0.cpp | 2 -- test/unit/function/cyl_bessel_j1.cpp | 2 -- test/unit/function/cyl_bessel_jn.cpp | 1 - test/unit/function/cyl_bessel_k0.cpp | 1 - test/unit/function/cyl_bessel_k1.cpp | 1 - test/unit/function/cyl_bessel_kn.cpp | 1 - test/unit/function/cyl_bessel_y0.cpp | 1 - test/unit/function/cyl_bessel_y1.cpp | 1 - test/unit/function/cyl_bessel_yn.cpp | 1 - 12 files changed, 16 deletions(-) diff --git a/test/unit/function/cyl_bessel_i0.cpp b/test/unit/function/cyl_bessel_i0.cpp index 0d8778ea..e9814629 100644 --- a/test/unit/function/cyl_bessel_i0.cpp +++ b/test/unit/function/cyl_bessel_i0.cpp @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i0 over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{1.3556565953650257e+00, 3.3014253671480853e+00, -2.5475601159592665e+00, 1.1681025625526764e+00, 2.0457564384951734e+00, -1.2532453407268573e+00, -4.7295944112409973e+00, -2.2840688135673504e+00}; std::array im{4.0914290638525896e+00, -5.6356545999116547e-01, -3.6064274257830387e+00, 1.6951171853616265e+00, -1.3703722677989172e+00, -1.2563658033049907e+00, -6.3364633643107848e-01, 3.7665729797577785e+00}; @@ -34,7 +33,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i0 over real" TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_i0(c), res, 1.0e-7) << i << " <- " << c << '\n'; } } - std::cout.precision(17); std::array re{1.1554018562666468e-01, 5.3044189227319072e-01, 3.5245310048347744e-01, 3.2527350287655354e-01, 2.4688936247411308e+00, -5.1146428618899231e-01, 1.6635701653597774e+00, -1.7452778188732947e+00}; diff --git a/test/unit/function/cyl_bessel_i1.cpp b/test/unit/function/cyl_bessel_i1.cpp index f0a20c3c..84f2bf87 100644 --- a/test/unit/function/cyl_bessel_i1.cpp +++ b/test/unit/function/cyl_bessel_i1.cpp @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i1 over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{-4.1836770228521267e+01, -4.1682207527601065e+01, -4.0619579092445747e+01, -4.4234864207453732e+01, -4.3110416350751848e+01, -4.3626638817587562e+01, -4.4580540281737981e+01, -4.1708122005195847e+01}; std::array im{-4.4724569983549614e+01, -4.6490400566267638e+01, -4.7317594550520738e+01, -4.1337120165516517e+01, -4.1276020002765804e+01, -4.6670524375443286e+01, -4.9467141247871425e+01, -4.1317766188944873e+01}; @@ -33,7 +32,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_i1 over real" TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_i1(c), res, 1.0e-7) << i << " <- " << c << '\n'; } } - std::cout.precision(17); std::array re{-2.0045932610658785e+00, 1.5610767253921298e-02, 1.2016964006673581e+00, -1.4746256231752448e+00, 3.8509464414269778e-01, 4.3119169358173060e-01, 1.8117761375602393e-01, 2.2535543048838065e+00}; std::array im{-1.3262936608790188e+00, 7.3358653580584687e-01, -7.7829729505626277e-01, 1.8039445371160530e+00, -3.9048140529175879e-01, 2.3733805921203750e+00, 4.1558677759532259e-01, -2.2418145868867645e-01}; diff --git a/test/unit/function/cyl_bessel_in.cpp b/test/unit/function/cyl_bessel_in.cpp index 42aa423c..96aea3b3 100644 --- a/test/unit/function/cyl_bessel_in.cpp +++ b/test/unit/function/cyl_bessel_in.cpp @@ -16,7 +16,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_in over real" { if constexpr(sizeof(T) == 8) { - std::cout.precision(17); std::array re{-4.2214088790474964e+00, -7.8055061036477191e-01, 3.4677730339081023e+00, -3.5357366395596843e+00, 2.6849636588863888e+00, 4.1415226198517496e+00, -4.3623200824031505e+00, 1.0341209529062589e-01}; std::array im{1.0956977560286330e+00, -7.5270747683236494e-01, -2.3663854475056922e+00, 2.1367696130628167e-01, 3.8676287626125005e+00, 1.2360189302485036e+00, -1.6272444761141791e+00, 2.3032771326096082e+00}; diff --git a/test/unit/function/cyl_bessel_j0.cpp b/test/unit/function/cyl_bessel_j0.cpp index 0302271b..581e33a0 100644 --- a/test/unit/function/cyl_bessel_j0.cpp +++ b/test/unit/function/cyl_bessel_j0.cpp @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j0 over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{ e_t(-9.8403264736233908e+01), e_t(8.6287324687322851e+01), e_t(2.8677183208357349e+01), e_t(7.5642931882481321e+01), e_t(9.6351793077594323e+01) , e_t(-7.8716695947540288e+01), e_t(-5.0019457963533974e+01), e_t(5.9030075032843811e+01), e_t(-7.5596124248257235e+00), e_t(8.5810875418338028e+01)}; @@ -37,7 +36,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j0 over real" TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_j0(im[i]), kyosu::real(kyosu::cyl_bessel_j0(kyosu::complex(im[i], e_t(0.0)))), 1.0e-4)<< re[i] << '\n'; } } - std::cout.precision(17); std::array re{2.175102683907268e+00, 3.007738475083643e-01, 3.571791197760752e+00, 1.916834892035037e+00, 2.309329011574615e+00, 1.116302757926106e+00, 3.454188251414261e+00, 2.504516854443353e+00}; std::array im{ 3.008362577802520e+00, 3.733804275914350e+00, 2.988312632706299e+00, 4.411110052057489e+00, 3.055249958745497e+00, 2.508994160634028e+00, 3.100976745146938e-01, 3.854880352108083e+00}; diff --git a/test/unit/function/cyl_bessel_j1.cpp b/test/unit/function/cyl_bessel_j1.cpp index 5e3ce470..2de8f892 100644 --- a/test/unit/function/cyl_bessel_j1.cpp +++ b/test/unit/function/cyl_bessel_j1.cpp @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j1 over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{ -9.8403264736233908e+01, 8.6287324687322851e+01, 2.8677183208357349e+01, 7.5642931882481321e+01, 9.6351793077594323e+01 , -7.8716695947540288e+01, -5.0019457963533974e+01, 5.9030075032843811e+01, -7.5596124248257235e+00, 8.5810875418338028e+01}; std::array im{-1.9271037994246587e+01, 7.9043790847200569e+01, -3.0599555605285602e+01, -7.1720588923445590e+01, -6.1110206702938363e+01 @@ -34,7 +33,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_j1 over real" TTS_RELATIVE_EQUAL(kyosu::cyl_bessel_j1(c), res, 1.0e-7) << i << " <- " << c << '\n'; } } - std::cout.precision(17); std::array re{2.175102683907268e+00, 3.007738475083643e-01, 3.571791197760752e+00, 1.916834892035037e+00, 2.309329011574615e+00, 1.116302757926106e+00, 3.454188251414261e+00, 2.504516854443353e+00}; std::array im{ 3.008362577802520e+00, 3.733804275914350e+00, 2.988312632706299e+00, 4.411110052057489e+00, 3.055249958745497e+00, 2.508994160634028e+00, 3.100976745146938e-01, 3.854880352108083e+00}; diff --git a/test/unit/function/cyl_bessel_jn.cpp b/test/unit/function/cyl_bessel_jn.cpp index a0c02310..a8c7c70c 100644 --- a/test/unit/function/cyl_bessel_jn.cpp +++ b/test/unit/function/cyl_bessel_jn.cpp @@ -16,7 +16,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_jn over real" { if constexpr(sizeof(T) == 8) { - std::cout.precision(17); std::array re{ -9.8403264736233908e+01, 8.6287324687322851e+01, 2.8677183208357349e+01, 7.5642931882481321e+01, 9.6351793077594323e+01 , -7.8716695947540288e+01, -5.0019457963533974e+01, 5.9030075032843811e+01, -7.5596124248257235e+00, 8.5810875418338028e+01}; std::array im{-1.9271037994246587e+01, 7.9043790847200569e+01, -3.0599555605285602e+01, -7.1720588923445590e+01, -6.1110206702938363e+01 diff --git a/test/unit/function/cyl_bessel_k0.cpp b/test/unit/function/cyl_bessel_k0.cpp index acbd42a9..b4ac046b 100644 --- a/test/unit/function/cyl_bessel_k0.cpp +++ b/test/unit/function/cyl_bessel_k0.cpp @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k0 over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{2.3584877349134130e+00, 9.0164385767932564e-01, -1.5393599740361523e+00, 1.5668095462351217e+00, 1.6179702389500128e+00, -5.5681137226619104e-01, -2.2990386896893691e+00, 1.9300383451287466e+00}; std::array im{-8.0071372624748560e-01, 1.9823595245805736e+00, -1.3978761922193099e+00, -1.4654991890321578e+00, -5.3416883666219606e-01, -1.9265301646062505e+00, -1.6132600067568537e+00, -2.3575615216282779e+00}; diff --git a/test/unit/function/cyl_bessel_k1.cpp b/test/unit/function/cyl_bessel_k1.cpp index fa2e5574..64afaefe 100644 --- a/test/unit/function/cyl_bessel_k1.cpp +++ b/test/unit/function/cyl_bessel_k1.cpp @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_k1 over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{-7.3954299670486945e-01, -3.7928780215327196e-02, -1.4509312401953633e+00, 2.2964197111982774e+00, -7.9699549838500194e-01, 5.4099472045802843e-01, -5.8957905551605738e-01, -5.1146527014834298e-01}; std::array im{-2.1345822747239773e+00, -3.3484907524396035e-01, 2.2665940322875286e+00, -1.0045032564810230e+00, -2.4392625060316648e+00, 2.6547418070570761e-01, -2.3352892270792931e+00, 1.3879256025736524e+00}; diff --git a/test/unit/function/cyl_bessel_kn.cpp b/test/unit/function/cyl_bessel_kn.cpp index 21d7c910..401e0417 100644 --- a/test/unit/function/cyl_bessel_kn.cpp +++ b/test/unit/function/cyl_bessel_kn.cpp @@ -16,7 +16,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_in over real" { if constexpr(sizeof(T) == 8) { - std::cout.precision(17); std::array re{1.4376454063874640e+00, -1.0673204704071071e-01, -1.0828666774506268e+00, -7.9765980000519210e-01, 4.2409189866361996e-01,7.6608496996954323e-01, -2.2849830522845975e+00, 4.3813679108628711e-01}; std::array im{-2.5175040147125849e-01, 6.2364116715495577e-01, 1.1653810042557073e+00, 1.3466103438716175e+00, -1.8438471344336298e+00, -2.0952595273684760e+00, -7.7012218588605907e-01, 1.3705948670934771e+00}; diff --git a/test/unit/function/cyl_bessel_y0.cpp b/test/unit/function/cyl_bessel_y0.cpp index 85085b41..7b8f373f 100644 --- a/test/unit/function/cyl_bessel_y0.cpp +++ b/test/unit/function/cyl_bessel_y0.cpp @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_y0 over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{2.4955917718518910e+00, -1.6982259765731111e+00, 1.2315182267355251e+00, 2.1258824905764118e+00, -2.1711533910447676e+00, 9.1742698952666013e-01, -2.5479731910545422e-02, 1.4924891723416751e+00}; std::array im{-5.0148281002365624e-01, 2.0253551112861641e+00, -5.0175898965530097e-01, -4.0455814438323523e-01, -2.1923016231568171e-01, -1.3185119199506317e+00, 1.0779142409536568e+00, 1.2528910061509635e-01}; diff --git a/test/unit/function/cyl_bessel_y1.cpp b/test/unit/function/cyl_bessel_y1.cpp index 11a37499..04593e67 100644 --- a/test/unit/function/cyl_bessel_y1.cpp +++ b/test/unit/function/cyl_bessel_y1.cpp @@ -17,7 +17,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_y1 over real" using e_t = eve::element_type_t; if constexpr(sizeof(e_t) == 8) { - std::cout.precision(17); std::array re{6.5758573026995004e-01, -2.4649387808542378e+00, 1.9205766369481703e+00, -6.9739731404600436e-01, 4.4532863650021348e-01, 2.4757649834958393e+00, -5.3435309039155832e-01, -2.4349028761454212e+00}; std::array im{1.2035170114779326e+00, -3.5932055457567835e-01, 8.3662582734559288e-01, 2.2469103234867700e+00, -9.1023149008937376e-01, 2.4818838321237546e-01, 1.5847819595769823e+00, -1.3967310652706204e+00}; diff --git a/test/unit/function/cyl_bessel_yn.cpp b/test/unit/function/cyl_bessel_yn.cpp index c9ad9e41..21e43ba2 100644 --- a/test/unit/function/cyl_bessel_yn.cpp +++ b/test/unit/function/cyl_bessel_yn.cpp @@ -16,7 +16,6 @@ TTS_CASE_WITH ( "Check kyosu::cyl_bessel_yn over real" { if constexpr(sizeof(T) == 8) { - std::cout.precision(17); std::array re{2.0309144503737393e+00, -1.5992515729134549e+00, 1.6982551054108397e+00, -1.3327549270623136e+00, 2.6585368919111341e-01, 7.3163244235853131e-01, -1.5491151507087859e+00, 2.4846599242220186e+00}; std::array im{-1.4137922349350807e+00, -1.7357125614937789e+00, 4.0762417134967421e-02, -1.4645033242641059e+00, 1.5914290999536518e+00, 2.3057243858262124e+00, -2.2075167685839387e+00, -5.7866716195750623e-01}; From 3459dce7a6539572dc4d29f95dc54e25f645a415 Mon Sep 17 00:00:00 2001 From: jtlap Date: Fri, 3 Nov 2023 12:42:32 +0100 Subject: [PATCH 16/17] spurious comments --- include/kyosu/types/impl/detail/bessely1.hpp | 3 -- test/doc/besselin.hpp | 33 -------------------- 2 files changed, 36 deletions(-) delete mode 100644 test/doc/besselin.hpp diff --git a/include/kyosu/types/impl/detail/bessely1.hpp b/include/kyosu/types/impl/detail/bessely1.hpp index 86d49380..d2705a28 100644 --- a/include/kyosu/types/impl/detail/bessely1.hpp +++ b/include/kyosu/types/impl/detail/bessely1.hpp @@ -21,10 +21,7 @@ namespace kyosu::_ using u_t = eve::underlying_type_t; auto twoopi = eve::two_o_pi(eve::as()); auto r1 = R(1, z); -// std::cout << "r1 " << r1 << std::endl; auto y0 = cyl_bessel_y0(z); -// std::cout << "y0 " << y0 << std::endl; -// std::cout << "j0 " << cyl_bessel_j0(z) << std::endl; 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/test/doc/besselin.hpp b/test/doc/besselin.hpp deleted file mode 100644 index 3e05cf88..00000000 --- a/test/doc/besselin.hpp +++ /dev/null @@ -1,33 +0,0 @@ -#include -#include -#include - -int main() -{ - using kyosu::cyl_bessel_j0; - std::cout.precision(16); - std::cout << std::scientific << std::endl; - std::cout << "1.0 " << " -> " << cyl_bessel_j0(1.0) << "\n"; - std::cout << "1+0i " << " -> " << cyl_bessel_j0(kyosu::complex(1.0, 0.0)) << "\n"; - std::cout << "15.0 " << " -> " << cyl_bessel_j0(15.0) << "\n"; - std::cout << "15+0i " << " -> " << cyl_bessel_j0(kyosu::complex(15.0, 0.0)) << "\n"; - std::cout << "40.0 " << " -> " << cyl_bessel_j0(40.0) << "\n"; - std::cout << "40+0i " << " -> " << cyl_bessel_j0(kyosu::complex(40.0, 0.0)) << "\n"; - std::cout << "60.0 " << " -> " << cyl_bessel_j0(60.0) << "\n"; - std::cout << "60+0i " << " -> " << cyl_bessel_j0(kyosu::complex(60.0, 0.0)) << "\n"; - - eve::wide> z(1.0, 15.0, 40.0, 60.0); - auto zz = kyosu::complex(z); - std::cout << z << " -> \n" << cyl_bessel_j0(z) << "\n"; - std::cout << zz << " -> \n" << cyl_bessel_j0(zz) << "\n"; - - std::cout << "1.0 " << " -> " << cyl_bessel_j0(1.0) << "\n"; - std::cout << "1.0+36.0i " << " -> " << cyl_bessel_j0(kyosu::complex(1.0, 36.0)) << "\n"; - std::cout << "1.0+36.0i+2.0j+1.5k " << " -> " << cyl_bessel_j0(kyosu::quaternion(1.0, 36.0, 2.0, 1.5)) << "\n"; - - - eve::wide> z1(1.0, 2.0, 40.0, 0.0), z2(36.0, 0.5, 0.0, 40.0); - auto z0 = kyosu::complex(z1, z2); - std::cout << z0 << "\n -> " << cyl_bessel_j0(z0) << "\n"; - return 0; -} From 0baf431b710ef5f1f5a2d74d86a22b7e445941df Mon Sep 17 00:00:00 2001 From: jtlap Date: Fri, 3 Nov 2023 14:32:17 +0100 Subject: [PATCH 17/17] rstest suppressed --- test/doc/rstest.cpp | 40 ---------------------------------------- 1 file changed, 40 deletions(-) delete mode 100644 test/doc/rstest.cpp diff --git a/test/doc/rstest.cpp b/test/doc/rstest.cpp deleted file mode 100644 index 0e5ebb73..00000000 --- a/test/doc/rstest.cpp +++ /dev/null @@ -1,40 +0,0 @@ -#include -#include -#include -//#include - -int main() -{ - std::cout.precision(16); - size_t n = 10; - auto z = kyosu::complex(1.0, 2.0); - using c_t = decltype(z); - auto rs = kyosu::_::Rs(n, z); - std::vector < c_t> rrs(n); - for(size_t i=0; i < n; ++i) - { - rrs[i] = kyosu::_::R(i+1, z); - std::cout << "rs " << rs[i] << std::endl; - std::cout << "rrs " << rrs[i] << std::endl; - std::cout << "j/jm1 " << kyosu::cyl_bessel_jn(i, z)/kyosu::cyl_bessel_jn(i+1, z)<< std::endl; - } - auto js = kyosu::_::Js(n, z); - for(size_t i=1; i < n; ++i) - { - std::cout << i << " -> js " << js[i] << std::endl; - std::cout << i << " -> jn " << kyosu::cyl_bessel_jn(i, z) << std::endl; - } -// std::vector js(n+1); -// for(size_t i=0; i < n+1; ++i) -// { -// js[i] = kyosu::cyl_bessel_jn(i, z); -// } -// for(size_t i=1; i < n; ++i) -// { -// std::cout << i << " -> rs " << rs[i] << std::endl; -// std::cout << i << " -> jsp1/rs " << js[i]/rs[i-1] << std::endl; -// std::cout << i << " -> js " << js[i]/js[i-1] << std::endl; -// } - - return 0; -}