Skip to content

Commit

Permalink
rg
Browse files Browse the repository at this point in the history
  • Loading branch information
@jtlap
jtlap committed Jan 15, 2025
1 parent edb764c commit d9730f1
Show file tree
Hide file tree
Showing 5 changed files with 292 additions and 0 deletions.
1 change: 1 addition & 0 deletions include/kyosu/functions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@
#include <kyosu/functions/ellint_rc.hpp>
#include <kyosu/functions/ellint_rd.hpp>
#include <kyosu/functions/ellint_rf.hpp>
#include <kyosu/functions/ellint_rg.hpp>
#include <kyosu/functions/ellint_rj.hpp>
#include <kyosu/functions/erf.hpp>
#include <kyosu/functions/erfcx.hpp>
Expand Down
104 changes: 104 additions & 0 deletions include/kyosu/functions/elint_rg.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,104 @@
//======================================================================================================================
/*
Kyosu - Complex Without Complexes
Copyright: KYOSU Contributors & Maintainers
SPDX-License-Identifier: BSL-1.0
*/
//======================================================================================================================
#pragma once
#include <eve/module/elliptic/ellint_rj.hpp>
#include <kyosu/details/callable.hpp>
#include <kyosu/functions/ellint_rf.hpp>
#include <kyosu/functions/ellint_rd.hpp>

namespace kyosu
{

template<typename Options>
struct ellint_rg_t : eve::elementwise_callable<ellint_rg_t, Options, eve::threshold_option>
{
template<eve::floating_value T0, eve::floating_value T1,
eve::floating_value T2>
constexpr KYOSU_FORCEINLINE
auto operator()(T0 a, T1 b, T2 c) const noexcept
{ return eve::ellint_rg(a, b, c); }

template<concepts::complex T0, concepts::complex T1, concepts::complex T2, >
constexpr KYOSU_FORCEINLINE
auto operator()(T0 a, T1 b, T2 c) const noexcept
requires(concepts::complex<T0> || concepts::complex<T1> || concepts::complex<T2>)
{ 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 <kyosu/kyosu.hpp>
//! @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<ellint_rg_t>;
//================================================================================================
//! @}
}
namespace kyosu::_
{
template<typename T, eve::callable_options O >
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<r_t>()), kyosu::sqrt(x*y/z));
}
}
110 changes: 110 additions & 0 deletions include/kyosu/functions/ellint_rg.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,110 @@
//======================================================================================================================
/*
Kyosu - Complex Without Complexes
Copyright: KYOSU Contributors & Maintainers
SPDX-License-Identifier: BSL-1.0
*/
//======================================================================================================================
#pragma once
#include <eve/module/elliptic/ellint_rj.hpp>
#include <kyosu/details/callable.hpp>
#include <kyosu/functions/ellint_rf.hpp>
#include <kyosu/functions/ellint_rd.hpp>

namespace kyosu
{

template<typename Options>
struct ellint_rg_t : eve::elementwise_callable<ellint_rg_t, Options, eve::threshold_option>
{
template<eve::floating_value T0, eve::floating_value T1,
eve::floating_value T2>
constexpr KYOSU_FORCEINLINE
auto operator()(T0 a, T1 b, T2 c) const noexcept
{ return eve::ellint_rg(a, b, c); }

template<concepts::complex T0, concepts::complex T1, concepts::complex T2>
constexpr KYOSU_FORCEINLINE
auto operator()(T0 a, T1 b, T2 c) const noexcept
requires(concepts::complex<T0> || concepts::complex<T1> || concepts::complex<T2>)
{ 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 <kyosu/kyosu.hpp>
//! @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<ellint_rg_t>;
//================================================================================================
//! @}
}
namespace kyosu::_
{
template<typename T, eve::callable_options O >
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<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<r_t>()), root);
}
}
38 changes: 38 additions & 0 deletions test/unit/function/ellint_rg.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,38 @@
//======================================================================================================================
/*
Kyosu - Complex Without Complexes
Copyright : KYOSU Contributors & Maintainers
SPDX-License-Identifier: BSL-1.0
*/
//======================================================================================================================
#include <iomanip>
#include <kyosu/kyosu.hpp>
#include <test.hpp>

TTS_CASE_TPL ( "Check ellint_rf over complex"
, kyosu::scalar_real_types
)
<typename T>(tts::type<T>)
{
auto pr = tts::prec<T>(1.0e-3, 1.0e-8);
using c_t = kyosu::complex_t<T>;
using a_t = std::array<T, 6 >;
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);
}
};
39 changes: 39 additions & 0 deletions test/unit/function/ellint_rg.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
//======================================================================================================================
/*
Kyosu - Complex Without Complexes
Copyright : KYOSU Contributors & Maintainers
SPDX-License-Identifier: BSL-1.0
*/
//======================================================================================================================
#include <iomanip>
#include <kyosu/kyosu.hpp>
#include <test.hpp>

TTS_CASE_TPL ( "Check ellint_rf over complex"
, kyosu::scalar_real_types
)
<typename T>(tts::type<T>)
{
auto pr = tts::prec<T>(1.0e-3, 1.0e-8);
using c_t = kyosu::complex_t<T>;
using a_t = std::array<T, 7 >;
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);
}
};

0 comments on commit d9730f1

Please sign in to comment.