| 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 <polyfloat/module/core/sum_of_prod.hpp> | ||
| 15 | |||
| 16 | namespace plf | ||
| 17 | { | ||
| 18 | |||
| 19 | template<typename Options> | ||
| 20 | struct neville_t | ||
| 21 | : plf::strict_tuple_callable<neville_t, eve::neville_t, Options, kahan_option, raw_option, pedantic_option> | ||
| 22 | { | ||
| 23 | template<typename... Ts> struct result : as_polyfloat_like<Ts...> | ||
| 24 | { | ||
| 25 | }; | ||
| 26 | |||
| 27 | template<concepts::polyfloat_like... Ts> | ||
| 28 | requires(eve::same_lanes_or_scalar<Ts...>) | ||
| 29 | 18 | POLYFLOAT_ABI typename result<Ts...>::type constexpr operator()(Ts... ts) const noexcept | |
| 30 | { | ||
| 31 | 18 | return POLYFLOAT_CALL(ts...); | |
| 32 | } | ||
| 33 | |||
| 34 | 18 | POLYFLOAT_CALLABLE_OBJECT(neville_t, neville_); | |
| 35 | }; | ||
| 36 | //====================================================================================================================== | ||
| 37 | //! @addtogroup functions | ||
| 38 | //! @{ | ||
| 39 | //! @var neville | ||
| 40 | //! @brief `callable` computing the value of the interpolation polynomial \f$p\f$ of degree n, satisfying: | ||
| 41 | //! \f$p(x_i) = y_i\f$ for \f$i = 0..n\f$ | ||
| 42 | //! | ||
| 43 | //! @groupheader{Header file} | ||
| 44 | //! | ||
| 45 | //! @code | ||
| 46 | //! #include <eve/module/core.hpp> | ||
| 47 | //! @endcode | ||
| 48 | //! | ||
| 49 | //! @groupheader{Callable Signatures} | ||
| 50 | //! | ||
| 51 | //! @code | ||
| 52 | //! namespace eve | ||
| 53 | //! { | ||
| 54 | //! // Regular overloads | ||
| 55 | //! constexpr auto neville(auto x, auto const& xs..., auto const& ys...) noexcept; // 1 | ||
| 56 | //! constexpr auto neville(auto x, eve::non_empty_product_type auto const& xx, eve::non_empty_product_type auto const& yy) noexcept; // 2 | ||
| 57 | //! } | ||
| 58 | //! @endcode | ||
| 59 | //! | ||
| 60 | //! **Parameters** | ||
| 61 | //! | ||
| 62 | //! * `xs ...`, `ys ...`: [real values](@ref eve::floating_value) arguments | ||
| 63 | //! * `xx`, `yy` : kumi::tuple of arguments | ||
| 64 | //! | ||
| 65 | //! **Return value** | ||
| 66 | //! | ||
| 67 | //! * 1. Return the value at x of the polynomial \f$p\f$ of degree n, satisfying \f$p(x_i) = y_i\f$ for \f$i = 0..n\f$, where n+1 is the common number of `xs`, `ys`. | ||
| 68 | //! * 2. Same as 1., using the tuple elements. | ||
| 69 | //! | ||
| 70 | //! @groupheader{External references} | ||
| 71 | //! * [Wikipedia neville's algorithm](https://en.wikipedia.org/wiki/Neville%27s_algorithm) | ||
| 72 | //! | ||
| 73 | //! @groupheader{Example} | ||
| 74 | //! | ||
| 75 | //! @godbolt{doc/core/neville.cpp} | ||
| 76 | //====================================================================================================================== | ||
| 77 | |||
| 78 | inline constexpr auto neville = eve::functor<neville_t>; | ||
| 79 | //====================================================================================================================== | ||
| 80 | //! @} | ||
| 81 | //====================================================================================================================== | ||
| 82 | } | ||
| 83 | |||
| 84 | namespace plf::_ | ||
| 85 | { | ||
| 86 | |||
| 87 | template<typename X, typename... XsYs, eve::callable_options O> | ||
| 88 | 18 | constexpr auto neville_(POLYFLOAT_DELAY(), O const& o, X x, XsYs... xsys) noexcept | |
| 89 | requires(sizeof...(XsYs) % 2 == 0) | ||
| 90 | { | ||
| 91 | using t_t = as_polyfloat_like_t<X, XsYs...>; | ||
| 92 | 18 | constexpr auto siz = sizeof...(XsYs); | |
| 93 | 18 | constexpr auto siz_2 = siz / 2; | |
| 94 | if constexpr (siz == 0) return eve::zero(eve::as<X>()); | ||
| 95 | else | ||
| 96 | { | ||
| 97 | 44 | auto cvt = [](auto a) { return plf::convert(a, eve::as<eve::element_type_t<t_t>>()); }; | |
| 98 | 18 | auto xsyst = eve::zip(t_t(cvt(xsys))...); | |
| 99 | if constexpr (siz == 2) return get<0>(xsyst); | ||
| 100 | else if constexpr (siz == 4) | ||
| 101 | { | ||
| 102 | 10 | auto x0 = get<0>(xsyst); | |
| 103 | 10 | auto x1 = get<1>(xsyst); | |
| 104 | 10 | auto y0 = get<2>(xsyst); | |
| 105 | 10 | auto y1 = get<3>(xsyst); | |
| 106 | 10 | return plf::sum_of_prod[o]((x - x1), y0, (x0 - x), y1) / (x0 - x1); | |
| 107 | } | ||
| 108 | else | ||
| 109 | { | ||
| 110 | 8 | std::array<t_t, siz> xy{t_t(xsys)...}; | |
| 111 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 4 times.
|
24 | for (size_t k = 1; k < siz_2; ++k) |
| 112 | { | ||
| 113 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 8 times.
|
40 | for (size_t i = 0, is2 = siz_2; i < siz_2 - k; ++i, ++is2) |
| 114 | { | ||
| 115 | 24 | xy[is2] = plf::sum_of_prod[o]((x - xy[i + k]), xy[is2], (xy[i] - x), xy[is2 + 1]) / (xy[i] - xy[i + k]); | |
| 116 | } | ||
| 117 | } | ||
| 118 | 8 | return xy[siz_2]; | |
| 119 | } | ||
| 120 | } | ||
| 121 | } | ||
| 122 | } | ||
| 123 |