Jacobi elliptic am and heuman_lambda

include/eve/module/core/decorator/core.hpp

@@ -15,10 +15,6 @@
 #include <eve/as.hpp>
 #include <cfenv>
 
-//======================================================================================================================
-// New option style  - TODO rename later without the '2'
-// almost done : saturated is waiting for convert
-//======================================================================================================================
 namespace eve
 {
   struct almost_mode          {};
@@ -27,16 +23,19 @@ namespace eve
   struct condon_shortley_mode {};
   struct cylindrical_mode     {};
   struct definitely_mode      {};
+  struct eccentric_mode       {};
   struct kind_1_mode          {};
   struct kind_2_mode          {};
   struct left_mode            {};
+  struct modular_mode         {};
   struct numeric_mode         {};
   struct p_kind_mode          {};
   struct promote_mode         {};
   struct q_kind_mode          {};
   struct right_mode           {};
   struct spherical_mode       {};
   struct successor_mode       {};
+  struct threshold_mode       {};
 
   struct upper_mode       {static constexpr int value = FE_UPWARD;     };
   struct lower_mode       {static constexpr int value = FE_DOWNWARD;   };
@@ -58,9 +57,11 @@ namespace eve
   [[maybe_unused]] inline constexpr auto condon_shortley  = ::rbr::flag( condon_shortley_mode{} );
   [[maybe_unused]] inline constexpr auto cylindrical      = ::rbr::flag( cylindrical_mode{}     );
   [[maybe_unused]] inline constexpr auto downward         = ::rbr::flag( downward_mode{}        );
+  [[maybe_unused]] inline constexpr auto eccentric        = ::rbr::flag( eccentric_mode{}            );
   [[maybe_unused]] inline constexpr auto kind_1           = ::rbr::flag( kind_1_mode{}          );
   [[maybe_unused]] inline constexpr auto kind_2           = ::rbr::flag( kind_2_mode{}          );
   [[maybe_unused]] inline constexpr auto left             = ::rbr::flag( left_mode{}            );
+  [[maybe_unused]] inline constexpr auto modular          = ::rbr::flag( modular_mode{}         );
   [[maybe_unused]] inline constexpr auto numeric          = ::rbr::flag( numeric_mode{}         );
   [[maybe_unused]] inline constexpr auto pedantic         = ::rbr::flag( pedantic_mode{}        );
   [[maybe_unused]] inline constexpr auto p_kind           = ::rbr::flag( p_kind_mode{}          );
@@ -83,9 +84,11 @@ namespace eve
   struct compensated_option     : detail::exact_option<compensated>     {};
   struct condon_shortley_option : detail::exact_option<condon_shortley> {};
   struct cylindrical_option     : detail::exact_option<cylindrical>     {};
+  struct eccentric_option       : detail::exact_option<eccentric>       {};
   struct kind_1_option          : detail::exact_option<kind_1>          {};
   struct kind_2_option          : detail::exact_option<kind_2>          {};
   struct left_option            : detail::exact_option<left>            {};
+  struct modular_option         : detail::exact_option<modular>         {};
   struct numeric_option         : detail::exact_option<numeric>         {};
   struct p_kind_option          : detail::exact_option<p_kind>          {};
   struct promote_option         : detail::exact_option<promote>         {};
@@ -153,4 +156,20 @@ namespace eve
 
     EVE_FORCEINLINE constexpr auto default_to(auto const& base) const { return base; }
   };
+
+  // New threshold option that carry a value
+  template<typename Value> struct threshold_t;
+
+  struct threshold_option
+  {
+    template<typename Value>
+    EVE_FORCEINLINE constexpr auto process(auto const& base, threshold_t<Value> const& opts) const
+    {
+      auto news = rbr::merge(options{opts}, base);
+      return options<decltype(news)>{news};
+    }
+
+    EVE_FORCEINLINE constexpr auto default_to(auto const& base) const { return base; }
+  };
+
 }

include/eve/module/core/detail/tolerance.hpp

@@ -115,4 +115,53 @@ namespace eve
   };
 
   inline constexpr definitely_t<default_tolerance> definitely = {};
