forked from sPHENIX-Collaboration/acts
/
dfe_poly.hpp
129 lines (118 loc) · 4.43 KB
/
dfe_poly.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
// SPDX-License-Identifier: MIT
// Copyright 2018-2019 Moritz Kiehn
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
/// \file
/// \brief Efficient evaluation of polynomial functions
/// \author Moritz Kiehn <msmk@cern.ch>
/// \date 2018-02-26
#include <array>
#include <iterator>
#include <type_traits>
#include <utility>
#pragma once
namespace dfe {
/// Evaluate a polynomial of arbitrary order.
///
/// \param x Where to evaluate the polynomial.
/// \param coeffs ReversibleContainer with n+1 coefficients.
///
/// The coefficients must be given in increasing order. E.g. a second order
/// polynomial has three coefficients c0, c1, c2 that define the function
///
/// f(x) = c0 + c1*x + c2*x^2
///
template<typename T, typename Container>
constexpr T
polynomial_val(const T& x, const Container& coeffs) {
// Use Horner's method to evaluate polynomial, i.e. expand
// f(x) = c0 + c1*x + c2*x^2 + c3*x^3
// to the following form
// f(x) = c0 + x * (c1 + x * (c2 + x * c3))
// to allow iterative computation with minimal number of operations
T value = x; // make sure dynamic-sized types, e.g. std::valarray, work
value = 0;
for (auto c = std::rbegin(coeffs); c != std::rend(coeffs); ++c) {
value = *c + x * value;
}
return value;
}
/// Evaluate the value and the derivative of a polynomial of arbitrary order.
///
/// \param x Where to evaluate the derivative.
/// \param coeffs ReversibleContainer with n+1 coefficients.
///
/// The coefficients must be given in increasing order. E.g. a second order
/// polynomial has three coefficients c0, c1, c2 that define the function
///
/// f(x) = c0 + c1*x + c2*x^2
///
/// and the derivative
///
/// df(x)/dx = 0 + c1 + 2*c2*x
///
template<typename T, typename Container>
constexpr std::pair<T, T>
polynomial_valder(const T& x, const Container& coeffs) {
// Use Horner's method to evaluate polynomial and its derivative at the
// same time
T q = x; // make sure dynamic-sized types, e.g. std::valarray, work
T p = x;
q = 0;
p = 0;
for (auto c = std::rbegin(coeffs); c != std::rend(coeffs); ++c) {
q = p + x * q;
p = *c + x * p;
}
return {p, q};
}
/// Evaluate the the derivative of a polynomial of arbitrary order.
///
/// \param x Where to evaluate the derivative.
/// \param coeffs ReversibleContainer with n+1 coefficients.
///
/// The coefficients must be given in increasing order. E.g. a second order
/// polynomial has three coefficients c0, c1, c2 that define the derivative
///
/// df(x)/dx = 0 + c1 + 2*c2*x
///
template<typename T, typename Container>
constexpr T
polynomial_der(const T& x, const Container& coeffs) {
return polynomial_valder(x, coeffs).second;
}
/// Evaluate a polynomial with an order fixed at compile time.
template<typename T, typename U>
constexpr auto
polynomial_val(const T& x, std::initializer_list<U> coeffs) {
return polynomial_val<T, std::initializer_list<U>>(x, coeffs);
}
/// Evaluate the derivative of a polynomial with an order fixed at compile time.
template<typename T, typename U>
constexpr auto
polynomial_der(const T& x, std::initializer_list<U> coeffs) {
return polynomial_der<T, std::initializer_list<U>>(x, coeffs);
}
/// Evaluate the derivative of a polynomial with an order fixed at compile time.
template<typename T, typename U>
constexpr auto
polynomial_valder(const T& x, std::initializer_list<U> coeffs) {
return polynomial_valder<T, std::initializer_list<U>>(x, coeffs);
}
} // namespace dfe