diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index b5cc0c9131..82de474a32 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -20,7 +20,7 @@ jobs: documentation: runs-on: ubuntu-latest container: - image: ghcr.io/jfalcou/compilers:v6 + image: ghcr.io/jfalcou/compilers:v7 strategy: fail-fast: false matrix: diff --git a/.github/workflows/integration.yml b/.github/workflows/integration.yml index ada807bb42..17ba85d6b0 100644 --- a/.github/workflows/integration.yml +++ b/.github/workflows/integration.yml @@ -21,7 +21,7 @@ jobs: install: runs-on: [ubuntu-latest] container: - image: ghcr.io/jfalcou/compilers:v6 + image: ghcr.io/jfalcou/compilers:v7 strategy: fail-fast: false @@ -103,4 +103,3 @@ jobs: cmake ../examples/multi-arch -G Ninja ninja multi-arch ./multi-arch - diff --git a/.github/workflows/random.yml b/.github/workflows/random.yml index c016ae4030..14df38f63f 100644 --- a/.github/workflows/random.yml +++ b/.github/workflows/random.yml @@ -20,7 +20,7 @@ jobs: random: runs-on: ubuntu-latest container: - image: ghcr.io/jfalcou/compilers:v6 + image: ghcr.io/jfalcou/compilers:v7 strategy: fail-fast: false matrix: diff --git a/.github/workflows/update.yml b/.github/workflows/update.yml index 0261ea72df..c0734bf4e6 100644 --- a/.github/workflows/update.yml +++ b/.github/workflows/update.yml @@ -13,7 +13,7 @@ jobs: generate-doc: runs-on: ubuntu-latest container: - image: ghcr.io/jfalcou/compilers:v6 + image: ghcr.io/jfalcou/compilers:v7 strategy: fail-fast: false steps: diff --git a/include/eve/module/math/regular/impl/lentz.hpp b/include/eve/module/math/regular/impl/lentz.hpp new file mode 100644 index 0000000000..0ee1c1be0a --- /dev/null +++ b/include/eve/module/math/regular/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 + struct is_pair : std::false_type{}; + + template + struct is_pair> : std::true_type{}; + + // + // lentz_b + // Evaluates: + // + // a1 + // b0 + ----------- + // a2 + // b1 + ----------- + // a3 + // b2 + ---------- + // ... + // + // Note that the first a0 returned by generator Gen is discarded. + // + + template + 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::value; + using tmp_t = decltype([v](){ if constexpr(pure_pair) return get<0>(v); else return v;}()); + using r_t = std::decay_t; + using u_t = underlying_type_t; + u_t tiny = 16*smallestposval(as()) ; + u_t terminator(eps <= 0 ? eve::abs(eps)*eve::eps(as()) : 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; + } + + // + // lentz_a + // Evaluates: + // + // a0 + // ------------ + // a1 + // b0 + ----------- + // a2 + // b1 + ---------- + // a3 + // b2 + ------- + // ... + // + // Note that the first a0 and b0 returned by generator Gen are both used. + // + + template + 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::value; + using tmp_t = decltype([v](){ if constexpr(pure_pair) return get<0>(v); else return v;}()); + using r_t = std::decay_t; + using u_t = underlying_type_t; + u_t tiny = 16*smallestposval(as()) ; + u_t terminator(eps <= 0 ? eve::abs(eps)*eve::eps(as()) : 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()); + 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; + } +} diff --git a/include/eve/module/math/regular/lentz_a.hpp b/include/eve/module/math/regular/lentz_a.hpp new file mode 100644 index 0000000000..0aeb9fbf17 --- /dev/null +++ b/include/eve/module/math/regular/lentz_a.hpp @@ -0,0 +1,67 @@ +//================================================================================================== +/* + EVE - Expressive Vector Engine + Copyright : EVE Project Contributors + SPDX-License-Identifier: BSL-1.0 +*/ +//================================================================================================== +#pragma once + +#include + +namespace eve +{ +//================================================================================================ +//! @addtogroup contfrac +//! @{ +//! @var lentz_a +//! @brief Implement the lentz scheme to evaluate continued fractions +//! +//! **Defined in header** +//! +//! @code +//! #include +//! @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/math/regular/lentz_a.cpp} +//! +//! @} +//================================================================================================ +EVE_MAKE_CALLABLE(lentz_a_, lentz_a); +} + +#include diff --git a/include/eve/module/math/regular/lentz_b.hpp b/include/eve/module/math/regular/lentz_b.hpp new file mode 100644 index 0000000000..3e80bf7f37 --- /dev/null +++ b/include/eve/module/math/regular/lentz_b.hpp @@ -0,0 +1,72 @@ +//================================================================================================== +/* + EVE - Expressive Vector Engine + Copyright : EVE Project Contributors + SPDX-License-Identifier: BSL-1.0 +*/ +//================================================================================================== +#pragma once + +#include + +namespace eve +{ +//================================================================================================ +//! @addtogroup contfrac +//! @{ +//! @var lentz_b +//! @brief Implement the lentz scheme to evaluate continued fractions +//! +//! **Defined in header** +//! +//! @code +//! #include +//! @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/math/regular/lentz_b.cpp} +//! +//! @} +//================================================================================================ +EVE_MAKE_CALLABLE(lentz_b_, lentz_b); +} + +#include diff --git a/include/eve/module/math/regular/math.hpp b/include/eve/module/math/regular/math.hpp index 88eb23b2e1..a5e891b161 100644 --- a/include/eve/module/math/regular/math.hpp +++ b/include/eve/module/math/regular/math.hpp @@ -59,6 +59,8 @@ #include #include #include +#include +#include #include #include #include diff --git a/test/doc/math/regular/lentz_a.cpp b/test/doc/math/regular/lentz_a.cpp new file mode 100644 index 0000000000..a07c76cb5c --- /dev/null +++ b/test/doc/math/regular/lentz_a.cpp @@ -0,0 +1,31 @@ +#include +#include +#include +#include +#include + +using w_t = eve::wide>; + +template +struct const_fraction +{ + typedef T result_type; + + result_type operator()() + { + return T{1, 2, 3, 4, 5, 6, 7, 8}; + } +}; + + +int main() +{ + w_t z{1, 2, 3, 4, 5, 6, 7, 8}; + auto ref = (-z+eve::sqrt(eve::sqr(z)+4))/2; + std::cout << "ref constant fracs are: " << ref << std::endl; + const_fraction func; + auto gr = eve::lentz_a(func,eve::eps(eve::as()), 100); + std::cout << " constant fracs are: " << gr << std::endl; + + return 0; +} diff --git a/test/doc/math/regular/lentz_b.cpp b/test/doc/math/regular/lentz_b.cpp new file mode 100644 index 0000000000..b830af30b7 --- /dev/null +++ b/test/doc/math/regular/lentz_b.cpp @@ -0,0 +1,44 @@ +#include +#include +#include + + +template +struct const_fraction +{ + auto operator()(){ return T{1, 2, 3, 4}; } +}; + +template +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 +T mytan(T a) +{ + tan_fraction fract(a); + return a/eve::lentz_b(fract, eve::eps(eve::as>()), 100); +} + +int main() +{ + using w_t = eve::wide>; + const_fraction 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()), 100); + std::cout << " constant fracs " << gr << std::endl; + w_t z{eve::pio_4(eve::as()), eve::pio_3(eve::as()), 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; +} diff --git a/test/unit/module/math/lentz_a.cpp b/test/unit/module/math/lentz_a.cpp new file mode 100644 index 0000000000..d655f45f62 --- /dev/null +++ b/test/unit/module/math/lentz_a.cpp @@ -0,0 +1,42 @@ +//================================================================================================== +/** + EVE - Expressive Vector Engine + Copyright : EVE Project Contributors + SPDX-License-Identifier: BSL-1.0 +**/ +//================================================================================================== +#include "test.hpp" +#include + +template +struct obj_const_fraction +{ + U z; + obj_const_fraction(U v) : z(v){}; + auto operator()() { return z; } +}; + +//================================================================================================== +//== lentz_a tests +//================================================================================================== +TTS_CASE_WITH("Check behavior of lentz_a on reals", + eve::test::simd::ieee_reals, + tts::generate(tts::randoms(1.0, 4.0) + ) + ) + (T const& a1) +{ + using u_t = eve::underlying_type_t; + auto eps = eve::eps(eve::as()); + size_t m = 1000u; + double tol = sizeof(u_t) == 4 ? 1.0e-2 : 1.0e-7; + using eve::lentz_a; + { + // frac with constant coefs + auto ref = (-a1+eve::sqrt(eve::sqr(a1)+4))/2; + obj_const_fraction f1(a1); + TTS_RELATIVE_EQUAL(lentz_a(f1, eps, m), ref, tol); + auto f2 = [a1](){return a1; }; + TTS_RELATIVE_EQUAL(lentz_a(f2, eps, m), ref, tol); + } +}; diff --git a/test/unit/module/math/lentz_b.cpp b/test/unit/module/math/lentz_b.cpp new file mode 100644 index 0000000000..2285d44fda --- /dev/null +++ b/test/unit/module/math/lentz_b.cpp @@ -0,0 +1,52 @@ +//================================================================================================== +/** + EVE - Expressive Vector Engine + Copyright : EVE Project Contributors + SPDX-License-Identifier: BSL-1.0 +**/ +//================================================================================================== +#include "test.hpp" +#include + +template +struct obj_const_fraction +{ + U z; + obj_const_fraction(U v) : z(v){}; + auto operator()() { return z; } +}; + +//================================================================================================== +//== lentz_b tests +//================================================================================================== +TTS_CASE_WITH("Check behavior of lentz_b on reals", + eve::test::simd::ieee_reals, + tts::generate(tts::randoms(0.0, eve::pio_4) + ,tts::randoms(1.0, 4.0) + ) + ) + (T const& z , T const& a1) +{ + using u_t = eve::underlying_type_t; + auto eps = eve::eps(eve::as()); + size_t m = 1000u; + double tol = sizeof(u_t) == 4 ? 1.0e-2 : 1.0e-7; + using eve::lentz_b; + { + // frac with constant coefs + auto ref = (a1+eve::sqrt(eve::sqr(a1)+4))/2; + obj_const_fraction f1(a1); + TTS_RELATIVE_EQUAL(lentz_b(f1, eps, m), ref, tol); + auto f2 = [a1](){return a1; }; + TTS_RELATIVE_EQUAL(lentz_b(f2, eps, m), ref, tol); + } + { + // tan test with lambda + auto ref = eve::tan(z); + auto a = -eve::sqr(z); + auto b = eve::mone(eve::as(z)); + auto func1 = [a, &b](){ b+= 2; return kumi::tuple{a, b}; }; + auto f = eve::lentz_b(func1, eps, m); + TTS_RELATIVE_EQUAL(z/f, ref, tol); + } +};