From d9730f101df4b361cd4978fd7eef47b5a237b342 Mon Sep 17 00:00:00 2001 From: jtlap Date: Wed, 15 Jan 2025 08:23:07 +0100 Subject: [PATCH] rg --- include/kyosu/functions.hpp | 1 + include/kyosu/functions/elint_rg.hpp | 104 ++++++++++++++++++++++++ include/kyosu/functions/ellint_rg.hpp | 110 ++++++++++++++++++++++++++ test/unit/function/ellint_rg.cpp | 38 +++++++++ test/unit/function/ellint_rg.hpp | 39 +++++++++ 5 files changed, 292 insertions(+) create mode 100644 include/kyosu/functions/elint_rg.hpp create mode 100644 include/kyosu/functions/ellint_rg.hpp create mode 100644 test/unit/function/ellint_rg.cpp create mode 100644 test/unit/function/ellint_rg.hpp diff --git a/include/kyosu/functions.hpp b/include/kyosu/functions.hpp index 2758a8b9..bc44840a 100644 --- a/include/kyosu/functions.hpp +++ b/include/kyosu/functions.hpp @@ -77,6 +77,7 @@ #include #include #include +#include #include #include #include diff --git a/include/kyosu/functions/elint_rg.hpp b/include/kyosu/functions/elint_rg.hpp new file mode 100644 index 00000000..0ec744bf --- /dev/null +++ b/include/kyosu/functions/elint_rg.hpp @@ -0,0 +1,104 @@ +//====================================================================================================================== +/* + 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(real auto x, real auto y, real auto z) noexcept; // 1 +//! constexpr auto ellint_rg(complex auto x, complex auto y, complex 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$\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$ is returned: +//! 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 + { + eve::swap_if(x < y, x, y); + eve::swap_if(x < z, x, z); + eve::swap_if(y > z, y, z); + // now all(x >= z) and all(z >= y) + return average(z*ellint_rf[o](x,y,z)-(x-z)*(y-z)*ellint_rd[o](x, y, z)*eve::third(as()), kyosu::sqrt(x*y/z)); + } +} diff --git a/include/kyosu/functions/ellint_rg.hpp b/include/kyosu/functions/ellint_rg.hpp new file mode 100644 index 00000000..667e8b34 --- /dev/null +++ b/include/kyosu/functions/ellint_rg.hpp @@ -0,0 +1,110 @@ +//====================================================================================================================== +/* + 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(real auto x, real auto y, real auto z) noexcept; // 1 +//! constexpr auto ellint_rg(complex auto x, complex auto y, complex 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$\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$ is returned: +//! 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/test/unit/function/ellint_rg.cpp b/test/unit/function/ellint_rg.cpp new file mode 100644 index 00000000..250c1805 --- /dev/null +++ b/test/unit/function/ellint_rg.cpp @@ -0,0 +1,38 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include +#include + +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_rg.hpp b/test/unit/function/ellint_rg.hpp new file mode 100644 index 00000000..cb0c4a11 --- /dev/null +++ b/test/unit/function/ellint_rg.hpp @@ -0,0 +1,39 @@ +//====================================================================================================================== +/* + Kyosu - Complex Without Complexes + Copyright : KYOSU Contributors & Maintainers + SPDX-License-Identifier: BSL-1.0 +*/ +//====================================================================================================================== +#include +#include +#include + +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, 2.0, -1.0, 0.0, 0.0}; + a_t ix= { 1.0, 0.0, 0.0, 0.0, 1.0, -1.0, 0.0}; + a_t ry= { 16.0, 3.0, 0.0, 3.0, 0.0, -1.0, 0.0.0796}; + a_t iy= { 0.0, -0.0, 1.0, 0.0, 1.0, 1.0, 0.0}; + a_t rz= { 16.0, 4.0, 0.0, 4.0, 0.0, 0.0, 4.0}; + a_t iz= { 0.0, 0.0, -1.0, 0.0, 0.0, 1.0, 0.0}; + + a_t re = {3.1415926535898, 1.7255030280692, 0.42360654239699, 0.44660591677018, 0.36023392184473, 1.0284 75809 0288}; + a_t im = {0.0, 0.0, 0.0, 0.70768352357515, 0.40348623401722, 0.0 }; +// 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_rg(x, y, z); + TTS_RELATIVE_EQUAL(rr, r, pr); + } +};