carlson elliptic

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"

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;
   }

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<Z>;
     auto twoopi = eve::two_o_pi(eve::as<u_t>());
     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<u_t>())), 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<Z>;
     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<Z, 0>();
-    return cb_yn(n, z, dummy);
+    return  _::cb_yn(n, z, dummy);
   }
 
   //===-------------------------------------------------------------------------------------------

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

include/kyosu/functions.hpp

@@ -38,6 +38,7 @@
 #include <kyosu/functions/acoth.hpp>
 #include <kyosu/functions/acotpi.hpp>
 #include <kyosu/functions/agd.hpp>
+#include <kyosu/functions/am.hpp>
 #include <kyosu/functions/atanpi.hpp>
 #include <kyosu/functions/acsc.hpp>
 #include <kyosu/functions/acscpi.hpp>
@@ -72,6 +73,12 @@
 #include <kyosu/functions/digamma.hpp>
 #include <kyosu/functions/dist.hpp>
 #include <kyosu/functions/dot.hpp>
+#include <kyosu/functions/ellint_fe.hpp>
+#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>
 #include <kyosu/functions/erfi.hpp>
@@ -115,6 +122,7 @@
 #include <kyosu/functions/is_pure.hpp>
 #include <kyosu/functions/is_real.hpp>
 #include <kyosu/functions/is_unitary.hpp>
+#include <kyosu/functions/jacobi_elliptic.hpp>
 #include <kyosu/functions/lambda.hpp>
 #include <kyosu/functions/lbeta.hpp>
 #include <kyosu/functions/ldiv.hpp>

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 <kyosu/details/callable.hpp>
+#include <kyosu/functions/jacobi_elliptic.hpp>
+#include <eve/module/elliptic/ellint_1.hpp>
+#include <eve/module/elliptic/am.hpp>
+#include <kyosu/functions/sqr.hpp>
+
+namespace kyosu
+{
+
+  template<typename Options>
+  struct am_t : eve::strict_elementwise_callable<am_t, Options, eve::modular_option,
+                                                 eve::eccentric_option, eve::threshold_option>
+  {
+    template<concepts::real T0, concepts::real T1>
+    constexpr KYOSU_FORCEINLINE
+    auto operator()(T0 a, T1 b) const noexcept
+    { return (*this)(kyosu::complex(a), b); }
+
+    template<concepts::complex T0, concepts::real T1>
+    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 <kyosu/functions.hpp>
+//!   @endcode
+//!
+//!   @groupheader{Callable Signatures}
+//!
+//!   @code
+//!   namespace kyosu
+//!   {
+//!      // Regular overload
+//!       constexpr Z       am( real<Z>   z, real<U> m)               noexcept; //1
+//!       constexpr auto    am(complex<Z> z, real<U> m)               noexcept; //1
+//!
+//!      //Semantic modifiers
+//!      constexpr Z    am[modular]( real<Z>      z, real<U> alpha)   noexcept; //1
+//!      constexpr auto am[modular](complex<Z>    z, real<U> alpha)   noexcept; //1
+//!      constexpr Z    am[eccentric]( real<Z>    z, real<U> m)       noexcept; //1
+//!      constexpr auto am[eccentric](complex<Z>  z, real<U> m)       noexcept; //1
+//!      constexpr auto am[threshold = tol]( real<Z>   z, real<U> k)  noexcept; //1
+//!      constexpr auto am[threshold = tol](complex<Z> z, real<U> 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<am_t>;
+//======================================================================================================================
+//! @}
+//======================================================================================================================
+}
+
+namespace kyosu::_
+{
+  template<typename Z, typename M, eve::callable_options O>
+  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;
+  }
+}

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**
 //!