+
+
+  //  ============================================================================================
+  // threshold
+
+  template<typename Value>
+  struct threshold_t
+  {
+    constexpr threshold_t() {}
+    constexpr explicit threshold_t(Value v) : value_(v) {}
+    constexpr threshold_t(threshold_t const& v) : value_(v.value_) {}
+    constexpr threshold_t& operator=(threshold_t const& v) { value_ = v.value_; return *this; }
+
+    /// ID type associated to the keyword
+    using id_type = threshold_mode;
+
+    template<typename T> static constexpr bool accept() { return true; }
+
+    std::ostream& show(std::ostream& os, auto const& v) const
+    {
+      return os << "threshold by " << v.value_;
+    }
+
+    using tag_type          = threshold_t<default_tolerance>;
+    using keyword_type      = threshold_t<default_tolerance>;
+    using stored_value_type = threshold_t<Value>;
+
+    template<eve::scalar_value Type>
+    constexpr auto operator=(Type v) const noexcept { return threshold_t<Type>{v}; }
+
+    template<typename T>
+    constexpr auto operator=(threshold_t<T> const& v) const noexcept
+    {
+      return v;
+    }
+
+    constexpr auto operator()(keyword_type const&) const noexcept { return *this; }
+
+    template<typename T> constexpr auto value(T const&) const
+    {
+      using type = element_type_t<T>;
+      if constexpr(std::same_as<Value,default_tolerance>) return eps(as<type>{});
+      else                                                return type{value_};
+    }
+
+    Value value_;
+  };
+
+  inline constexpr threshold_t<default_tolerance> threshold = {};
 }

include/eve/module/elliptic/am.hpp

