GCC Code Coverage Report

Directory: include/polyfloat/
File: include/polyfloat/module/core/horner.hpp
Date: 2026-11-09 15:28:50
Exec Total Coverage
Lines: 28 28 100.0%
Functions: 180 180 100.0%
Branches: 1 1 100.0%
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 horner_t
25 : plf::strict_tuple_callable<horner_t, eve::horner_t, Options, kahan_option, raw_option, pedantic_option>
26 {
27 template<concepts::polyfloat_like... Ts>
28 // requires(eve::same_lanes_or_scalar<Ts...>)
29 3422 POLYFLOAT_ABI as_polyfloat_like_t<Ts...> constexpr operator()(Ts... ts) const noexcept
30 {
31 3422 return POLYFLOAT_CALL(ts...);
32 }
33
34 template<concepts::polyfloat_like X, eve::non_empty_product_type Tup>
35 requires(eve::same_lanes_or_scalar_tuple<Tup> && !concepts::polyfloat_like<Tup>)
36 3378 POLYFLOAT_ABI constexpr as_polyfloat_like_t<X, Tup> operator()(X const& x, Tup const& t) const noexcept
37 requires(kumi::size_v<Tup> >= 1)
38 {
39 3378 return POLYFLOAT_CALL(x, t);
40 }
41
42 10898 POLYFLOAT_CALLABLE_OBJECT(horner_t, horner_);
43 };
44 //======================================================================================================================
45 //! @addtogroup functions
46 //! @{
47 //! @var horner
48 //! @brief Implement the horner scheme to evaluate polynomials with coefficients
49 //! in decreasing power order
50 //!
51 //! @groupheader{Header file}
52 //!
53 //! @code
54 //! #include <polyfloat/functions.hpp>
55 //! @endcode
56 //!
57 //! @groupheader{Callable Signatures}
58 //!
59 //! @code
60 //! namespace polyfloat
61 //! {
62 //! @endcode
63 //!
64 //! @groupheader{Callable Signatures}
65 //!
66 //! @code
67 //! namespace eve
68 //! {
69 //! // Regular overloads
70 //! constexpr auto horner(floating_value auto x, value auto ...ci) noexcept; // 1
71 //! constexpr auto horner(floating_value auto x, eve::coefficients auto tci) noexcept; // 2
72 //!
73 //! // Lanes masking
74 //! constexpr auto horner[conditional_expr auto c](/*any of the above overloads*/) noexcept; // 3
75 //! constexpr auto horner[logical_value auto m](/*any of the above overloads*/) noexcept; // 3
76 //!
77 //! // Semantic options
78 //! constexpr auto horner[pedantic](/*any of the above overloads*/) noexcept; // 4
79 //! constexpr auto horner[kahan](/*any of the above overloads*/) noexcept; // 5
80 //! }
81 //! @endcode
82 //!
83 //! **Parameters**
84 //!
85 //! * `x`: [evaluation point floating value](@ref eve::floating_value) arguments.
86 //! * `ci...`: [floating values](@ref eve::floating_value) polynom coefficients in decreasing power order,
87 //! Note that the values of the `ci` are not necessarily floating but the non floating ones
88 //! are to be scalar
89 //! * `tci`: eve::coefficients is a specialization of `product_type` of floating values,
90 //! used to avoid possible ambiguities
91 //! * `c`: `conditional_expr` masking the operation.
92 //! * `m`: [Logical value](@ref eve::logical_value) masking the operation.
93 //!
94 //! **Return value**
95 //!
96 //! If \f$(c_i)_{0\le i\le n-1}\f$ denotes the coefficients of the polynomial by decreasing
97 //! power order, the Horner scheme evaluates the polynom \f$p\f$ at \f$x\f$ by :
98 //! \f$\qquad\qquad\displaystyle p(x) = (((c_0x+c_1)x+ ... )x + c_{n-1})\f$
99 //!
100 //! 1. The value of the polynom at `x` is returned.
101 //! 2. Same as the call with the elements of the tuple.
102 //! 3. [The operation is performed conditionnaly](@ref conditional).
103 //! 4. `fma[pedantic]` instead of `fma` is used in internal computations.
104 //! 5. a Kahan like compensated algorithm is used to enhance accuracy.
105 //!
106 //! @note If the coefficients are simd values of width N, this means you simultaneously
107 //! compute the values of N polynomials.
108 //! * If x is scalar, the polynomials are all computed at the same point
109 //! * If x is simd, the nth polynomial is computed on the nth value of x
110 //!
111 //! @groupheader{External references}
112 //! * [Wikipedia](https://en.wikipedia.org/wiki/Horner's_method)
113 //!
114 //!
115 //! @groupheader{Example}
116 //!
117 //! @godbolt{doc/core/horner.cpp}
118 //======================================================================================================================
119
120 inline constexpr auto horner = eve::functor<horner_t>;
121 //======================================================================================================================
122 //! @}
123 //======================================================================================================================
124 }
125
126 namespace plf::_
127 {
128
129 template<typename X, eve::callable_options O> constexpr auto horner_(POLYFLOAT_DELAY(), O const&, X) noexcept
130 {
131 return plf::zero(eve::as<X>());
132 }
133
134 template<typename X, typename C, typename... Cs, eve::callable_options O>
135 3422 constexpr auto horner_(POLYFLOAT_DELAY(), O const&, X xx, C c, Cs... cs) noexcept
136 {
137 using r_t = as_polyfloat_like_t<X, C, Cs...>;
138 // using u_t = eve::as_element<r_t>; //r_t::ptype_t;
139 8 auto cvt = [](auto a) { return plf::convert(a, eve::as_element<r_t>{}); };
140 3422 constexpr auto N = sizeof...(Cs);
141 if constexpr (N == 0) return plf::convert(c, eve::as_element<r_t>{});
142 else if constexpr (O::contains(kahan))
143 {
144 using a_t = std::array<r_t, N>;
145 4 a_t err;
146 4 auto i = 0;
147 4 auto s = cvt(c);
148 4 auto x = cvt(xx);
149 8 auto step = [&s, &err, x, &i](auto a) {
150 8 auto [pi, epi] = plf::two_prod(s, x);
151 8 auto [si, esi] = plf::two_add(pi, a);
152 4 s = si;
153 4 err[i] = epi + esi;
154 4 ++i;
155 4 return s;
156 };
157 4 ((s = step(cvt(cs))), ...);
158 using tup_t = kumi::result::generate_t<N, decltype([](std::size_t) { return r_t(); })>;
159 4 auto t = std::bit_cast<tup_t, a_t>(err);
160 4 return s + plf::horner(x, coefficients(t));
161 }
162 else
163 {
164 3418 auto x = plf::convert(xx, eve::as_element<r_t>{});
165 3418 r_t that(c);
166 3418 ((that = fma[pedantic](that, x, plf::convert(cs, eve::as_element<r_t>{}))), ...);
167 3418 return that;
168 }
169 }
170
171 template<typename X, eve::product_type Tuple, eve::callable_options O>
172 3376 constexpr auto horner_(POLYFLOAT_DELAY(), O const& o, X x, eve::coefficients<Tuple> const& tup) noexcept
173 {
174 if constexpr (Tuple::size() == 0) return eve::zero(as(x));
175 else
176 {
177 using r_t = as_polyfloat_like_t<X, Tuple>;
178
1/1
✓ Branch 2 taken 2050 times.
 9465 return kumi::apply([&](auto... m) { return plf::horner[o](x, convert(m, eve::as_element<r_t>())...); }, tup);
179 }
180 }
181
182 template<typename X, eve::_::range R, eve::callable_options O>
183 constexpr auto horner_(POLYFLOAT_DELAY(), O const& o, X xx, R const& r) noexcept
184 {
185 using r_t = as_polyfloat_like<X, typename R::value_type>;
186 using u_t = eve::element_type_t<r_t>;
187 auto cvt = [](auto a) { return plf::convert(a, as<u_t>()); };
188 auto x = cvt(xx);
189 auto cur = std::begin(r);
190 auto last = std::end(r);
191 if (last == cur) return zero(as<r_t>());
192 else if (std::distance(cur, last) == 1) return cvt(*cur);
193 else
194 {
195 using std::advance;
196 auto that = cv_t(*cur);
197 auto step = [&](auto th, auto arg) { return fma[pedantic][o](x, th, arg); };
198 for (advance(cur, 1); cur != last; advance(cur, 1)) that = step(that, cv_t(*cur));
199 return that;
200 }
201 }
202 }
203