Contfrac

include/eve/module/contfrac.hpp

@@ -0,0 +1,23 @@
+//==================================================================================================
+/*
+  EVE - Expressive Vector Engine
+  Copyright : EVE Project Contributors
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#pragma once
+
+//==================================================================================================
+//! @addtogroup functions
+//! @{
+//! @defgroup contfrac Continuous fractions
+//!
+//! @brief  This module provides an simd compatible implementation of the lentz algorithm.
+//!         to evaluate continuous fractions.
+//!
+//! **Convenience header:** @code{.cpp} #include <eve/module/contfrac.hpp> @endcode
+//!
+//! @}
+//==================================================================================================
+#include <eve/module/contfrac/lentz_a.hpp>
+#include <eve/module/contfrac/lentz_b.hpp>

include/eve/module/contfrac/impl/lentz.hpp

@@ -0,0 +1,135 @@
+//==================================================================================================
+/*
+  EVE - Expressive Vector Engine
+  Copyright : EVE Project Contributors
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#pragma once
+
+namespace eve::detail
+{
+
+  template <typename T>
+  struct is_pair : public std::false_type{};
+
+  template <typename T, typename U>
+  struct is_pair<kumi::tuple<T,U>> : public std::true_type{};
+
+  //
+  // continued_fraction_b
+  // Evaluates:
+  //
+  //               a1
+  //     b0 + -----------
+  //                   a2
+  //          b1 + -----------
+  //                        a3
+  //               b2 + ----------
+  //                       ...
+  //
+  // Note that the first a0 returned by generator Gen is discarded.
+  //
+
+  template <typename Gen, typename U>
+  EVE_FORCEINLINE constexpr auto lentz_b_(EVE_SUPPORTS(cpu_), Gen g, const U& eps, size_t max_terms)
+    noexcept
+  {
+    using eve::abs;
+    auto v = g();
+    constexpr auto pure_pair = is_pair<decltype(v)>::value;
+    using tmp_t = decltype([v](){ if constexpr(pure_pair) return get<0>(v); else  return v;}());
+    using r_t = std::decay_t<tmp_t>;
+    using u_t =   underlying_type_t<r_t>;
+    u_t tiny = 16*smallestposval(as<u_t>()) ;
+    u_t terminator(eps <= 0 ? eve::abs(eps)*eve::eps(as<u_t>()) : u_t(eps));
+    size_t counter(max_terms);
+
+    r_t f{};
+    if constexpr(pure_pair) f = get<1>(v); else f = v;
+    f = if_else(is_eqz(f), tiny, f);
+    r_t C = f;
+    r_t D{0};
+    r_t delta{};
+
+    do{
+      v = g();
+      if constexpr(pure_pair)
+      {
+        D = fam(get<1>(v), get<0>(v), D);
+        C = fam(get<1>(v), get<0>(v), rec(C));
+      }
+      else
+      {
+        D = v+D;
+        C = v+rec(C);
+      }
+      D = if_else(is_eqz(D), tiny, D);
+      C = if_else(is_eqz(C), tiny, C);
+      D = rec(D);
+      delta = C*D;
+      f *= delta;
+    } while (any(abs(dec(delta)) > terminator) && --counter);
+    return f;
+  }
+
+  //
+  // continued_fraction_a
+  // Evaluates:
+  //
+  //            a0
+  //      ------------
+  //                a1
+  //      b0 + -----------
+  //                    a2
+  //           b1 + ----------
+  //                        a3
+  //                 b2 + -------
+  //                        ...
+  //
+  // Note that the first a0 and b0 returned by generator Gen are both used.
+  //
+
+  template <typename Gen, typename U>
+  EVE_FORCEINLINE auto lentz_a_(EVE_SUPPORTS(cpu_), Gen g, const U& eps, size_t max_terms)  noexcept
+  {
+    using eve::abs;
+    auto v = g();
+    constexpr auto pure_pair = is_pair<decltype(v)>::value;
+    using tmp_t = decltype([v](){ if constexpr(pure_pair) return get<0>(v); else  return v;}());
+    using r_t = std::decay_t<tmp_t>;
+    using u_t =   underlying_type_t<r_t>;
+    u_t tiny = 16*smallestposval(as<u_t>()) ;
+    u_t terminator(eps <= 0 ? eve::abs(eps)*eve::eps(as<u_t>()) : u_t(eps));
+    size_t counter(max_terms);
+
+    r_t f{};
+    if constexpr(pure_pair) f = get<1>(v); else f = v;
+    f = if_else(is_eqz(f), tiny, f);
+    r_t a0{};
+    if constexpr(pure_pair) a0 = get<0>(v); else a0 = one(as<r_t>());
+    a0 = if_else(is_eqz(a0), tiny, a0);
+    auto C = f;
+    r_t  D{0};
+    auto delta(D);
+    do{
+      v = g();
+      if constexpr(pure_pair)
+      {
+        D = fam(get<1>(v), get<0>(v), D);
+        C = fam(get<1>(v), get<0>(v), rec(C));
+      }
+      else
+      {
+        D = v+D;
+        C = v+rec(C);
+      }
+      D = if_else(is_eqz(D), tiny, D);
+      C = if_else(is_eqz(C), tiny, C);
+      D = rec(D);
+      delta = C*D;
+      f *= delta;
+    } while (eve::any(abs(dec(delta)) > terminator) && --counter);
+    return a0/f;
+  }
+}

include/eve/module/contfrac/lentz_a.hpp

@@ -0,0 +1,67 @@
+//==================================================================================================
+/*
+  EVE - Expressive Vector Engine
+  Copyright : EVE Project Contributors
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#pragma once
+
+#include <eve/detail/overload.hpp>
+
+namespace eve
+{
+//================================================================================================
+//! @addtogroup contfrac
+//! @{
+//!   @var lentz_a
+//!   @brief Implement the lentz scheme to evaluate continued fractions
+//!
+//!   **Defined in header**
+//!
+//!   @code
+//!   #include <eve/module/contfrac.hpp>
+//!   @endcode
+//!
+//!   @groupheader{Callable Signatures}
+//!
+//!   @code
+//!   namespace eve
+//!   {
+//!     template< typename Gen, eve::scalar_ordered_value T> auto lentz_a(Gen g, const T& tol, size_t & max_terms) noexcept;
+//!   }
+//!   @endcode
+//!
+//!   **Parameters**
+//!
+//!     * `g`   : generator function.
+//!     * `tol` : tolerance value. If negative the effective tolerance will be abs(tol)*eve::eps(as(< u_t>)
+//!               where u_t is the underlying floating type  associated to the return type of the invocable g.
+//!     * `max_terms` : no more than max_terms calls to the generator will be made,
+//!
+//!   The generator type should be an invocable which supports the following operations:
+//!   *  The call to g()  returns a floating value or a pair (kumi::tuple) of such.
+//!      Each time this operator is called then the next pair of a and b values has to be returned, 
+//!      or, if result_type is not a pair type, then the next b value
+//!      has to be returned and all the a values are assumed to be equal to one.
+//!
+//!   *  In all the continued fraction evaluation functions the effective tol parameter
+//!      is the relative precision desired in the result,
+//!      The evaluation of the fraction will continue until the last term evaluated leaves
+//!      the relative error in the result less than tolerance or the max_terms iteration is reached.
+//!
+//!   **Return value**
+//!
+//!     The value of the continued fraction is returned.
+//!     \f$\displaystyle \frac{a_0}{b_0+\frac{a_1}{b_1+\frac{a_2}{b_2+\cdots\vphantom{\frac{1}{1}} }}}\f$
+//!
+//!   @groupheader{Example}
+//!
+//!   @godbolt{doc/polynomial/regular/lentz_a.cpp}
+//!
+//! @}
+//================================================================================================
+EVE_MAKE_CALLABLE(lentz_a_, lentz_a);
+}
+
+#include <eve/module/contfrac/impl/lentz.hpp>

include/eve/module/contfrac/lentz_b.hpp

@@ -0,0 +1,72 @@
+//==================================================================================================
+/*
+  EVE - Expressive Vector Engine
+  Copyright : EVE Project Contributors
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#pragma once
+
+#include <eve/detail/overload.hpp>
+
+namespace eve
+{
+//================================================================================================
+//! @addtogroup contfrac
+//! @{
+//!   @var lentz_b
+//!   @brief Implement the lentz scheme to evaluate continued fractions
+//!
+//!   **Defined in header**
+//!
+//!   @code
+//!   #include <eve/module/contfrac.hpp>
+//!   @endcode
+//!
+//!   @groupheader{Callable Signatures}
+//!
+//!   @code
+//!   namespace eve
+//!   {
+//!     template< typename Gen, eve::floating_value T> auto lentz_b(Gen g, const T& tol, size_t & max_terms) noexcept;
+//!   }
+//!   @endcode
+//!
+//!   **Parameters**
+//!
+//!     * `g`   : generator function.
+//!     * `tol` : tolerance value. If negative the effective tolerance will be abs(tol)*eve::eps(as(< u_t>)
+//!               where u_t is the underlying floating type associated to the return type of the invocable g.
+//!     * `max_terms` : no more than max_terms calls to the generator will be made.
+//!
+//!   The generator type should be an invocable  which supports the following operations:
+//!   *  The call to g()  returns a floating  value or a pair (kumi::tuple) of such.
+//!      Each time this operator is called then the next pair of a and b values  has to be returned,
+//!      or, if result_type is not a pair type, then the next b value
+//!       has to be returned and all the a values are assumed to be equal to one.
+//!
+//!   * In all the continued fraction evaluation functions the  effective tol parameter
+//!     is the relative precision desired in the result,
+//!     The evaluation of the fraction will continue until the last term evaluated leaves
+//!     the relative error in the result less than tolerance or the max_terms iteration is reached.
+//!
+//!   **Return value**
+//!
+//!     The value of the continued fraction is returned.
+//!     \f$\displaystyle b_0+\frac{a_1}{b_1+\frac{a_2}{b_2+\frac{a_3}{b_3+\cdots\vphantom{\frac{1}{1}} }}}\f$
+//!
+//!     Note that the the first a value (a0) generated is not used here.
+//!
+//!   @note the implementation is largely inspired by the boost/math/fraction one, with less requirements on the invocable.
+//!         Peculiarly lambda functions can be used.
+//!
+//!   @groupheader{Example}
+//!
+//!   @godbolt{doc/polynomial/regular/lentz_b.cpp}
+//!
+//! @}
+//================================================================================================
+EVE_MAKE_CALLABLE(lentz_b_, lentz_b);
+}
+
+#include <eve/module/contfrac/impl/lentz.hpp>

test/doc/CMakeLists.txt

@@ -101,3 +101,10 @@ glob_unit("doc" ${doc_root} "special/regular/*.cpp"  )
 add_custom_target(doc.traits.exe        )
 add_dependencies(doc.exe doc.traits.exe )
 glob_unit("doc" ${doc_root} "traits/*.cpp")
+
+##==================================================================================================
+## GLOB and process traits doc tests
+##==================================================================================================
+add_custom_target(doc.contfrac.exe        )
+add_dependencies(doc.exe doc.contfrac.exe )
+glob_unit("doc" ${doc_root} "contfrac/*.cpp")

test/doc/contfrac/lentz_a.cpp

@@ -0,0 +1,31 @@
+#include <eve/module/contfrac.hpp>
+#include <eve/wide.hpp>
+#include <iostream>
+#include <list>
+#include <vector>
+
+using w_t = eve::wide<double, eve::fixed<8>>;
+
+template <class T>
+struct const_fraction
+{
+   typedef T result_type;
+
+   result_type operator()()
+   {
+     return T{1, 2, 3, 4, 5, 6, 7, 8};
+   }
+};
+
+
+int main()
+{
+  const_fraction<w_t> func;
+  auto gr = eve::lentz_a(func,eve::eps(eve::as<double>()), 100);
+  std::cout << "The golden ratio is: " << gr << std::endl;
+  w_t z{1, 2, 3, 4, 5, 6, 7, 8};
+  auto ref = (-z+eve::sqrt(eve::sqr(z)+4))/2;
+  std::cout << "ref golden ratio is: " << ref << std::endl;
+
+  return 0;
+}

test/doc/contfrac/lentz_b.cpp

@@ -0,0 +1,45 @@
+#include <eve/module/contfrac.hpp>
+#include <eve/module/math.hpp>
+#include <eve/wide.hpp>
+#include <iostream>
+
+
+template <class T>
+struct const_fraction
+{
+  auto operator()(){ return T{1, 2, 3, 4}; }
+};
+
+template <typename T>
+struct tan_fraction
+{
+  T a, b;
+  tan_fraction(T v) : a(-v * v), b(-1) {}
+  auto operator()()
+  {
+    b += T(2);
+    return kumi::tuple{a, b};
+  }
+};
+
+template <class T>
+T mytan(T a)
+{
+  tan_fraction<T> fract(a);
+  return a/eve::lentz_b(fract, eve::eps(eve::as<eve::underlying_type_t<T>>()), 100);
+}
+
+int main()
+{
+  using w_t = eve::wide<double, eve::fixed<4>>;
+  const_fraction<w_t> func;
+  w_t zz{1, 2, 3, 4};
+  std::cout << "ref constant fracs " << (zz+eve::sqrt(eve::sqr(zz)+4))/2 << std::endl;
+  auto gr = eve::lentz_b(func,eve::eps(eve::as<double>()), 100);
+  std::cout << "    constant fracs " << gr << std::endl;
+  w_t z{eve::pio_4(eve::as<double>()), eve::pio_3(eve::as<double>()), 0.0, 1.0};
+  std::cout << "frac tan(" << z << ") is: " << mytan(z) << std::endl;
+  std::cout << "ref  tan(" << z << ") is: " << eve::tan(z) << std::endl;
+
+  return 0;
+}

test/unit/CMakeLists.txt

@@ -112,4 +112,8 @@ add_custom_target(unit.special.exe         )
 add_dependencies(unit.exe unit.special.exe )
 glob_unit("unit" ${unit_root} "module/special/*.cpp" )
 
+add_custom_target(unit.contfrac.exe         )
+add_dependencies(unit.exe unit.contfrac.exe )
+glob_unit("unit" ${unit_root} "module/contfrac/*.cpp" )
+
 endif()