@@ -0,0 +1,157 @@
+//==================================================================================================
+/**
+  EVE - Expressive Vector Engine
+  Copyright : EVE Contributors & Maintainers
+  SPDX-License-Identifier: MIT
+**/
+//==================================================================================================
+#pragma once
+
+#include <eve/arch.hpp>
+#include <eve/traits/overload.hpp>
+#include <eve/module/core/decorator/core.hpp>
+#include <eve/module/core.hpp>
+#include <eve/module/math.hpp>
+#include <eve/traits/common_value.hpp>
+
+namespace eve
+{
+  template<typename Options>
+  struct am_t : elementwise_callable<am_t, Options, threshold_option, modular_option, eccentric_option>
+  {
+    template<eve::floating_value T>
+    constexpr EVE_FORCEINLINE
+    T operator()(T a) const noexcept { return EVE_DISPATCH_CALL(a); }
+
+    template<eve::floating_value T0, eve::floating_value T1>
+    requires (same_lanes_or_scalar<T0, T1>)
+    constexpr EVE_FORCEINLINE
+    eve::common_value_t<T0, T1> operator()(T0 a, T1 b) const noexcept
+    { return EVE_DISPATCH_CALL(a, b); }
+
+    EVE_CALLABLE_OBJECT(am_t, am_);
+  };
+
+//================================================================================================
+//! @addtogroup elliptic
+//! @{
+//!   @var am
+//!   @brief `elementwise_callable` object computing the  Jacobi's Amplitude function.
+//!
+//!   @groupheader{Header file}
+//!
+//!   @code
+//!   #include <eve/module/elliptic.hpp>
+//!   @endcode
+//!
+//!   @groupheader{Callable Signatures}
+//!
+//!   @code
+//!   namespace eve
+//!   {
+//!      // Regular overload
+//!      constexpr auto am(floating_value auto u, floating_value auto x)                 noexcept; // 1
+//!
+//!      //Semantic modifiers
+//!      constexpr auto am[modular](floating_value auto u, floating_value auto alpha)    noexcept; // 1
+//!      constexpr auto am[eccentric(floating_value auto u, floating_value auto k)       noexcept; // 1
+//!      constexpr auto am[threshold = tol](floating_value auto u, floating_value auto x) noexcept; // 1
+//!
+//!      // Lanes masking
+//!      constexpr auto am[conditional_expr auto c](/*any of the above overloads*/)       noexcept; // 2
+//!      constexpr auto am[logical_value autolm](/*any of the above overloads*/)          noexcept; // 2
+//!   }
+//!   @endcode
+//!
+//!   **Parameters**
+//!
+//!     * `u`: argument.
+//!     * `x`: amplitude parameter (\f$0\le m\le 1).
+//!     * `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).
+//!      2. [The operation is performed conditionally](@ref conditional)
+//!
+//! @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 (give the right option).
+
+//!  @groupheader{External references}
+//!   *  [C++ standard reference: am](https://en.cppreference.com/w/cpp/numeric/special_functions/am)
+//!   *  [DLMF: Jacobi zeta](https://dlmf.nist.gov/22.16)
+//!   *  [Wolfram MathWorld: Jacobi Zeta Function](https://mathworld.wolfram.com/JacobiZetaFunction.html)
+//!   *  [Boost: jacobi_zeta](https://beta.boost.org/doc/libs/1_84_0/libs/math/doc/html/math_toolkit/ellint/jacobi_zeta.html)
+//!
+//!  @groupheader{Example}
+//!  @godbolt{doc/elliptic/am.cpp}
+//================================================================================================
+  inline constexpr auto am = functor<am_t>;
+//================================================================================================
+//! @}
+//================================================================================================
+
+  namespace detail
+  {
+    template<floating_value T, callable_options O>
+    T am_(EVE_REQUIRES(cpu_), O const& o, T u, T x) noexcept
+    {
+      auto tol = [&](){
+        if constexpr (O::contains(threshold)) return o[threshold].value(x);
+        else  return eve::epsilon(x);
+      }();
+      x =  abs(x);
+      if (O::contains(modular)) x = sin(x);
+      else if (O::contains(eccentric)) x = sqrt(x);
+
+      auto xx = eve::abs(x);
+      auto xxisone = xx == one(as(x));
+      auto k = if_else(xxisone, T(0.5), xx);
+      auto am_kernel = [u, k, tol]()
+        {
+          constexpr int N = 30; //10 may be sufficent
+          std::array<T, N+1> a, g, c;
+          a[0] = one(as(u));
+          g[0] = sqrt(oneminus(sqr(k)));
+          c[0] = k;
+
+          // Perform the sequence of Gaussian transformations of arithmetic and
+          // geometric means of successive arithmetic and geometric means until
+          // the two means converge to a common mean (upto asked threshold)
+          // starting with a = 1 and g = k', which were set above.
+
+          auto two_n =  one(as(u));
+          int n = 0;
+          for (; n < N; n++)
+          {
+            if (eve::all(eve::is_not_greater_equal(abs(c[n]), tol))) break;
+            two_n += two_n;
+            a[n+1] = average(a[n], g[n]);
+            g[n+1] = sqrt(a[n] * g[n]);
+            c[n+1] = average(a[n], -g[n]);
+          }
+
+          // Prepare for the inverse transformation of phi = x * cm.
+          auto phi = two_n * a[n] * u;
+
+          // Perform backward substitution
+          for (; n > 0; --n)
+          {
+            phi = average(phi, asin(c[n]*sin(phi)/a[n]) );
+          }
+          return if_else(is_infinite(u) && is_eqz(k), u, phi);
+        };
+
+      T r;
+      if (eve::any(xxisone)) r = fms(T(2), atan(exp(u)), pio_2(as(xx)));
+      if (eve::all(xxisone)) return r;
+      return if_else(xxisone, r, am_kernel());
+    }
+  }
+}

include/eve/module/elliptic/elliptic.hpp

@@ -15,3 +15,7 @@
 #include <eve/module/elliptic/ellint_rf.hpp>
 #include <eve/module/elliptic/ellint_rg.hpp>
 #include <eve/module/elliptic/ellint_rj.hpp>
+#include <eve/module/elliptic/am.hpp>
+#include <eve/module/elliptic/jacobi_elliptic.hpp>
+#include <eve/module/elliptic/jacobi_zeta.hpp>
+#include <eve/module/elliptic/heuman_lambda.hpp>