GCC Code Coverage Report

Line	Branch	Exec	Source
1			//======================================================================================================================
2			/*
3			POLYFLOAT - Extended precision floating points
4			Copyright : POLYFLOAT Contributors & Maintainers
5			SPDX-License-Identifier: BSL-1.0
6			*/
7			//======================================================================================================================
8			#pragma once
9
10			#include <polyfloat/details/callable.hpp>
11			#include <polyfloat/types/concepts.hpp>
12			#include <polyfloat/types/traits.hpp>
13			#include <type_traits>
14			#include <eve/concept/range.hpp>
15			#include <polyfloat/module/core/convert.hpp>
16			#include <polyfloat/module/core/two_add.hpp>
17			#include <polyfloat/module/core/two_prod.hpp>
18			#include <eve/module/math.hpp>
19
20			namespace plf
21			{
22
23			template<typename Options>
24			struct reverse_horner_t : plf::strict_tuple_callable<reverse_horner_t,
25			eve::reverse_horner_t,
26			Options,
27			kahan_option,
28			raw_option,
29			pedantic_option>
30			{
31			template<typename... Ts> struct result : as_polyfloat_like<Ts...>
32			{
33			};
34
35			template<concepts::polyfloat_like... Ts>
36			requires(eve::same_lanes_or_scalar<Ts...>)
37		10	POLYFLOAT_ABI typename result<Ts...>::type constexpr operator()(Ts... ts) const noexcept
38			{
39		10	return POLYFLOAT_CALL(ts...);
40			}
41
42			template<concepts::polyfloat_like X, eve::non_empty_product_type Tup>
43			requires(eve::same_lanes_or_scalar_tuple<Tup> && !concepts::polyfloat_like<Tup>)
44		3360	POLYFLOAT_ABI constexpr as_polyfloat_like_t<X, Tup> operator()(X const& x, Tup const& t) const noexcept
45			requires(kumi::size_v<Tup> >= 1)
46			{
47		3360	return POLYFLOAT_CALL(x, t);
48			}
49
50		3368	POLYFLOAT_CALLABLE_OBJECT(reverse_horner_t, reverse_horner_);
51			};
52			//======================================================================================================================
53			//! @addtogroup functions
54			//! @{
55			//! @var reverse_horner
56			//! @brief Implement the reverse_horner scheme to evaluate polynomials with coefficients
57			//! in decreasing power order
58			//!
59			//! @groupheader{Header file}
60			//!
61			//! @code
62			//! #include <polyfloat/functions.hpp>
63			//! @endcode
64			//!
65			//! @groupheader{Callable Signatures}
66			//!
67			//! @code
68			//! namespace polyfloat
69			//! {
70			//! @endcode
71			//!
72			//! @groupheader{Callable Signatures}
73			//!
74			//! @code
75			//! namespace eve
76			//! {
77			//! // Regular overloads
78			//! constexpr auto reverse_horner(floating_value auto x, value auto ...ci) noexcept; // 1
79			//! constexpr auto reverse_horner(floating_value auto x, eve::coefficients auto tci) noexcept; // 2
80			//!
81			//! // Lanes masking
82			//! constexpr auto reverse_horner[conditional_expr auto c](/*any of the above overloads*/) noexcept; // 3
83			//! constexpr auto reverse_horner[logical_value auto m](/*any of the above overloads*/) noexcept; // 3
84			//!
85			//! // Semantic options
86			//! constexpr auto reverse_horner[pedantic](/*any of the above overloads*/) noexcept; // 4
87			//! constexpr auto reverse_horner[kahan](/*any of the above overloads*/) noexcept; // 5
88			//! }
89			//! @endcode
90			//!
91			//! **Parameters**
92			//!
93			//! * `x`: [evaluation point floating value](@ref eve::floating_value) arguments.
94			//! * `ci...`: [floating values](@ref eve::floating_value) polynom coefficients in decreasing power order,
95			//! Note that the values of the `ci` are not necessarily floating but the non floating ones
96			//! are to be scalar
97			//! * `tci`: eve::coefficients is a specialization of `product_type` of floating values,
98			//! used to avoid possible ambiguities
99			//! * `c`: `conditional_expr` masking the operation.
100			//! * `m`: [Logical value](@ref eve::logical_value) masking the operation.
101			//!
102			//! **Return value**
103			//!
104			//! If \f$(c_i)_{0\le i\le n-1}\f$ denotes the coefficients of the polynomial by decreasing
105			//! power order, the Reverse_Horner scheme evaluates the polynom \f$p\f$ at \f$x\f$ by :
106			//! \f$\qquad\qquad\displaystyle p(x) = (((c_0x+c_1)x+ ... )x + c_{n-1})\f$
107			//!
108			//! 1. The value of the polynom at `x` is returned.
109			//! 2. Same as the call with the elements of the tuple.
110			//! 3. [The operation is performed conditionnaly](@ref conditional).
111			//! 4. `fma[pedantic]` instead of `fma` is used in internal computations.
112			//! 5. a Kahan like compensated algorithm is used to enhance accuracy.
113			//!
114			//! @note If the coefficients are simd values of width N, this means you simultaneously
115			//! compute the values of N polynomials.
116			//! * If x is scalar, the polynomials are all computed at the same point
117			//! * If x is simd, the nth polynomial is computed on the nth value of x
118			//!
119			//! @groupheader{External references}
120			//! * [Wikipedia](https://en.wikipedia.org/wiki/Reverse_Horner's_method)
121			//!
122			//!
123			//! @groupheader{Example}
124			//!
125			//! @godbolt{doc/core/reverse_horner.cpp}
126			//======================================================================================================================
127
128			inline constexpr auto reverse_horner = eve::functor<reverse_horner_t>;
129			//======================================================================================================================
130			//! @}
131			//======================================================================================================================
132			}
133
134			namespace plf::_
135			{
136
137			template<typename X, eve::callable_options O> constexpr auto reverse_horner_(POLYFLOAT_DELAY(), O const&, X) noexcept
138			{
139			return plf::zero(eve::as<X>());
140			}
141
142			template<typename X, typename C, typename... Cs, eve::callable_options O>
143		10	constexpr auto reverse_horner_(POLYFLOAT_DELAY(), O const& o, X xx, C c0, Cs... cs) noexcept
144			{
145			using r_t = as_polyfloat_like_t<X, C, Cs...>;
146		10	constexpr auto N = sizeof...(Cs);
147			if constexpr (N == 0) return plf::convert(c0, eve::as_element<r_t>{});
148			else
149			{
150		10	auto x = r_t(xx);
151			using t_t = kumi::result::fill_t<sizeof...(cs) + 1, r_t>;
152		10	t_t c{r_t{c0}, r_t{cs}...};
153		10	return reverse_horner[o](x, coefficients<t_t>(c));
154			}
155			}
156
157			template<typename X, eve::product_type Tuple, eve::callable_options O>
158		3358	constexpr auto reverse_horner_(POLYFLOAT_DELAY(), O const& o, X x, eve::coefficients<Tuple> const& tup) noexcept
159			{
160		6716	return horner[o](x, eve::coefficients(kumi::reverse(tup)));
161			}
162
163			template<typename X, eve::_::range R, eve::callable_options O>
164			constexpr auto reverse_horner_(POLYFLOAT_DELAY(), O const& o, X xx, R const& r) noexcept
165			{
166			using r_t = as_polyfloat_like<X, typename R::value_type>;
167			using u_t = eve::element_type_t<r_t>;
168			auto cvt = [](auto a) { return plf::convert(a, as<u_t>()); };
169			auto x = cvt(xx);
170			auto cur = std::rbegin(r);
171			auto first = std::rend(r);
172			if (first == cur) return zero(as<r_t>());
173			else if (std::distance(cur, first) == 1) return cvt(*cur);
174			else
175			{
176			auto dfma = fma[o];
177			auto that = cvt(zero(as<r_t>()));
178			auto step = [&](auto th, auto arg) { return dfma(x, th, arg); };
179			for (; cur != first; ++cur) that = step(that, cvt(*cur));
180			return that;
181			}
182			}
183			}
184