diff --git a/cmake/dependencies.cmake b/cmake/dependencies.cmake index 11096f42..919175b9 100644 --- a/cmake/dependencies.cmake +++ b/cmake/dependencies.cmake @@ -38,7 +38,7 @@ CPMAddPackage ( NAME TTS GITHUB_REPOSITORY jfalcou/tts "TTS_QUIET ON" ) CPMAddPackage ( NAME EVE GITHUB_REPOSITORY jfalcou/eve - GIT_TAG 955556c0e679c6861e16a01f0cf00b963ece00f0 + GIT_TAG main OPTIONS "EVE_BUILD_TEST OFF" "EVE_BUILD_BENCHMARKS OFF" "EVE_BUILD_DOCUMENTATION OFF" diff --git a/include/kyosu/details/bessel/bessel_utils.hpp b/include/kyosu/details/bessel/bessel_utils.hpp index 440560a4..19f4f522 100644 --- a/include/kyosu/details/bessel/bessel_utils.hpp +++ b/include/kyosu/details/bessel/bessel_utils.hpp @@ -95,7 +95,7 @@ namespace kyosu::_ D = kyosu::rec(D); delta = C*D; f = f * delta; - } while(eve::any((kyosu::abs(kyosu::dec(delta)) > terminator) && --counter)); + } while(eve::any((kyosu::abs(kyosu::dec(delta)) > terminator)) && --counter); g.reset(); return f; } diff --git a/include/kyosu/details/bessel/besseln/cb_jyn.hpp b/include/kyosu/details/bessel/besseln/cb_jyn.hpp index e6af48e2..904fce0c 100644 --- a/include/kyosu/details/bessel/besseln/cb_jyn.hpp +++ b/include/kyosu/details/bessel/besseln/cb_jyn.hpp @@ -304,7 +304,7 @@ namespace kyosu::_ if (eve::any(rzlt0)) { auto sgn1 = eve::if_else(izgt0, u_t(1), u_t(-1)); - r = if_else(rzlt0, (r+2*muli(sgn1*cb_j0(z))), r); + r = if_else(rzlt0, (r+2*muli(sgn1*_::cb_j0(z))), r); } return r; } @@ -318,8 +318,8 @@ namespace kyosu::_ using u_t = eve::underlying_type_t; auto twoopi = eve::two_o_pi(eve::as()); auto r1 = R(1, z); - auto y0 = cb_y0(z); - auto recs1 = rec(r1)-twoopi/(z*cb_j0(z)*y0); + auto y0 = kyosu::_::cb_y0(z); + auto recs1 = rec(r1)-twoopi/(z* kyosu::_::cb_j0(z)*y0); return if_else(is_eqz(z), complex(eve::minf(eve::as())), y0*recs1); } @@ -331,15 +331,15 @@ namespace kyosu::_ { if ( is_eqz(nn) ) { - return cb_j0(z); + return kyosu::_::cb_j0(z); } else if (nn == 1) { - return cb_j1(z); + return kyosu::_::cb_j1(z); } else if ( nn == -1 ) { - return -cb_j1(z); + return - kyosu::_::cb_j1(z); } else { @@ -351,8 +351,8 @@ namespace kyosu::_ auto srz = eve::signnz(real(z)); z *= srz; - auto j0 = cb_j0(z); - auto j1 = cb_j1(z); + auto j0 = kyosu::_::cb_j0(z); + auto j1 = kyosu::_::cb_j1(z); auto forward = [n, j0, j1](auto z){ auto b0 = j0; @@ -419,14 +419,14 @@ namespace kyosu::_ using u_t = eve::underlying_type_t; if (n <= 1) { - cjv[0] = cb_j0(z); - cyv[0] = cb_y0(z); + cjv[0] = kyosu::_::cb_j0(z); + cyv[0] = kyosu::_::cb_y0(z); if ( is_eqz(nn) ) { return kumi::tuple{cjv[0], cyv[0]}; } - cjv[1] = cb_j1(z); - cyv[1] = cb_y1(z); + cjv[1] = kyosu::_::cb_j1(z); + cyv[1] = kyosu::_::cb_y1(z); if (nn == 1) { return kumi::tuple{cjv[1], cyv[1]}; @@ -445,10 +445,10 @@ namespace kyosu::_ auto izgt0 = eve::is_gtz(imag(z)); z = if_else(rzle0, -z, z);//real(z) is now positive auto rz = rec(z); - cjv[0] = cb_j0(z); - cjv[1] = cb_j1(z); - cyv[0] = cb_y0(z); - cyv[1] = cb_y1(z); + cjv[0] = kyosu::_::cb_j0(z); + cjv[1] = kyosu::_::cb_j1(z); + cyv[0] = kyosu::_::cb_y0(z); + cyv[1] = kyosu::_::cb_y1(z); auto forwardj = [n, rz, &cjv](auto z){ auto bkm2 = cjv[0]; @@ -607,7 +607,7 @@ namespace kyosu::_ auto cb_yn(N n, Z z) noexcept { auto dummy = std::span(); - return cb_yn(n, z, dummy); + return _::cb_yn(n, z, dummy); } //===------------------------------------------------------------------------------------------- diff --git a/include/kyosu/details/bessel/besselr/cb_ikr.hpp b/include/kyosu/details/bessel/besselr/cb_ikr.hpp index 42071d9f..3040ce16 100644 --- a/include/kyosu/details/bessel/besselr/cb_ikr.hpp +++ b/include/kyosu/details/bessel/besselr/cb_ikr.hpp @@ -115,7 +115,7 @@ namespace kyosu::_ auto argzpos = eve::is_gtz(kyosu::arg(z)); z = if_else(argzpos, conj(z), z); auto fac = exp_ipi(-v/2); - auto r = fac*cb_jr(v,kyosu::muli(z)); + auto r = fac*_::cb_jr(v,kyosu::muli(z)); return if_else(argzpos, conj(r), r); } else diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index a6e83cf1..bc44840a 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -38,6 +38,7 @@ #include #include #include +#include #include #include #include @@ -72,6 +73,12 @@ #include #include #include +#include +#include +#include +#include +#include +#include #include #include #include @@ -115,6 +122,7 @@ #include #include #include +#include #include #include #include diff --git a/include/kyosu/functions/am.hpp b/include/kyosu/functions/am.hpp new file mode 100644 index 00000000..ad04a111 --- /dev/null +++ b/include/kyosu/functions/am.hpp @@ -0,0 +1,128 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include +#include +#include + +namespace kyosu +{ + + template + struct am_t : eve::strict_elementwise_callable + { + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b) const noexcept + { return (*this)(kyosu::complex(a), b); } + + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b) const noexcept + { return KYOSU_CALL(a, b); } + + + KYOSU_CALLABLE_OBJECT(am_t, am_); + }; + +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var am +//! @brief Computes Jacobi's Amplitude function. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! // Regular overload +//! constexpr Z am( real z, real m) noexcept; //1 +//! constexpr auto am(complex z, real m) noexcept; //1 +//! +//! //Semantic modifiers +//! constexpr Z am[modular]( real z, real alpha) noexcept; //1 +//! constexpr auto am[modular](complex z, real alpha) noexcept; //1 +//! constexpr Z am[eccentric]( real z, real m) noexcept; //1 +//! constexpr auto am[eccentric](complex z, real m) noexcept; //1 +//! constexpr auto am[threshold = tol]( real z, real k) noexcept; //1 +//! constexpr auto am[threshold = tol](complex z, real k) noexcept; //1 +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `u`: argument. +//! * `m`: amplitude parameter (\f$0\le m\le 1). +//! * `alpha `: modular angle in radian. +//! * `k`: elliptic modulus (eccentricity) . +//! * `tol': accuracy tolerance (by defaut [epsilon](@ref eve::epsilon). +//! * `c`: [Conditional expression](@ref eve::conditional_expr) masking the operation. +//! * `l`: [Logical value](@ref eve::logical_value) masking the operation. +//! +//! +//! **Return value** +//! +//! 1. return the jacobian amplitude function. Take care that the meaning of the second parameters +//! depends on the option used (see note below). +//! +//! @note +//! * \f$\alpha\f$ is named the modular angle given in radian (modular option). +//! * \f$ k = \sin\alpha \f$ is named the elliptic modulus or eccentricity (eccentric option). +//! * \f$ m = k^2 = \sin^2\alpha\f$ is named the parameter (no option). +//! Each of the above three quantities is completely determined by any of the others (given that they are non-negative). +//! Thus, they can be used interchangeably up to roundings errors by giving the right option. +//! +//! @groupheader{External references} +//! * [C++ standard reference: am](https://en.cppreference.com/w/cpp/numeric/special_functions/am) +//! * [DLMF: Jacobi Amplitude](https://dlmf.nist.gov/22.16) +//! * [Wolfram MathWorld: Jacobi Amplitude](https://mathworld.wolfram.com/JacobiAmplitude.html) +//! * [Wikipedia: Jacobi elliptic functions](https://en.wikipedia.org/wiki/Jacobi_elliptic_functions) +//! +//! @groupheader{Example} +//! @godbolt{doc/am.cpp} +//====================================================================================================================== + inline constexpr auto am = eve::functor; +//====================================================================================================================== +//! @} +//====================================================================================================================== +} + +namespace kyosu::_ +{ + template + KYOSU_FORCEINLINE constexpr auto am_(KYOSU_DELAY(), O const& o, Z u, M m) noexcept + { + auto tol = [&](){ + if constexpr (O::contains(eve::threshold)) return o[eve::threshold].value(m); + else return eve::epsilon(eve::maximum(eve::abs(m))); + }(); + m = eve::abs(m); + if (O::contains(eve::modular)) m = eve::sin(m); + else if (O::contains(eve::eccentric)) m = eve::sqrt(m); + if (eve::all(is_real(u))) return kyosu::complex(eve::am(kyosu::real(u), m)); + + auto [phi, psi] = u; + auto [s,c,d] = eve::jacobi_elliptic[eve::threshold = tol](phi, m); + auto mc = eve::sqrt(eve::oneminus(eve::sqr(m))); + auto [s1,c1,d1] = eve::jacobi_elliptic[eve::threshold = tol](psi,mc); + auto x0 = s*d1; + auto x1 = c*c1; + auto y = s1*d; + auto km = eve::ellint_1(m); + auto nx = eve::floor ((phi+2*km )/(4*km)); + return kyosu::complex(eve::atan2[eve::pedantic](x0,x1), eve::atanh(y)) + eve::two_pi(as(phi))*nx; + } +} diff --git a/include/kyosu/functions/average.hpp b/include/kyosu/functions/average.hpp index 233eaf8a..59729b2d 100644 --- a/include/kyosu/functions/average.hpp +++ b/include/kyosu/functions/average.hpp @@ -49,7 +49,7 @@ namespace kyosu //! //! **Parameters** //! -//! * `z0`, `z1...`: Values to process. Can be a mix of complex and real floating values and complex values. +//! * `z0`, `z1...`: Values to process. Can be a mix of complex and real floating values. //! //! **Return value** //! diff --git a/include/kyosu/functions/ellint_fe.hpp b/include/kyosu/functions/ellint_fe.hpp new file mode 100644 index 00000000..8b4ecdb0 --- /dev/null +++ b/include/kyosu/functions/ellint_fe.hpp @@ -0,0 +1,163 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include +#include + +namespace kyosu +{ + + template + struct ellint_fe_t : eve::strict_elementwise_callable + { + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b) const noexcept + { return eve::zip(eve::ellint_1(a, b), eve::ellint_2(a, b)); } + + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b) const noexcept -> decltype(eve::zip(T0(), T0())) + { return KYOSU_CALL(a, b); } + + + KYOSU_CALLABLE_OBJECT(ellint_fe_t, ellint_fe_); + }; + +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var ellint_fe +//! @brief Computes Jacobi's Amplitude function. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! // Regular overload +//! constexpr Z ellint_fe( real z, real m) noexcept; //1 +//! constexpr auto ellint_fe(complex z, real m) noexcept; //1 +//! +//! //Semantic modifiers +//! constexpr Z ellint_fe[modular]( real z, real alpha) noexcept; //1 +//! constexpr auto ellint_fe[modular](complex z, real alpha) noexcept; //1 +//! constexpr Z ellint_fe[eccentric]( real z, real m) noexcept; //1 +//! constexpr auto ellint_fe[eccentric](complex z, real m) noexcept; //1 +//! constexpr auto ellint_fe[threshold = tol]( real z, real k) noexcept; //1 +//! constexpr auto ellint_fe[threshold = tol](complex z, real k) noexcept; //1 +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `u`: argument. +//! * `m`: amplitude parameter (\f$0\le m\le 1). +//! * `alpha `: modular angle in radian. +//! * `k`: elliptic modulus (eccentricity) . +//! * `tol': accuracy tolerance (by defaut [epsilon](@ref eve::epsilon). +//! * `c`: [Conditional expression](@ref eve::conditional_expr) masking the operation. +//! * `l`: [Logical value](@ref eve::logical_value) masking the operation. +//! +//! +//! **Return value** +//! +//! 1. return the elliptic incomplete functions or first and second kind. Tafe care that the meaning of the second parameters +//! depends on the option used (see note below). +//! +//! @note +//! * \f$\alpha\f$ is named the modular angle given in radian (modular option). +//! * \f$ k = \sin\alpha \f$ is named the elliptic modulus or eccentricity (eccentric option). +//! * \f$ m = k^2 = \sin^2\alpha\f$ is named the parameter (no option). +//! Each of the above three quantities is completely determined by any of the others (given that they are non-negative). +//! Thus, they can be used interchangeably up to roundings errors by giving the right option. +//! +//! @groupheader{External references} +//! * [DLMF: Legendre-F¢s Integrals](https://dlmf.nist.gov/19.2#ii)-A +//! * [Wolfram MathWorld: Jacobi Amplitude](https://functions.wolfram.com/EllipticIntegrals/EllipticF/) +//! * [Wolfram MathWorld: Jacobi Amplitude](https://functions.wolfram.com/EllipticIntegrals/EllipticE/) +//! * [Wikipedia: Elliptic integral](https://en.wikipedia.org/wiki/Elliptic_integral) +//! +//! @groupheader{Example} +//! @godbolt{doc/ellint_fe.cpp} +//====================================================================================================================== + inline constexpr auto ellint_fe = eve::functor; +//====================================================================================================================== +//! @} +//====================================================================================================================== +} + +namespace kyosu::_ +{ + template + constexpr auto ellint_fe_(KYOSU_DELAY(), O const& o, Z u, M m) noexcept -> decltype(eve::zip(Z(), Z())) + { + if constexpr(!std::same_as, eve::underlying_type_t>) + return ellint_fe[o](u, eve::convert(m, as>())); + else + { + m = eve::abs(m); + if (O::contains(eve::modular)) m = eve::sin(m); + else if (O::contains(eve::eccentric)) m = eve::sqrt(m); + auto [phi, psi] = u; + + if (eve::all(is_real(u))) return eve::zip(kyosu::complex(eve::ellint_1[o](phi, m)), + kyosu::complex(eve::ellint_2[o](phi, m))); + auto m2 = eve::sqr(m); + auto thresh = eve::eps(eve::as(phi)); + phi = if_else(eve::abs(phi) < thresh, thresh, phi); //avoiding singularity at 0 + + auto b = -(eve::sqr(eve::cot(phi)) + m2*eve::sqr(eve::sinh(psi)*eve::csc(phi))-1+m2); //*eve::half(eve::as(phi)); + auto c = -(1-m2)*eve::sqr(eve::cot(phi)); + auto X1 = -b/2+eve::sqrt(eve::sqr(b)/4-c); + auto lambda = eve::acot(eve::sqrt(X1)); + auto mu = eve::atan( rec(m)*eve::sqrt(eve::dec(eve::sqr(tan(phi)*eve::cot(lambda))))); + + // change of variables taking into account periodicity ceil to the right + lambda = eve::sign_alternate(eve::floor(2*phi*eve::inv_pi(eve::as(phi))))*lambda+ + eve::pi(eve::as(phi))*eve::ceil(phi/eve::pi(eve::as(phi))-eve::half(eve::as(phi))+eve::eps(eve::as(phi))); + mu = eve::sign(psi)*mu; + auto mc = eve::sqrt(eve::oneminus(m2)); + lambda = if_else(is_real(u), phi, lambda); + auto f1 = eve::ellint_1[o](lambda, m); + auto e1 = eve::ellint_2[o](lambda, m); + auto f2 = eve::ellint_1[o](mu, mc); + auto e2 = eve::ellint_2[o](mu, mc); + f1 = if_else(is_imag(u), zero, f1); + e2 = if_else(is_eqz(mu), zero, e2); + + auto f = kyosu::complex(f1, f2); + + auto [sin_lam, cos_lam] = eve::sincos(lambda); + auto [sin_mu , cos_mu ] = eve::sincos(mu); + auto sin_mu2 = eve::sqr(sin_mu); + auto sin_lam2 = eve::sqr(sin_lam); + auto b1 = m2*sin_lam*cos_lam*sin_mu2*eve::sqrt(eve::oneminus(m2*sin_lam2)); + auto b2 = sin_mu*cos_mu*(1-m2*sin_lam2)*eve::sqrt(1-oneminus(m2)*sin_mu2); + auto b3 = eve::sqr(cos_mu) + m2*sin_lam2*sin_mu2; + auto e = kyosu::complex(b1, b2)/b3; + e += kyosu::complex(e1, if_else(is_real(u), zero, f2-e2)); + kyosu::real(e) = if_else(is_imag(u), zero, real(e)); + + auto test = eve::maxabs(phi, psi) < eve::eps(eve::as(kyosu::real(u))); + if (eve::any(test)) + { + f = if_else(test, u, f); + e = if_else(test, u, e); + } + return eve::zip(f, e); + } + } +} diff --git a/include/kyosu/functions/ellint_rc.hpp b/include/kyosu/functions/ellint_rc.hpp new file mode 100644 index 00000000..b18c4442 --- /dev/null +++ b/include/kyosu/functions/ellint_rc.hpp @@ -0,0 +1,143 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include +#include +#include +#include +#include +#include + +namespace kyosu +{ + + template + struct ellint_rc_t : eve::elementwise_callable + { + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b) const noexcept + { return eve::ellint_rc(a, b); } + + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b) const noexcept + requires(concepts::complex || concepts::complex) + { return KYOSU_CALL(a, b); } + + KYOSU_CALLABLE_OBJECT(ellint_rc_t, ellint_rc_); + }; + + +//================================================================================================ +//! @addtogroup functions +//! @{ +//! @var ellint_rc +//! @brief Computes the Carlson's elliptic integral +//! \f$ \mathbf{R}_\mathbf{C}(x, y) = +//! \frac12 \int_{0}^{\infty} \scriptstyle(t+x)^{-1/2}(t+y)]^{-1}\;\mathrm{d}t\f$. +//! +//! @groupheader{Header file} +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace eve +//! { +//! // Regular overload +//! constexpr auto ellint_rc(auto x, auto y) noexcept; // 1 +//! +//! // Lanes masking +//! constexpr auto ellint_rc[conditional_expr auto c](auto x, auto y) noexcept; // 2 +//! constexpr auto ellint_rc[logical_value auto m](auto x, auto y) noexcept; // 2 +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `x`, `y`: Can be a mix of complex and real floating values. `x`, `y` must be non zero +//! and have phase less in magnitude than \f$\pi\f$, with the exception that `x` may be 0. +//! * `c`: [Conditional expression](@ref eve::conditional_expr) masking the operation. +//! * `m`: [Logical value](@ref eve::logical_value) masking the operation. +//! +//! **Return value** +//! +//! 1. the value of the Carlson degenerate elliptic integral: +//! \f$\mathbf{R}_\mathbf{C}(x, y) = \frac12 \int_{0}^{\infty} +//! \scriptstyle(t+x)^{-1/2}(t+y)^{-1}\scriptstyle\;\mathrm{d}t\f$ is returned. +//! 2. [The operation is performed conditionally](@ref conditional) +//! +//! @groupheader{External references} +//! * [DLMF: Elliptic Integrals](https://dlmf.nist.gov/19.2) +//! * [Wolfram MathWorld: Elliptic Integral](https://mathworld.wolfram.com/CarlsonEllipticIntegrals.html) +//! +//! @groupheader{Example} +//! @godbolt{doc/ellint_rc.cpp} +//================================================================================================ + inline constexpr auto ellint_rc = eve::functor; +//================================================================================================ +//! @} +//================================================================================================ +} +namespace kyosu::_ +{ + template + constexpr auto ellint_rc_(KYOSU_DELAY(), O const& o, T xx, T yy) noexcept + { + using r_t = eve::underlying_type_t; + auto tol = [&](){ + if constexpr (O::contains(eve::threshold)) return o[eve::threshold].value(xx); + else return eve::eps(eve::as()); + }(); + auto internal_rc = [tol](auto x, auto y){ + T xn = x; + T yn = y; + T an = kyosu::average(x, y, y); + T a0 = an; + r_t epsi = eve::pow_abs(3 *tol, -eve::rec(r_t(8))); + auto q = epsi * kyosu::abs(an - xn); + r_t fn = eve::one(eve::as()); + + // duplication + auto hf = half(as(real(x))); + for( unsigned k = 1 ; k < 30; ++k ) + { + T root_xn = kyosu::sqrt(xn); + T root_yn = kyosu::sqrt(yn); + root_xn = if_else(real(root_xn) < 0, -root_xn, root_xn); + root_yn = if_else(real(root_yn) < 0, -root_yn, root_yn); + auto lambda = 2*root_xn*root_yn + yn; + an = kyosu::average(an, lambda) * hf; + xn = kyosu::average(xn, lambda) * hf; + yn = kyosu::average(yn, lambda) * hf; + q *= r_t(0.25); + fn *= r_t(4); + if( eve::all(q < kyosu::abs(an)) ) break; + } + auto s = (y-a0)/(fn*an); + // Taylor series expansion to the 7th order + return kyosu::horner(s, r_t(9.0/8.0), r_t(159.0/208.0), r_t(9.0/22.0), r_t(3.0/8.0), r_t(1.0/7.0), r_t(3.0/10.0), r_t(0.0), r_t(1.0))/kyosu::sqrt(an); + }; + auto test = eve::any(is_real(yy) && eve::is_ltz(kyosu::real(yy)));// compute Cauchy principal value on the cut + if (eve::any(test)) + { + auto fac = kyosu::sqrt(xx/(xx-yy)); + auto yr = if_else(test, yy, zero); + auto xr = if_else(test, xx, zero); + return kyosu::if_else(eve::is_lez(kyosu::real(yy)), internal_rc(xr-yr, -yr)*fac, internal_rc(xx, yy)); + } + else + return internal_rc(xx, yy); + } +} diff --git a/include/kyosu/functions/ellint_rd.hpp b/include/kyosu/functions/ellint_rd.hpp new file mode 100644 index 00000000..dd6aa2d8 --- /dev/null +++ b/include/kyosu/functions/ellint_rd.hpp @@ -0,0 +1,169 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include +#include +#include +#include + +namespace kyosu +{ + + template + struct ellint_rd_t : eve::elementwise_callable + { + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b, T2 c) const noexcept + { return eve::ellint_rd(a, b, c); } + + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b, T2 c) const noexcept + requires(concepts::complex || concepts::complex || concepts::complex) + { return KYOSU_CALL(a, b, c); } + + KYOSU_CALLABLE_OBJECT(ellint_rd_t, ellint_rd_); + }; + +//================================================================================================ +//! @addtogroup functions +//! @{ +//! @var ellint_rd +//! @brief Computes the Carlson's elliptic integral +//! \f$ \mathbf{R}_\mathbf{D}(x, y) = \frac32 \int_{0}^{\infty} +//! \scriptstyle[(t+x)(t+y)]^{-1/2}(t+z)^{-3/2}\;\mathrm{d}t\f$. +//! +//! @groupheader{Header file} +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace eve +//! { +//! // Regular overload +//! constexpr auto ellint_rd( auto x, auto y, auto z) noexcept; // 1 +//! +//! // semantic modifier +//! constexpr auto ellint_rd[threshold = tol](auto x, auto y, auto z) noexcept; // 1 +//! +//! // Lanes masking +//! constexpr auto ellint_rd[conditional_expr auto c](auto x, auto y, auto z) noexcept; // 2 +//! constexpr auto ellint_rd[logical_value auto m](auto x, auto y, auto z) noexcept; // 2 +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `x`, `y`, `z`: Can be a mix of complex and real floating values. `z` must be non zero +//! * `p`: [floating real arguments](@ref eve::floating_value). +//! * `c`: [Conditional expression](@ref eve::conditional_expr) masking the operation. +//! * `m`: [Logical value](@ref eve::logical_value) masking the operation. +//! +//! **Return value** +//! +//! 1. the value of the \f$\mathbf{R}_\mathbf{D}\f$ Carlson elliptic integral: +//! \f$ \frac32 \int_{0}^{\infty}\scriptstyle[(t+x)(t+y)]^{-1/2}(t+z)^{-3/2}\;\mathrm{d}t\f$ +//! is returned with relative error less in magnitude than tol +//! (tol default to [eps](@ref eve::eps)). +//! The integral is well defined if `x`, `y`, `z` lie in the +//! complex plane cut along the nonpositive real axis, +//! with the exception that at z must be non 0 +//! 2. [The operation is performed conditionnaly](@ref conditional) +//! +//! @groupheader{External references} +//! * [DLMF: Elliptic Integrals](https://dlmf.nist.gov/19.2) +//! * [Wolfram MathWorld: Elliptic Integral](https://mathworld.wolfram.com/CarlsonEllipticIntegrals.html) +//! +//! @groupheader{Example} +//! @godbolt{doc/ellint_rd.cpp} +//================================================================================================ + inline constexpr auto ellint_rd = eve::functor; +//================================================================================================ +//! @} +} +namespace kyosu::_ +{ + template + constexpr auto ellint_rd_(KYOSU_DELAY(), O const& o, T x, T y, T z) noexcept + { + using r_t = eve::underlying_type_t; + constexpr r_t hf = half(as()); + constexpr r_t qtr = r_t(0.25); + auto tol = [&](){ + if constexpr (O::contains(eve::threshold)) return o[eve::threshold].value(x); + else return eve::eps(eve::as(real(x))); + }(); + auto xn = x; + auto yn = y; + auto zn = z; + auto an = kyosu::average(x, y, z, z, z); + auto a0 = an; + auto epsi = eve::pow_abs(qtr*tol, -eve::rec(r_t(6))); + auto q = epsi *kyosu::maxabs(an - x, an - y, an - z)*r_t(1.2); + r_t fn(one(as())); + auto rd_sum(zero(as(x))); + + // duplication + for(unsigned k = 0 ; k < 30; ++k ) + { + T rx = kyosu::sqrt(xn); + T ry = kyosu::sqrt(yn); + T rz = kyosu::sqrt(zn); + rx = if_else(eve::is_ltz(real(rx)), -rx, rx); + ry = if_else(eve::is_ltz(real(ry)), -ry, ry); + rz = if_else(eve::is_ltz(real(rz)), -rz, rz); + T lambda = rx * ry + rx * rz + ry * rz; + rd_sum += fn / (rz * (zn + lambda)); + an = average(an, lambda) * hf; + xn = average(xn, lambda) * hf; + yn = average(yn, lambda) * hf; + zn = average(zn, lambda) * hf; + q *= qtr; + fn *= qtr; + if( eve::all(q < kyosu::abs(an)) ) break; + } + + T invan = kyosu::rec(an); + T xx = fn*(a0 - x)*invan; + T yy = fn*(a0 - y)*invan; + T zz = -(xx + yy)/3; + T xxyy = xx*yy; + T zz2 = sqr(zz); + T e2 = fnma(T(6), zz2, xxyy); + T e3 = (3*xxyy - 8*zz2)*zz; + T e4 = 3*(xxyy - zz2)*zz2; + T e5 = xxyy*zz2*zz; + + constexpr r_t c0 = r_t(-3/14.0); + constexpr r_t c1 = r_t(1/6.0); + constexpr r_t c2 = r_t(9/88.0); + constexpr r_t c3 = r_t(-3/22.0); + constexpr r_t c4 = r_t(-9/52.0); + constexpr r_t c5 = r_t(3/26.0); + constexpr r_t c6 = r_t(-1/16.0); + constexpr r_t c7 = r_t(3/40.0); + constexpr r_t c8 = r_t(3/20.0); + constexpr r_t c9 = r_t(45/272.0); + constexpr r_t c10 = r_t(-9/68.0); + + T e22 = sqr(e2); + T result = fn*invan*sqrt(invan) + *(1 + e2*c0 + e3*c1 + e22*c2 + e4*c3 + e2*(e3*c4 + e22*c6) + e5*c5 + + sqr(e3)*c7 + e2*e4*c8 + c9*e22*e3 + (e3*e4 + e2*e5)*c10); + + return r_t(3)*rd_sum+result; + } +} diff --git a/include/kyosu/functions/ellint_rf.hpp b/include/kyosu/functions/ellint_rf.hpp new file mode 100644 index 00000000..c643065e --- /dev/null +++ b/include/kyosu/functions/ellint_rf.hpp @@ -0,0 +1,151 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include +#include +#include + +namespace kyosu +{ + + template + struct ellint_rf_t : eve::elementwise_callable + { + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b, T2 c) const noexcept + { return eve::ellint_rf(a, b, c); } + + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b, T2 c) const noexcept + requires(concepts::complex || concepts::complex || concepts::complex) + { return KYOSU_CALL(a, b, c); } + + KYOSU_CALLABLE_OBJECT(ellint_rf_t, ellint_rf_); + }; + + +//================================================================================================ +//! @addtogroup functions +//! @{ +//! @var ellint_rf +//! @brief Computes the Carlson's elliptic integral +//! \f$ \mathbf{R}_\mathbf{F}(x, y) = +//! \frac32 \int_{0}^{\infty} \scriptstyle[(t+x)(t+y)(t+z)]^{-1/2}\scriptstyle\;\mathrm{d}t\f$. +//! +//! @groupheader{Header file} +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace eve +//! { +//! // Regular overload +//! constexpr auto ellint_rf(auto x, auto y, auto z) noexcept; // 1 +//! +//! // Lanes masking +//! constexpr auto ellint_rf[conditional_expr auto c](auto x, auto y, auto z) noexcept; // 2 +//! constexpr auto ellint_rf[logical_value auto m](auto x, auto y, auto z) noexcept; // 2 +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `x`, `y`, `z`: Can be a mix of complex and real floating values. +//! * `c`: [Conditional expression](@ref eve::conditional_expr) masking the operation. +//! * `m`: [Logical value](@ref eve::logical_value) masking the operation. +//! +//! **Return value** +//! +//! 1. the value of the \f$\mathbf{R}_\mathbf{F}\f$ Carlson elliptic integral: +//! \f$ \frac32 \int_{0}^{\infty} \scriptstyle[(t+x)(t+y)(t+z)]^{-1/2}\scriptstyle\;\mathrm{d}t\f$. +//! is returned. +//! The integral is well defined if `x`, `y`, `z` lie in the +//! complex plane cut along the nonpositive real axis, +//! with the exception that one of `x`, `y`, `z` must be non 0 +//! 2. [The operation is performed conditionnaly](@ref conditional) +//! +//! @groupheader{External references} +//! * [DLMF: Elliptic Integral](https://dlmf.nist.gov/19.2) +//! * [Wolfram MathWorld: Elliptic Integral](https://mathworld.wolfram.com/CarlsonEllipticIntegrals.html) +//! +//! @groupheader{Example} +//! @godbolt{doc/elliptic/ellint_rf.cpp} +//================================================================================================ + inline constexpr auto ellint_rf = eve::functor; +//================================================================================================ +//! @} +//================================================================================================ +} +namespace kyosu::_ +{ + template + constexpr auto ellint_rf_(KYOSU_DELAY(), O const& o, T x, T y, T z) noexcept + { + using r_t = eve::underlying_type_t; + auto tol = [&](){ + if constexpr (O::contains(eve::threshold)) return o[eve::threshold].value(x); + else return eve::eps(eve::as()); + }(); + T xn = x; + T yn = y; + T zn = z; + T an = kyosu::average(x, y, z); + T a0 = an; + r_t epsi = eve::pow_abs(3 *tol, -eve::rec(r_t(6))); + auto q = epsi * kyosu::maxabs(an - xn, an - yn, an - zn); + r_t fn = eve::one(eve::as()); + + // duplication + r_t hf = half(as()); + for(unsigned k = 1 ; k < 30; ++k ) + { + auto root_x = kyosu::sqrt(xn); + auto root_y = kyosu::sqrt(yn); + auto root_z = kyosu::sqrt(zn); + root_x = if_else(real(root_x) < 0, -root_x, root_x); + root_y = if_else(real(root_y) < 0, -root_y, root_y); + root_z = if_else(real(root_z) < 0, -root_z, root_z); + auto lambda = root_x * root_y + root_x * root_z + root_y * root_z; + an = kyosu::average(an, lambda) * hf; + xn = kyosu::average(xn, lambda) * hf; + yn = kyosu::average(yn, lambda) * hf; + zn = kyosu::average(zn, lambda) * hf; + q *= r_t(0.25); + fn *= r_t(4); + if( eve::all(q < kyosu::abs(an)) ) break; + } + auto denom = kyosu::rec(an * fn); + auto xx = (a0 - x) * denom; + auto yy = (a0 - y) * denom; + auto zz = -xx - yy; + + // Taylor series expansion to the 7th order + auto p = xx * yy; + auto e2 = kyosu::fnma(zz, zz, p); + auto e3 = p * zz; + constexpr r_t c0 = r_t(1/14.0); + constexpr r_t c1 = r_t(3/104.0); + constexpr r_t c2 = r_t(-1/10.0); + constexpr r_t c4 = r_t(1/24.0); + constexpr r_t c5 = r_t(-3/44.0); + constexpr r_t c6 = r_t(-5/208.0); + constexpr r_t c7 = r_t(-1/16.0); + return (kyosu::fma(e3, + kyosu::fma(e3, c1, c0), + kyosu::fma(e2, (c2 + e3 * c5 + e2 * (c4 + e2 * c6 + e3 * c7)), kyosu::one(as(x))))) + / kyosu::sqrt(an); + } +} diff --git a/include/kyosu/functions/ellint_rg.hpp b/include/kyosu/functions/ellint_rg.hpp new file mode 100644 index 00000000..66b2f2df --- /dev/null +++ b/include/kyosu/functions/ellint_rg.hpp @@ -0,0 +1,109 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include +#include + +namespace kyosu +{ + + template + struct ellint_rg_t : eve::elementwise_callable + { + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b, T2 c) const noexcept + { return eve::ellint_rg(a, b, c); } + + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b, T2 c) const noexcept + requires(concepts::complex || concepts::complex || concepts::complex) + { return KYOSU_CALL(a, b, c); } + + KYOSU_CALLABLE_OBJECT(ellint_rg_t, ellint_rg_); + }; + +//================================================================================================ +//! @addtogroup functions +//! @{ +//! @var ellint_rg +//! @brief Computes the Carlson's elliptic integral +//! \f$ \mathbf{R}_\mathbf{G}(x, y) = \frac1{4\pi} \int_{0}^{2\pi}\int_{0}^{\pi} +//! \scriptstyle\sqrt{x\sin^2\theta\cos^2\phi +//! +y\sin^2\theta\sin^2\phi +//! +z\cos^2\theta} \scriptstyle\;\mathrm{d}\theta\;\mathrm{d}\phi\f$ +//! +//! @groupheader{Header file} +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace eve +//! { +//! // Regular overload +//! constexpr auto ellint_rg(auto x, auto y, auto z) noexcept; // 1 +//! +//! // Lanes masking +//! constexpr auto ellint_rg[conditional_expr auto c](/*all previous overloads*/) noexcept; // 2 +//! constexpr auto ellint_rg[logical_value auto m](/*all previous overloads*/) noexcept; // 2 +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `x`, `y`, `z`: complex or real arguments. +//! * `c`: [Conditional expression](@ref eve::conditional_expr) masking the operation. +//! * `m`: [Logical value](@ref eve::logical_value) masking the operation. +//! +//! **Return value** +//! +//! 1. the value of the \f$\mathbf{R}_\mathbf{G}\f$ Carlson elliptic integral: +//! \f$\frac14 \int_{0}^{\infty}[(t+x)(t+y)(t+z)]^{-1/2}\left(\frac{x}{t+x}+\frac{y}{t+y}+\frac{z}{t+z}\right)t\;\mathrm{d}t\f$ +//! is returned. +//! All of x, y, z may be 0 and those that are nonzero must lie in the complex plane cut +//! along the nonpositive real axis +//! 2. [The operation is performed conditionnaly](@ref conditional) +//! +//! @groupheader{External references} +//! * [DLMF: Elliptic Integrals](https://dlmf.nist.gov/19.2) +//! * [Wolfram MathWorld: Elliptic Integral](https://mathworld.wolfram.com/CarlsonEllipticIntegrals.html) +//! +//! +//! @groupheader{Example} +//! @godbolt{doc/ellint_rg.cpp} +//================================================================================================ + inline constexpr auto ellint_rg = eve::functor; +//================================================================================================ +//! @} +} +namespace kyosu::_ +{ + template + constexpr auto ellint_rg_(KYOSU_DELAY(), O const& o, T x, T y, T z) noexcept + { + auto ax = kyosu::abs(x); + auto ay = kyosu::abs(y); + auto az = kyosu::abs(z); + using r_t = eve::underlying_type_t; + eve::swap_if(ax < ay, x, y); + eve::swap_if(ax < az, x, z); + eve::swap_if(ay > az, y, z); + // now all(x >= z) and all(z >= y) + auto root = kyosu::sqrt(x*y/z); + root = if_else(eve::is_ltz(imag(root)), -root, root); + return average(z*ellint_rf[o](x,y,z)-(x-z)*(y-z)*ellint_rd[o](x, y, z)*eve::third(as()), root); + } +} diff --git a/include/kyosu/functions/ellint_rj.hpp b/include/kyosu/functions/ellint_rj.hpp new file mode 100644 index 00000000..f4bb9f81 --- /dev/null +++ b/include/kyosu/functions/ellint_rj.hpp @@ -0,0 +1,168 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include +#include +#include +#include +#include + +namespace kyosu +{ + + template + struct ellint_rj_t : eve::elementwise_callable + { + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b, T2 c, T3 d) const noexcept + { return eve::ellint_rj(a, b, c, d); } + + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b, T2 c, T3 d) const noexcept + requires(concepts::complex || concepts::complex || concepts::complex|| concepts::complex) + { return KYOSU_CALL(a, b, c, d); } + + KYOSU_CALLABLE_OBJECT(ellint_rj_t, ellint_rj_); + }; + +//================================================================================================ +//! @addtogroup functions +//! @{ +//! @var ellint_rj +//! @brief Computes the Carlson's elliptic integral +//! \f$ \mathbf{R}_\mathbf{J}(x, y) = \frac32 \int_{0}^{\infty} +//! \scriptstyle(t+p)^{-1}[(t+x)(t+y)(t+z)]^{-1/2}\scriptstyle\;\mathrm{d}t\f$. +//! +//! @groupheader{Header file} +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace eve +//! { +//! // Regular overload +//! constexpr auto ellint_rj(auto x, auto y, auto z, auto p) noexcept; // 1 +//! +//! // semantic modifier +//! constexpr auto ellint_rj[threshold = tol](auto x, auto y, auto z, auto p) noexcept; // 1 +//! +//! // Lanes masking +//! constexpr auto ellint_rj[conditional_expr auto c](auto x, auto y, auto z, auto p) noexcept; // 2 +//! constexpr auto ellint_rj[logical_value auto m](auto x, auto y, auto z, auto p) noexcept; // 2 +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `x`, `y`, `z`: [floating real arguments](@ref eve::floating_value). +//! * `p`: [floating real arguments](@ref eve::floating_value). +//! * `c`: [Conditional expression](@ref eve::conditional_expr) masking the operation. +//! * `m`: [Logical value](@ref eve::logical_value) masking the operation. +//! +//! **Return value** +//! +//! 1. the value of the \f$\mathbf{R}_\mathbf{J}\f$ Carlson elliptic integral: +//! \f$ \frac32 \int_{0}^{\infty} +//! \scriptstyle(t+p)^{-1}[(t+x)(t+y)(t+z)]^{-1/2}\;\mathrm{d}t\f$ is returned with relative error less in magnitude than tol +//! (default to [eps](@ref eve::eps)), +//! The integral is well defined if `x`, `y`, `z` lie in the +//! complex plane cut along the nonpositive real axis, +//! with the exception that at at most one of `x`, `y`, `z` can be 0. +//! 2. [The operation is performed conditionnaly](@ref conditional) +//! +//! @groupheader{External references} +//! * [DLMF: Elliptic Integrals](https://dlmf.nist.gov/19.2) +//! * [Wolfram MathWorld: Elliptic Integral](https://mathworld.wolfram.com/CarlsonEllipticIntegrals.html) +//! +//! @groupheader{Example} +//! @godbolt{doc/elliptic/ellint_rj.cpp} +//================================================================================================ + inline constexpr auto ellint_rj = eve::functor; +//================================================================================================ +//! @} +} +namespace kyosu::_ +{ + template + constexpr auto ellint_rj_(KYOSU_DELAY(), O const& o, T x, T y, T z, T p) noexcept + { + using r_t = eve::underlying_type_t; + auto tol = [&](){ + if constexpr (O::contains(eve::threshold)) return o[eve::threshold].value(x); + else return eve::eps(eve::as(real(x))); + }(); + auto xn = x; + auto yn = y; + auto zn = z; + auto pn = p; + + auto an = kyosu::average(x, y, z, p, p); + auto a0 = an; + auto delta = (p - x) * (p - y) * (p - z); + auto q = eve::pow_abs(tol/4, -r_t(1)/8)*kyosu::maxabs(an - x, an - y, an - z, an - p); + + r_t fmn(one(as())); // 4^-n + auto rc_sum(zero(as(x))); + + for( unsigned n = 0; n < 30; ++n ) + { + auto rx = kyosu::sqrt(xn); + auto ry = kyosu::sqrt(yn); + auto rz = kyosu::sqrt(zn); + auto rp = kyosu::sqrt(pn); + rx = if_else(eve::is_ltz(real(rx)), -rx, rx); + ry = if_else(eve::is_ltz(real(ry)), -ry, ry); + rz = if_else(eve::is_ltz(real(rz)), -rz, rz); + rp = if_else(eve::is_ltz(real(rp)), -rp, rp); + auto dn = (rp + rx) * (rp + ry) * (rp + rz); + auto en = delta / dn; + en /= dn; + auto lambda = fma(rx, ry, fma(rx, rz, ry * rz)); + rc_sum += fmn / dn * ellint_rc(T(1), inc(en)); + + an = (an + lambda) * r_t(0.25); // / 4; + fmn /= 4; + + if( eve::all(fmn * q < kyosu::abs(an)) ) break; + + xn = (xn + lambda) * T(0.25); // / 4; + yn = (yn + lambda) * T(0.25); // / 4; + zn = (zn + lambda) * T(0.25); // / 4; + pn = (pn + lambda) * T(0.25); // / 4; + delta *= T(0.015625); // /= 64; + } + auto fmninvan = fmn * kyosu::rec(an); + auto xx = (a0 - x) * fmninvan; + auto yy = (a0 - y) * fmninvan; + auto zz = (a0 - z) * fmninvan; + p = -(xx + yy + zz) / 2; + auto p2 = sqr(p); + auto xxyy = xx * yy; + auto e2 = xxyy + kyosu::fma(xx, zz, kyosu::fms(yy, zz, 3 * p2)); + auto p3 = p * p2; + auto xxyyzz = xxyy * zz; + auto e3 = xxyyzz + 2 * e2 * p + 4 * p3; + auto e4 = (2 * xxyyzz + e2 * p + 3 * p3) * p; + auto e5 = xxyyzz * p2; + auto e22 = sqr(e2); + auto result = fmninvan /kyosu::sqrt(an) + * (1 - 3 * e2 / 14 + e3 / 6 + 9 * e22 / 88 - 3 * e4 / 22 - 9 * e2 * e3 / 52 + + 3 * e5 / 26 - e2 * e22 / 16 + 3 * e3 * e3 / 40 + 3 * e2 * e4 / 20 + + 45 * e22 * e3 / 272 - 9 * (e3 * e4 + e2 * e5) / 68); + result += 6 * rc_sum; + return result; + } +} diff --git a/include/kyosu/functions/fma.hpp b/include/kyosu/functions/fma.hpp index 30f3cdc7..7287299f 100644 --- a/include/kyosu/functions/fma.hpp +++ b/include/kyosu/functions/fma.hpp @@ -12,7 +12,7 @@ namespace kyosu { template - struct fma_t : eve::strict_elementwise_callable + struct fma_t : eve::elementwise_callable { template requires(concepts::cayley_dickson || concepts::cayley_dickson || concepts::cayley_dickson) diff --git a/include/kyosu/functions/fnma.hpp b/include/kyosu/functions/fnma.hpp index 25b786e8..242a34b8 100644 --- a/include/kyosu/functions/fnma.hpp +++ b/include/kyosu/functions/fnma.hpp @@ -12,7 +12,7 @@ namespace kyosu { template - struct fnma_t : eve::strict_elementwise_callable + struct fnma_t : eve::elementwise_callable { template requires(concepts::cayley_dickson || concepts::cayley_dickson || concepts::cayley_dickson) diff --git a/include/kyosu/functions/jacobi_elliptic.hpp b/include/kyosu/functions/jacobi_elliptic.hpp new file mode 100644 index 00000000..bcba129d --- /dev/null +++ b/include/kyosu/functions/jacobi_elliptic.hpp @@ -0,0 +1,125 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright: KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#pragma once +#include +#include + +namespace kyosu +{ + + template + struct jacobi_elliptic_t : eve::strict_elementwise_callable + { + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b) const noexcept + { return (*this)(kyosu::complex(a), b); } + + template + constexpr KYOSU_FORCEINLINE + auto operator()(T0 a, T1 b) const noexcept + { return KYOSU_CALL(a, b); } + + + KYOSU_CALLABLE_OBJECT(jacobi_elliptic_t, jacobi_elliptic_); + }; + +//====================================================================================================================== +//! @addtogroup functions +//! @{ +//! @var jacobi_elliptic +//! @brief Computes Jacobi's Amplitude function. +//! +//! @code +//! #include +//! @endcode +//! +//! @groupheader{Callable Signatures} +//! +//! @code +//! namespace kyosu +//! { +//! // Regular overload +//! template constexpr Z jacobi_elliptic(real m, real z) noexcept; //1 +//! template constexpr auto jacobi_elliptic(real m, complex z) noexcept; //1 +//! +//! //Semantic modifiers +//! template constexpr Z jacobi_elliptic[modular](real alpha, real z) noexcept; //1 +//! template constexpr auto jacobi_elliptic[modular](real alpha, complex z) noexcept; //1 +//! template constexpr Z jacobi_elliptic[eccentric](real k, real z) noexcept; //1 +//! template constexpr auto jacobi_elliptic[eccentric](real k, complex z) noexcept; //1 +//! template constexpr auto jacobi_elliptic[threshold = tol](real t, real z) noexcept; //1 +//! template constexpr auto jacobi_elliptic[threshold = tol](real t, complex z) noexcept; //1 +//! } +//! @endcode +//! +//! **Parameters** +//! +//! * `u`: argument. +//! * `m`: amplitude parameter (\f$0\le m\le 1). +//! * `alpha `: modular angle in radian. +//! * `k`: elliptic modulus (eccentricity) . +//! * `tol': accuracy tolerance (by defaut [epsilon](@ref eve::epsilon). +//! * `c`: [Conditional expression](@ref eve::conditional_expr) masking the operation. +//! * `l`: [Logical value](@ref eve::logical_value) masking the operation. +//! +//! +//! **Return value** +//! +//! 1. return the jacobian amplitude function. Take care that the meaning of the second parameters +//! depends on the option used (see note below). +//! +//! @note +//! * \f$\alpha\f$ is named the modular angle given in radian (modular option). +//! * \f$ k = \sin\alpha \f$ is named the elliptic modulus or eccentricity (eccentric option). +//! * \f$ m = k^2 = \sin^2\alpha\f$ is named the parameter (no option). +//! Each of the above three quantities is completely determined by any of the others (given that they are non-negative). +//! Thus, they can be used interchangeably up to roundings erreors by giving the right option. +//! +//! @groupheader{External references} +//! * [C++ standard reference: am](https://en.cppreference.com/w/cpp/numeric/special_functions/am) +//! * [DLMF: Jacobi Amplitude](https://dlmf.nist.gov/22.16) +//! * [Wolfram MathWorld: Jacobi Amplitude](https://mathworld.wolfram.com/JacobiAmplitude.html) +//! * [Wikipedia: Jacobi elliptic functions](https://en.wikipedia.org/wiki/Jacobi_elliptic_functions) +//! +//! @groupheader{Example} +//! @godbolt{doc/jacobi_elliptic.cpp} +//====================================================================================================================== + inline constexpr auto jacobi_elliptic = eve::functor; +//====================================================================================================================== +//! @} +//====================================================================================================================== +} + +namespace kyosu::_ +{ + template + KYOSU_FORCEINLINE constexpr auto jacobi_elliptic_(KYOSU_DELAY(), O const& o, Z u, M m) noexcept + { + auto tol = [&](){ + if constexpr (O::contains(eve::threshold)) return o[eve::threshold].value(m); + else return eve::epsilon(eve::maximum(eve::abs(m))); + }(); + m = eve::abs(m); + if (O::contains(eve::modular)) m = eve::sin(m); + else if (O::contains(eve::eccentric)) m = eve::sqrt(m); + auto [phi, psi] = u; + auto [s,c,d] = eve::jacobi_elliptic[eve::threshold = tol](phi, m); + if (eve::all(is_real(u))) return eve::zip(kyosu::complex(s), kyosu::complex(c), kyosu::complex(d)); + auto m2 = eve::sqr(m); + auto mc = eve::sqrt(eve::oneminus(m2)); + auto [s1,c1,d1] = eve::jacobi_elliptic[eve::threshold = tol](psi,mc); + auto idelta = kyosu::rec(eve::fam(eve::sqr(c1), m2, eve::sqr(s*s1))); + auto ds1 = d*s1; + auto sn = complex(s*d1, c*ds1*c1)*idelta; + auto cn = complex(c*c1, -s*ds1*d1)*idelta; + auto dn = complex(d*c1*d1, -m2*s*c*s1)*idelta; + return eve::zip(sn, cn, dn); + } +} diff --git a/include/kyosu/functions/oneminus.hpp b/include/kyosu/functions/oneminus.hpp index a1c9657d..480df018 100644 --- a/include/kyosu/functions/oneminus.hpp +++ b/include/kyosu/functions/oneminus.hpp @@ -9,7 +9,6 @@ #include #include #include -#include namespace kyosu { diff --git a/include/kyosu/functions/reverse_horner.hpp b/include/kyosu/functions/reverse_horner.hpp index f3e47dfe..6d2d7051 100644 --- a/include/kyosu/functions/reverse_horner.hpp +++ b/include/kyosu/functions/reverse_horner.hpp @@ -92,7 +92,7 @@ namespace kyosu //! //! * `coefs...` : real or cayley-dickson arguments. The coefficients by increasing power order //! -//! * `tup` : kumi tuple containing The coefficients by decreasing power order. +//! * `tup` : kumi tuple containing The coefficients by increasing power order. //! //! **Return value** //! diff --git a/include/kyosu/functions/sqrt.hpp b/include/kyosu/functions/sqrt.hpp index 6f4e28d2..99793cbe 100644 --- a/include/kyosu/functions/sqrt.hpp +++ b/include/kyosu/functions/sqrt.hpp @@ -8,6 +8,7 @@ #pragma once #include #include +#include namespace kyosu { @@ -120,7 +121,7 @@ namespace kyosu::_ , Z(rr1, ii1) , Z(ii1, rr1) ); - res = if_else(is_pure(z), eve::sqrt_2(eve::as(r))*Z( eve::half(eve::as(r)), eve::half(eve::as(r)))*eve::sqrt(iz), res); + res = if_else(is_pure(z), eve::sqrt_2(eve::as(r))*Z( eve::half(eve::as(r)), eve::half(eve::as(r)))*eve::sqrt(y), res); if (eve::any(is_not_finite(z))) { res = kyosu::if_else(rz == eve::minf(eve::as(rz)) diff --git a/test/doc/ellint_fe.cpp b/test/doc/ellint_fe.cpp new file mode 100644 index 00000000..7c92ae7e --- /dev/null +++ b/test/doc/ellint_fe.cpp @@ -0,0 +1,18 @@ +#include +#include +#include +#include + +int main() +{ + std::cout<< std::setprecision(16) << std::endl; + using w_t = eve::wide>; + auto z = kyosu::complex(w_t(2.0, 1.0, 1.0, 0.0 ), w_t(1.0, 1.0, 0.0, 1.0)); + auto m = w_t(0.5); + auto [f, e] = kyosu::ellint_fe(z, eve::sqrt(m)); + std::cout << " z " << z << std::endl; + std::cout << " m " << m << std::endl; + std::cout << "f " << f << std::endl; + std::cout << "e " << e << std::endl; + return 0; +} diff --git a/test/doc/ellint_rc.cpp b/test/doc/ellint_rc.cpp new file mode 100644 index 00000000..520c6fb6 --- /dev/null +++ b/test/doc/ellint_rc.cpp @@ -0,0 +1,21 @@ +// revision 1 +#include +#include + +using wide_t = eve::wide >; +using r_t = double; +wide_t ref1 = { 3.0, 2.0, 1.0, 0.5}; +wide_t imf1 = { 2.0, -1.0, -5.0, 0.0}; +wide_t ref2 = { 0.0, 1.0, 2.0, 3.0}; +wide_t imf2 = { 1.0f , -4.0, -2.0, 0.0}; +auto pf = kyosu::complex_t(ref1, imf1); +auto qf = kyosu::complex_t(ref2, imf2); + +int main() +{ + std::cout << "<- pf = " << pf << "\n"; + std::cout << "<- qf = " << qf << "\n"; + + std::cout << "-> ellint_rc(pf, qf) = " << kyosu::ellint_rc(pf, qf) << "\n"; + +} diff --git a/test/doc/ellint_rd.cpp b/test/doc/ellint_rd.cpp new file mode 100644 index 00000000..a703660a --- /dev/null +++ b/test/doc/ellint_rd.cpp @@ -0,0 +1,26 @@ +// revision 1 +#include +#include + +using wide_t = eve::wide >; +using r_t = double; +wide_t re1 = { 3.0, 2.0, 1.0, 0.5}; +wide_t im1 = { 2.0, -1.0, -5.0, 0.0}; +wide_t re2 = { 0.0, 1.0, 2.0, 3.0}; +wide_t im2 = { 1.0 , -4.0, -2.0, 0.0}; +wide_t re3 = { 0.1, -1.0, 2.0, 4.0}; +wide_t im3 = { 2.0 , -4.0, -3.0, 0.0}; +auto p = kyosu::complex_t(re1, im1); +auto q = kyosu::complex_t(re2, im2); +auto r = kyosu::complex_t(re3, im3); + +int main() +{ + std::cout << "<- p = " << p << "\n"; + std::cout << "<- q = " << q << "\n"; + std::cout << "<- r = " << r << "\n"; + + std::cout << "-> ellint_rd(p, q, r) = " << kyosu::ellint_rd(p, q, r) << "\n"; + std::cout << "-> ellint_rd[ignore_last(2)](p, q, r)= " << kyosu::ellint_rd[eve::ignore_last(2)](p, q, r) << "\n"; + std::cout << "-> ellint_rd[q != 4.0](p, q, r) = " << kyosu::ellint_rd[q != 4.0](p, q, r) << "\n"; +} diff --git a/test/doc/ellint_rf.cpp b/test/doc/ellint_rf.cpp new file mode 100644 index 00000000..37851bc7 --- /dev/null +++ b/test/doc/ellint_rf.cpp @@ -0,0 +1,26 @@ +// revision 1 +#include +#include + +using wide_t = eve::wide >; +using r_t = double; +wide_t re1 = { 3.0, 2.0, 1.0, 0.5}; +wide_t im1 = { 2.0, -1.0, -5.0, 0.0}; +wide_t re2 = { 0.0, 1.0, 2.0, 3.0}; +wide_t im2 = { 1.0 , -4.0, -2.0, 0.0}; +wide_t re3 = { 0.1, -1.0, 2.0, 4.0}; +wide_t im3 = { 2.0 , -4.0, -3.0, 0.0}; +auto p = kyosu::complex_t(re1, im1); +auto q = kyosu::complex_t(re2, im2); +auto r = kyosu::complex_t(re3, im3); + +int main() +{ + std::cout << "<- p = " << p << "\n"; + std::cout << "<- q = " << q << "\n"; + std::cout << "<- r = " << r << "\n"; + + std::cout << "-> ellint_rf(p, q, r) = " << kyosu::ellint_rf(p, q, r) << "\n"; + std::cout << "-> ellint_rf[ignore_last(2)](p, q, r)= " << kyosu::ellint_rf[eve::ignore_last(2)](p, q, r) << "\n"; + std::cout << "-> ellint_rf[q != 4.0](p, q, r) = " << kyosu::ellint_rf[q != 4.0](p, q, r) << "\n"; +} diff --git a/test/doc/ellint_rg.cpp b/test/doc/ellint_rg.cpp new file mode 100644 index 00000000..b6588795 --- /dev/null +++ b/test/doc/ellint_rg.cpp @@ -0,0 +1,26 @@ +// revision 1 +#include +#include + +using wide_t = eve::wide >; +using r_t = double; +wide_t re1 = { 3.0, 2.0, 1.0, 0.5}; +wide_t im1 = { 2.0, -1.0, -5.0, 0.0}; +wide_t re2 = { 0.0, 1.0, 4.0, 3.0}; +wide_t im2 = { 1.0 , -4.0, 0.0, 0.0}; +wide_t re3 = { 0.1, -1.0, 2.0, 4.0}; +wide_t im3 = { 2.0 , -4.0, -3.0, 0.0}; +auto p = kyosu::complex_t(re1, im1); +auto q = kyosu::complex_t(re2, im2); +auto r = kyosu::complex_t(re3, im3); + +int main() +{ + std::cout << "<- p = " << p << "\n"; + std::cout << "<- q = " << q << "\n"; + std::cout << "<- r = " << r << "\n"; + + std::cout << "-> ellint_rg(p, q, r) = " << kyosu::ellint_rg(p, q, r) << "\n"; + std::cout << "-> ellint_rg[ignore_last(2)](p, q, r)= " << kyosu::ellint_rg[eve::ignore_last(2)](p, q, r) << "\n"; + std::cout << "-> ellint_rg[q != 4.0](p, q, r) = " << kyosu::ellint_rg[q != 4.0](p, q, r) << "\n"; +} diff --git a/test/doc/ellint_rj.cpp b/test/doc/ellint_rj.cpp new file mode 100644 index 00000000..49dbd2ea --- /dev/null +++ b/test/doc/ellint_rj.cpp @@ -0,0 +1,31 @@ +// revision 1 +#include +#include + +using wide_t = eve::wide >; +using r_t = double; +wide_t re1 = { 3.0, 2.0, 1.0, 0.5}; +wide_t im1 = { 2.0, -1.0, -5.0, 0.0}; +wide_t re2 = { 0.0, 1.0, 2.0, 3.0}; +wide_t im2 = { 1.0 , -4.0, -2.0, 0.0}; +wide_t re3 = { 0.1, -1.0, 2.0, 4.0}; +wide_t im3 = { 2.0 , -4.0, -3.0, 0.0}; +wide_t re4 = { 0.1, -1.5, 2.3, 0.5}; +wide_t im4 = { 0.0 , -1.0, -2.0, 0.0}; +auto p = kyosu::complex_t(re1, im1); +auto q = kyosu::complex_t(re2, im2); +auto r = kyosu::complex_t(re3, im3); +auto s = kyosu::complex_t(re4, im4); + +int main() +{ + std::cout << "<- p = " << p << "\n"; + std::cout << "<- q = " << q << "\n"; + std::cout << "<- r = " << r << "\n"; + std::cout << "<- s = " << s << "\n"; + + + std::cout << "-> ellint_rj(p, q, r, s) = " << kyosu::ellint_rj(p, q, r, s) << "\n"; + std::cout << "-> ellint_rj[ignore_last(2)](p, q, r, s)= " << kyosu::ellint_rj[eve::ignore_last(2)](p, q, r, s) << "\n"; + std::cout << "-> ellint_rj[q != 4.0](p, q, r, s) = " << kyosu::ellint_rj[q != 4.0](p, q, r, s) << "\n"; +} diff --git a/test/doc/jacobi_elliptic.cpp b/test/doc/jacobi_elliptic.cpp new file mode 100644 index 00000000..1a93a24e --- /dev/null +++ b/test/doc/jacobi_elliptic.cpp @@ -0,0 +1,19 @@ +#include +#include +#include +#include + +int main() +{ + std::cout<< std::setprecision(16) << std::endl; + using w_t = eve::wide>; + auto z = kyosu::complex(w_t(2.0, 1.0), w_t(0.0, 1.0)); + auto m = w_t(0.5, 0.7); + auto [sn, cn, dn] = kyosu::jacobi_elliptic(z, m); + std::cout << " z " << z << std::endl; + std::cout << " m " << m << std::endl; + std::cout << "sn(z, m) " << sn << std::endl; + std::cout << "cn(z, m) " << cn << std::endl; + std::cout << "dn(z, m) " << dn << std::endl; + return 0; +} diff --git a/test/unit/function/am.cpp b/test/unit/function/am.cpp new file mode 100755 index 00000000..4da0bbbf --- /dev/null +++ b/test/unit/function/am.cpp @@ -0,0 +1,31 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include +#include + +TTS_CASE_TPL ( "Check jacobi_elliptic over complex" + , kyosu::scalar_real_types + ) + (tts::type) +{ + auto pr = tts::prec(4.0e-2, 1.0e-8); + using c_t = kyosu::complex_t; + using a_t = std::array; + a_t phi= {0.000000000000000e+00, 1.346396851538483e+00, 2.692793703076966e+00, 4.039190554615448e+00, 5.385587406153931e+00, 6.731984257692414e+00, 8.078381109230897e+00, }; + a_t m= {0.000000000000000e+00, 1.428571428571428e-01, 2.857142857142857e-01, 4.285714285714285e-01, 5.714285714285714e-01, 7.142857142857143e-01, 8.571428571428572e-01, }; + a_t ramp = { 0, 1.696717861493621, 0.853380674698701, 5.842387616977985, 4.900661768634335, 4.798157643305010, 5.063923053256403}; + a_t iamp = {8.078381109097574, -1.545860408368032, -1.126101384553287, -0.192035497540286, 0.569586142084441, 0.611135912286602, 0}; + + for(int i=0; i < 6; ++i) + { + c_t z(phi[i], phi[6-i]); + auto amp = kyosu::am(z, eve::sqrt(m[i])); + TTS_RELATIVE_EQUAL(amp, kyosu::complex(ramp[i], iamp[i]), pr); + } +}; diff --git a/test/unit/function/ellint_fe.cpp b/test/unit/function/ellint_fe.cpp new file mode 100644 index 00000000..626283cc --- /dev/null +++ b/test/unit/function/ellint_fe.cpp @@ -0,0 +1,34 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include +#include + +TTS_CASE_TPL ( "Check ellint_fe over complex" + , kyosu::scalar_real_types + ) + (tts::type) +{ + auto pr = tts::prec(4.0e-2, 1.0e-8); + using c_t = kyosu::complex_t; + using a_t = std::array; + constexpr auto nan = eve::nan(eve::as()); + a_t phi= {0.000000000000000e+00, 1.346396851538483e+00, 2.692793703076966e+00, 4.039190554615448e+00, 5.385587406153931e+00, 6.731984257692414e+00, 8.078381109230897e+00, }; + a_t m= {0.000000000000000e+00, 1.428571428571428e-01, 2.857142857142857e-01, 4.285714285714285e-01, 5.714285714285714e-01, 7.142857142857143e-01, 8.571428571428572e-01, }; + a_t ree = {nan, 1.545486188403233e+02, -2.240669689136373e+01, 1.731898389351810e+01, 8.395951912807700e-01, 5.720419655812918e+00, 5.779502973552437e+00}; + a_t ime = {nan, 3.654685952505582e+01, 5.340773457856197e+01, 1.218092118477447e+01, 3.853460993937446e+00, 1.539669180714825e+00, 0.0}; + a_t ref = {0.0, 6.150211481292479e-03, 3.403816928620261e+00, 3.637799767891663e+00, 7.539639717553287e+00, 8.632555723502625e+00, 1.262536872229223e+01}; + a_t imf = {nan, 2.410171764691075e+00, 2.081529547642794e+00, 1.886274604714848e+00, 1.685592709331894e+00, 1.157716506821809e+00, 0.0}; + for(int i=0; i <= 6; ++i) + { + c_t z(phi[i], phi[6-i]); + auto [f, e] = kyosu::ellint_fe(z, eve::sqrt(m[i])); + TTS_RELATIVE_EQUAL(f, kyosu::complex(ref[i], imf[i]), pr); + TTS_RELATIVE_EQUAL(e, kyosu::complex(ree[i], ime[i]), pr); + } +}; diff --git a/test/unit/function/ellint_rc.cpp b/test/unit/function/ellint_rc.cpp new file mode 100644 index 00000000..2773b7fc --- /dev/null +++ b/test/unit/function/ellint_rc.cpp @@ -0,0 +1,37 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include +#include +// the values tested are those contained in the original Carlson article +// computation of elliptic integrals, Bille C. Carlson arXiv:math/9409227 + +TTS_CASE_TPL ( "Check ellint_rf over complex" + , kyosu::scalar_real_types + ) + (tts::type) +{ + auto pr = tts::prec(4.0e-3, 1.0e-8); + using c_t = kyosu::complex_t; + using a_t = std::array; + a_t rx= { 0.0, 2.25, 0.0, 0.0, 0.25, 0.0}; + a_t ix= { 0.0, 0.0, 0.0, -1.0, 0.0, 1.0}; + a_t ry= { 0.25, 2.0, 0.0, 0.0, -2.0, -1.0}; + a_t iy= { 0.0, 0.0, 1.0, 1.0, 0.0, 0.0}; + a_t re = {3.1415926535898, 0.69314718055995, 1.1107207345396, 1.2260849569072, 0.23104906018665, 0.77778596920447}; + a_t im = {0.0, 0.0, -1.1107207345396, -0.34471136988768, 0.0, 0.19832484993429}; + if (sizeof(T) == 8) + for(int i=0; i < 6; ++i) + { + c_t x(rx[i], ix[i]); + c_t y(ry[i], iy[i]); + c_t r(re[i], im[i]); + auto rr = kyosu::ellint_rc(x, y); + TTS_RELATIVE_EQUAL(rr, r, pr); + } +}; diff --git a/test/unit/function/ellint_rd.cpp b/test/unit/function/ellint_rd.cpp new file mode 100644 index 00000000..59e1b882 --- /dev/null +++ b/test/unit/function/ellint_rd.cpp @@ -0,0 +1,39 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include +#include +// the values tested are those contained in the original Carlson article +// computation of elliptic integrals, Bille C. Carlson arXiv:math/9409227 + +TTS_CASE_TPL ( "Check ellint_rj over complex" + , kyosu::scalar_real_types + ) + (tts::type) +{ + auto pr = tts::prec(1.0e-3, 1.0e-8); + using c_t = kyosu::complex_t; + using a_t = std::array; + a_t rx= { 0.0, 2.0, 0.0, 0.0, 0.0 , -2.0}; + a_t ix= { 0.0, 0.0, 1.0, 0.0, 0.0 , -1.0}; + a_t ry= { 2.0, 3.0, 0.0, 0.0, -1.0 , 0.0}; + a_t iy= { 0.0, 0.0, -1.0, 1.0, 1.0 , -1.0}; + a_t rz= { 1.0, 4.0, 2.0, 0.0, 0.0 , -1.0}; + a_t iz= { 0.0, 0.0, 0.0, -1.0, 1.0 , 1.0}; + a_t re = {1.7972103521034, 0.16510527294261, 0.65933854154220, 1.2708196271910, -1.8577235439239, 1.8249027393704}; + a_t im = {0.0, 0.0, 0.0, 2.7811120159521, -0.96193450888839, -1.2218475784827}; + for(int i=0; i < 6; ++i) + { + c_t x(rx[i], ix[i]); + c_t y(ry[i], iy[i]); + c_t z(rz[i], iz[i]); + c_t r(re[i], im[i]); + auto rr = kyosu::ellint_rd(x, y, z); + TTS_RELATIVE_EQUAL(rr, r, pr); + } +}; diff --git a/test/unit/function/ellint_rf.cpp b/test/unit/function/ellint_rf.cpp new file mode 100644 index 00000000..b21a2c85 --- /dev/null +++ b/test/unit/function/ellint_rf.cpp @@ -0,0 +1,41 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include +#include +// the values tested are those contained in the original Carlson article +// computation of elliptic integrals, Bille C. Carlson arXiv:math/9409227 + +TTS_CASE_TPL ( "Check ellint_rf over complex" + , kyosu::scalar_real_types + ) + (tts::type) +{ + auto pr = tts::prec(4.0e-2, 1.0e-8); + using c_t = kyosu::complex_t; + using a_t = std::array; + a_t rx= {-1.0, 1.0, 0.0, 2.0, 0.0, -1.0, 0.0}; + a_t ix= { 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0}; + a_t ry= { 0.0, 2.0, 0.0, 3.0, 0.0, 0.0, 0.0}; + a_t iy= { 1.0, -0.0, -1.0, 0.0, -1.0, 1.0, 1.0}; + a_t rz= { 0.0, 0.0, 0.0, 4.0, 2.0, 1.0, 0.0}; + a_t iz= { 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 1.0}; + + a_t re = {0.79612586584234, 1.3110287771461, 1.8540746773014, 0.58408284167715, 1.0441445654064, 0.93912050218619, 1.1107207345396}; + a_t im = {-1.2138566698365, 0.0, 0.0, 0.0, 0.0, -0.5329625201863, -1.1107207345396}; + if (sizeof(T) == 8) + for(int i=0; i <= 6; ++i) + { + c_t x(rx[i], ix[i]); + c_t y(ry[i], iy[i]); + c_t z(rz[i], iz[i]); + c_t r(re[i], im[i]); + auto rr = kyosu::ellint_rf(x, y, z); + TTS_RELATIVE_EQUAL(rr, r, pr); + } +}; diff --git a/test/unit/function/ellint_rg.cpp b/test/unit/function/ellint_rg.cpp new file mode 100644 index 00000000..fd28deb8 --- /dev/null +++ b/test/unit/function/ellint_rg.cpp @@ -0,0 +1,40 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include +#include +// the values tested are those contained in the original Carlson article +// computation of elliptic integrals, Bille C. Carlson arXiv:math/9409227 + +TTS_CASE_TPL ( "Check ellint_rf over complex" + , kyosu::scalar_real_types + ) + (tts::type) +{ + auto pr = tts::prec(1.0e-3, 1.0e-8); + using c_t = kyosu::complex_t; + using a_t = std::array; + a_t rx= { 0.0, 2.0, 0.0, -1.0, 0.0, 0.0}; + a_t ix= { 0.0, 0.0, 0.0, 1.0, -1.0, 0.0}; + a_t ry= { 16.0, 3.0, 0.0, 0.0, -1.0, 0.0796}; + a_t iy= { 0.0, 0.0, 1.0, 1.0, 1.0, 0.0}; + a_t rz= { 16.0, 4.0, 0.0, 0.0, 0.0, 4.0}; + a_t iz= { 0.0, 0.0, -1.0, 0.0, 1.0, 0.0}; + + a_t re = {3.1415926535898, 1.7255030280692, 0.42360654239699, 0.44660591677018, 0.36023392184473, 1.0284758090288}; + a_t im = {0.0, 0.0, 0.0, 0.70768352357515, 0.40348623401722, 0.0 }; + for(int i=0; i <= 5; ++i) + { + c_t x(rx[i], ix[i]); + c_t y(ry[i], iy[i]); + c_t z(rz[i], iz[i]); + c_t r(re[i], im[i]); + auto rr = kyosu::ellint_rg(x, y, z); + TTS_RELATIVE_EQUAL(rr, r, pr); + } +}; diff --git a/test/unit/function/ellint_rj.cpp b/test/unit/function/ellint_rj.cpp new file mode 100644 index 00000000..258fd56b --- /dev/null +++ b/test/unit/function/ellint_rj.cpp @@ -0,0 +1,42 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include +#include +// the values tested are those contained in the original Carlson article +// computation of elliptic integrals, Bille C. Carlson arXiv:math/9409227 + +TTS_CASE_TPL ( "Check ellint_rj over complex" + , kyosu::scalar_real_types + ) + (tts::type) +{ + auto pr = tts::prec(4.0e-3, 1.0e-8); + using c_t = kyosu::complex_t; + using a_t = std::array; + a_t rx= { 0.0, 2.0, 2.0, 0.0, -1.0, 0.0, -1.0, -1.0}; + a_t ix= { 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0}; + a_t ry= { 1.0, 3.0, 3.0, 0.0, -1.0, 0.0, -1.0, -2.0}; + a_t iy= { 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0, -1.0}; + a_t rz= { 2.0, 4.0, 4.0, 0.0, 1.0, 0.0, 1.0, 0.0}; + a_t iz= { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0}; + a_t rp= { 3.0, 5.0, -1.0, 2.0, 2.0, 1.0, -3.0, -1.0}; + a_t ip= { 0.0, 0.0, 1.0, 0.0, 0.0, -1.0, 1.0, 1.0}; + a_t re = {0.77688623778582, 0.14297579667157, 0.13613945827771, 1.6490011662711, 0.94148358841220, 1.8260115229009, -0.61127970812028 , 1.8249027393704}; + a_t im = {0.0, 0.0, -0.38207561624427, 0.0, 0.0 , 1.2290661908643, -1.0684038390007, -1.2218475784827}; + for(int i=0; i < 8; ++i) + { + c_t x(rx[i], ix[i]); + c_t y(ry[i], iy[i]); + c_t z(rz[i], iz[i]); + c_t p(rp[i], ip[i]); + c_t r(re[i], im[i]); + auto rr = kyosu::ellint_rj(x, y, z, p); + TTS_RELATIVE_EQUAL(rr, r, pr); + } +}; diff --git a/test/unit/function/jacobi_elliptic.cpp b/test/unit/function/jacobi_elliptic.cpp new file mode 100644 index 00000000..53172cc1 --- /dev/null +++ b/test/unit/function/jacobi_elliptic.cpp @@ -0,0 +1,37 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include +#include + +TTS_CASE_TPL ( "Check jacobi_elliptic over complex" + , kyosu::scalar_real_types + ) + (tts::type) +{ + auto pr = tts::prec(4.0e-2, 1.0e-8); + using c_t = kyosu::complex_t; + using a_t = std::array; + a_t phi= {0.000000000000000e+00, 1.346396851538483e+00, 2.692793703076966e+00, 4.039190554615448e+00, 5.385587406153931e+00, 6.731984257692414e+00, 8.078381109230897e+00, }; + a_t m= {0.000000000000000e+00, 1.428571428571428e-01, 2.857142857142857e-01, 4.285714285714285e-01, 5.714285714285714e-01, 7.142857142857143e-01, 8.571428571428572e-01, }; + a_t rsn = {0.000000000000000e+00, 2.433149110823434e+00, 1.283941137443174e+00, -4.345523719909591e-01, -1.146031599989070e+00, -1.188244631240457e+00, -9.388455774794550e-01}; + a_t isn = {1.612004689955156e+03, 2.812489928480895e-01, -9.070417680308518e-01, -1.747485822956206e-01, 1.124636778394672e-01, 5.567227280212359e-02, 0.000000000000000e+00}; + a_t rcn = {1.612005000127921e+03, -3.080155826452643e-01, 1.120246188909950e+00, 9.211391985456892e-01, 2.183526381145390e-01, 1.021647932356633e-01, 3.443384696011595e-01}; + a_t icn = {-0.000000000000000e+00, 2.221708171357144e+00, 1.039582415796653e+00, -8.243858372166833e-02, 5.902696196755409e-01, 6.475056344849072e-01, 0.000000000000000e+00}; + a_t rdn = {1.000000000000000e+00, -4.592122381685098e-01, -9.426829656030340e-01, -9.660706784861853e-01, -5.256883155283983e-01, 2.102487586797839e-01, 4.944569734110357e-01}; + a_t idn = {-0.000000000000000e+00, 2.128865416496776e-01, -3.529707325786163e-01, 3.368760217314167e-02, -1.401014086514837e-01, 2.247415316214260e-01, 0.000000000000000e+00}; + + for(int i=0; i < 7; ++i) + { + c_t z(phi[i], phi[6-i]); + auto [sn, cn, dn] = kyosu::jacobi_elliptic(z, eve::sqrt(m[i])); + TTS_RELATIVE_EQUAL(sn, kyosu::complex(rsn[i], isn[i]), pr); + TTS_RELATIVE_EQUAL(cn, kyosu::complex(rcn[i], icn[i]), pr); + TTS_RELATIVE_EQUAL(dn, kyosu::complex(rdn[i], idn[i]), pr); + } +};