include/dfe/dfe_poly.hpp

0001 // SPDX-License-Identifier: MIT
0002 // Copyright 2018-2019 Moritz Kiehn
0003 //
0004 // Permission is hereby granted, free of charge, to any person obtaining a copy
0005 // of this software and associated documentation files (the "Software"), to deal
0006 // in the Software without restriction, including without limitation the rights
0007 // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
0008 // copies of the Software, and to permit persons to whom the Software is
0009 // furnished to do so, subject to the following conditions:
0010 //
0011 // The above copyright notice and this permission notice shall be included in
0012 // all copies or substantial portions of the Software.
0013 //
0014 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
0015 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0016 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
0017 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0018 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
0019 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
0020 // SOFTWARE.
0021
0022 /// \file
0023 /// \brief   Efficient evaluation of polynomial functions
0024 /// \author  Moritz Kiehn <msmk@cern.ch>
0025 /// \date    2018-02-26
0026
0027 #include <array>
0028 #include <iterator>
0029 #include <type_traits>
0030 #include <utility>
0031
0032 #pragma once
0033
0034 namespace dfe {
0035
0036 /// Evaluate a polynomial of arbitrary order.
0037 ///
0038 /// \param x      Where to evaluate the polynomial.
0039 /// \param coeffs ReversibleContainer with n+1 coefficients.
0040 ///
0041 /// The coefficients must be given in increasing order. E.g. a second order
0042 /// polynomial has three coefficients c0, c1, c2 that define the function
0043 ///
0044 ///     f(x) = c0 + c1*x + c2*x^2
0045 ///
0046 template<typename T, typename Container>
0047 constexpr T
0048 polynomial_val(const T& x, const Container& coeffs) {
0049   // Use Horner's method to evaluate polynomial, i.e. expand
0050   //   f(x) = c0 + c1*x + c2*x^2 + c3*x^3
0051   // to the following form
0052   //   f(x) = c0 + x * (c1 + x * (c2 + x * c3))
0053   // to allow iterative computation with minimal number of operations
0054   T value = x; // make sure dynamic-sized types, e.g. std::valarray, work
0055   value = 0;
0056   for (auto c = std::rbegin(coeffs); c != std::rend(coeffs); ++c) {
0057     value = *c + x * value;
0058   }
0059   return value;
0060 }
0061
0062 /// Evaluate the value and the derivative of a polynomial of arbitrary order.
0063 ///
0064 /// \param x      Where to evaluate the derivative.
0065 /// \param coeffs ReversibleContainer with n+1 coefficients.
0066 ///
0067 /// The coefficients must be given in increasing order. E.g. a second order
0068 /// polynomial has three coefficients c0, c1, c2 that define the function
0069 ///
0070 ///     f(x) = c0 + c1*x + c2*x^2
0071 ///
0072 /// and the derivative
0073 ///
0074 ///     df(x)/dx = 0 + c1 + 2*c2*x
0075 ///
0076 template<typename T, typename Container>
0077 constexpr std::pair<T, T>
0078 polynomial_valder(const T& x, const Container& coeffs) {
0079   // Use Horner's method to evaluate polynomial and its derivative at the
0080   // same time
0081   T q = x; // make sure dynamic-sized types, e.g. std::valarray, work
0082   T p = x;
0083   q = 0;
0084   p = 0;
0085   for (auto c = std::rbegin(coeffs); c != std::rend(coeffs); ++c) {
0086     q = p + x * q;
0087     p = *c + x * p;
0088   }
0089   return {p, q};
0090 }
0091
0092 /// Evaluate the the derivative of a polynomial of arbitrary order.
0093 ///
0094 /// \param x      Where to evaluate the derivative.
0095 /// \param coeffs ReversibleContainer with n+1 coefficients.
0096 ///
0097 /// The coefficients must be given in increasing order. E.g. a second order
0098 /// polynomial has three coefficients c0, c1, c2 that define the derivative
0099 ///
0100 ///     df(x)/dx = 0 + c1 + 2*c2*x
0101 ///
0102 template<typename T, typename Container>
0103 constexpr T
0104 polynomial_der(const T& x, const Container& coeffs) {
0105   return polynomial_valder(x, coeffs).second;
0106 }
0107
0108 /// Evaluate a polynomial with an order fixed at compile time.
0109 template<typename T, typename U>
0110 constexpr auto
0111 polynomial_val(const T& x, std::initializer_list<U> coeffs) {
0112   return polynomial_val<T, std::initializer_list<U>>(x, coeffs);
0113 }
0114
0115 /// Evaluate the derivative of a polynomial with an order fixed at compile time.
0116 template<typename T, typename U>
0117 constexpr auto
0118 polynomial_der(const T& x, std::initializer_list<U> coeffs) {
0119   return polynomial_der<T, std::initializer_list<U>>(x, coeffs);
0120 }
0121
0122 /// Evaluate the derivative of a polynomial with an order fixed at compile time.
0123 template<typename T, typename U>
0124 constexpr auto
0125 polynomial_valder(const T& x, std::initializer_list<U> coeffs) {
0126   return polynomial_valder<T, std::initializer_list<U>>(x, coeffs);
0127 }
0128
0129 } // namespace dfe