Utilities/detail/Polynomials.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include "Acts/Utilities/ArrayHelpers.hpp"
0012 #include "Acts/Utilities/MathHelpers.hpp"
0013
0014 #include <cassert>
0015
0016 namespace Acts::detail {
0017 /// @brief Transforms the coefficients of a n-th degree polynomial
0018 ///        into the one corresponding to its D-th derivative
0019 ///        The result are (N-D) non-vanishing coefficients
0020 /// @tparam D: Order of the derivative to calculate
0021 /// @tparam N: Order of the original polynomial
0022 /// @param coeffs: Reference to the polynomial's coefficients
0023 template <std::size_t D, std::size_t N>
0024 constexpr std::array<double, N - D> derivativeCoefficients(
0025     const std::array<double, N>& coeffs) {
0026   static_assert(N > D, "Coefficients trivially collapse to 0.");
0027   if constexpr (D == 0) {
0028     return coeffs;
0029   }
0030   std::array<double, N - D> newCoeffs{filledArray<double, N - D>(0.)};
0031   for (std::size_t i = 0; i < N - D; ++i) {
0032     newCoeffs[i] = product(i + 1ul, i + D) * coeffs[i + D];
0033   }
0034   return newCoeffs;
0035 }
0036
0037 /// @brief Evaluates a polynomial with degree (N-1) at domain value x
0038 /// @param x: Value where the polynomial is to be evaluated
0039 /// @param coeff: Polynomial coefficients, where the i-th entry
0040 ///               is associated to the x^{i} monomial
0041 template <std::size_t N>
0042 constexpr double polynomialSum(const double x,
0043                                const std::array<double, N>& coeffs) {
0044   double result{0.};
0045   double y{1.};
0046   for (unsigned k = 0; k < N; ++k) {
0047     if (abs(coeffs[k]) > std::numeric_limits<double>::epsilon()) {
0048       result += coeffs[k] * y;
0049     }
0050     y *= x;
0051   }
0052   return result;
0053 }
0054
0055 /// @brief Evaluates the D-th derivative of a polynomial with degree (N-1)
0056 ///        at domain value x
0057 /// @tparam D: Order of the derivative to calculate
0058 /// @tparam N: Order of the original polynomial
0059 /// @param x: Value where the polynomial is to be evaluated
0060 /// @param coeff: Polynomial coefficients, where the i-th entry
0061 ///               is associated to the x^{i} monomial
0062 template <std::size_t D, std::size_t N>
0063 constexpr double derivativeSum(const double x,
0064                                const std::array<double, N>& coeffs) {
0065   if constexpr (N > D) {
0066     return polynomialSum(x, derivativeCoefficients<D, N>(coeffs));
0067   }
0068   return 0.;
0069 }
0070
0071 namespace Legendre {
0072 /// @brief Calculates the n-th coefficient of the legendre polynomial series
0073 ///        (cf. https://en.wikipedia.org/wiki/Legendre_polynomials)
0074 /// @param l: Order of the legendre polynomial
0075 /// @param k: Coefficient inside the polynomial representation
0076 constexpr double coeff(const unsigned l, const unsigned k) {
0077   assert(k <= l);
0078   if ((k % 2) != (l % 2)) {
0079     return 0.;
0080   } else if (k > 1) {
0081     const double a_k = -(1. * (l - k + 2) * (l + k - 1)) /
0082                        (1. * (k * (k - 1))) * coeff(l, k - 2);
0083     return a_k;
0084   }
0085   unsigned fl = (l - l % 2) / 2;
0086   unsigned binom = binomial(l, fl) * binomial(2 * l - 2 * fl, l);
0087   return ((fl % 2) != 0 ? -1. : 1.) * pow(0.5, l) * (1. * binom);
0088 }
0089 /// @brief Assembles the coefficients of the L-th legendre polynomial
0090 ///        and returns them via an array
0091 /// @tparam L: Order of the legendre polynomial
0092 template <unsigned L>
0093 constexpr std::array<double, L + 1> coefficients() {
0094   std::array<double, L + 1> allCoeffs{filledArray<double, L + 1>(0.)};
0095   for (unsigned k = (L % 2); k <= L; k += 2) {
0096     allCoeffs[k] = coeff(L, k);
0097   }
0098   return allCoeffs;
0099 }
0100 /// @brief Evaluates the Legendre polynomial
0101 /// @param x: Domain value to evaluate
0102 /// @param l: Order of the Legendre polynomial
0103 /// @param d: Order of the derivative
0104 constexpr double evaluate(const double x, const unsigned l, unsigned d = 0u) {
0105   double sum{0.};
0106   for (unsigned k = l % 2; k <= l; k += 2u) {
0107     if (k >= d) {
0108       sum +=
0109           pow(x, k - d) * coeff(l, k) * (d > 0u ? product(1u + k - d, k) : 1u);
0110     }
0111   }
0112   return sum;
0113 }
0114 }  // namespace Legendre
0115
0116 /// @brief Implementation of the chebychev polynomials of the first & second kind
0117 ///        (c.f. https://en.wikipedia.org/wiki/Chebyshev_polynomials)
0118 namespace Chebychev {
0119 //// @brief Calculates the k-th coefficient of the Chebychev polynomial of the first kind (T_{n})
0120 /// @param n: Order of the polynomial
0121 /// @param k: Coefficient mapped to the monomial n-2k
0122 constexpr double coeffTn(const unsigned n, const unsigned k) {
0123   assert(n >= 2 * k);
0124   if (n == 0) {
0125     return 1.;
0126   }
0127   const double sign = (k % 2 == 1 ? -1. : 1.);
0128   const double t_k = sign * static_cast<double>(factorial(n - k - 1)) /
0129                      (static_cast<double>(factorial(k)) *
0130                       static_cast<double>(factorial(n - 2 * k))) *
0131                      static_cast<double>(n) *
0132                      pow(2., static_cast<int>(n - 2 * k - 1));
0133   return t_k;
0134 }
0135 /// @brief Collects the coefficients of the n-th Chebychev polynomial
0136 ///        of the first kind and returns them via an array
0137 /// @tparam Tn: Order of the Chebychev 1-st kind polynomial
0138 template <unsigned Tn>
0139 constexpr std::array<double, Tn + 1> coeffientsFirstKind() {
0140   std::array<double, Tn + 1> coeffs{filledArray<double, Tn + 1>(0.)};
0141   for (unsigned k = 0; 2 * k <= Tn; ++k) {
0142     coeffs[Tn - 2 * k] = coeffTn(Tn, k);
0143   }
0144   return coeffs;
0145 }
0146 /// @brief Evaluates the chebychev polynomial of the first kind
0147 /// @param x: Domain value to evaluate
0148 /// @param n: Order of the chebychev polynomial
0149 /// @param d: Order of the derivative
0150 constexpr double evalFirstKind(const double x, const unsigned n,
0151                                const unsigned d = 0u) {
0152   double result{0.};
0153   for (unsigned k = 0u; 2u * k + d <= n; ++k) {
0154     result += coeffTn(n, k) * pow(x, n - 2u * k - d) *
0155               (d > 0 ? product(n - 2u * k - d + 1u, n - 2u * k) : 1);
0156   }
0157   return result;
0158 }
0159 /// @brief Calculates the k-th coefficient of the Chebychev polynomial of the second kind (U_{n})
0160 /// @param n: Order of the polynomial
0161 /// @param k: Coefficient mapped to the monomial n-2k
0162 constexpr double coeffUn(const unsigned n, const unsigned k) {
0163   assert(n >= 2u * k);
0164   if (n == 0u) {
0165     return 1.;
0166   }
0167   const double sign = (k % 2 == 1u ? -1. : 1.);
0168   const double u_k = sign * binomial(n - k, k) * pow(2., n - 2u * k);
0169   return u_k;
0170 }
0171
0172 /// @brief Collects the coefficients of the n-th Chebychev polynomial
0173 ///        of the second kind and returns them via an array
0174 /// @tparam Un: Order of the Chebychev 2-nd kind polynomial
0175 template <unsigned Un>
0176 constexpr std::array<double, Un + 1> coeffientsSecondKind() {
0177   std::array<double, Un + 1> coeffs{filledArray<double, Un + 1>(0.)};
0178   for (unsigned k = 0u; 2u * k <= Un; ++k) {
0179     coeffs[Un - 2u * k] = coeffUn(Un, k);
0180   }
0181   return coeffs;
0182 }
0183 /// @brief Evaluates the chebychev polynomial of the second kind
0184 /// @param x: Domain value to evaluate
0185 /// @param n: Order of the chebychev polynomial
0186 /// @param d: Order of the derivative
0187 constexpr double evalSecondKind(const double x, const unsigned n,
0188                                 const unsigned d) {
0189   double result{0.};
0190   for (unsigned k = 0u; 2u * k + d <= n; ++k) {
0191     result += coeffUn(n, k) * pow(x, n - 2u * k - d) *
0192               (d > 0u ? product(n - 2u * k - d + 1u, n - 2u * k) : 1u);
0193   }
0194   return result;
0195 }
0196 }  // namespace Chebychev
0197
0198 /// @brief Helper macros to setup the evaluation of the n-th orthogonal
0199 ///        polynomial and of its derivatives. The  coefficients of the
0200 ///        n-th polynomial and of the first two derivatives
0201 ///        are precalculated at compile time. For higher order derivatives
0202 ///        a runtime evaluation function is used
0203 /// @param order: Explicit order of the orthogonal polynomial
0204 /// @param derivative: Requested derivative from the function call
0205 /// @param coeffGenerator: Templated function which returns the coefficients
0206 ///                        as an array
0207 #define SETUP_POLY_ORDER(order, derivative, coeffGenerator) \
0208   case order: {                                             \
0209     constexpr auto polyCoeffs = coeffGenerator<order>();    \
0210     switch (derivative) {                                   \
0211       case 0u:                                              \
0212         return derivativeSum<0>(x, polyCoeffs);             \
0213       case 1u:                                              \
0214         return derivativeSum<1>(x, polyCoeffs);             \
0215       case 2u:                                              \
0216         return derivativeSum<2>(x, polyCoeffs);             \
0217       default:                                              \
0218         break;                                              \
0219     }                                                       \
0220     break;                                                  \
0221   }
0222 /// @brief Helper macro to write a function to evaluate an orthogonal
0223 ///        polynomial. The first 10 polynomial terms together with their
0224 ///        first two derivatives are compiled as constexpr expressions
0225 ///        for the remainders the fall back runTime evaluation function is used
0226 /// @param polyFuncName: Name of the final function
0227 /// @param coeffGenerator: Name of the templated function over the polynomial order
0228 ///                        generating the coefficients of the n-th orthogonal
0229 ///                        polynomial
0230 /// @param runTimeEval: Evaluation function with the signature
0231 ///                       double runTimeEval(const double x,
0232 ///                                          const unsigned order,
0233 ///                                          onst unsigned derivative)
0234 ///                     which is used as fallback if higher order derivatives or
0235 ///                     higher order polynomials are requested from the user.
0236 #define EVALUATE_POLYNOMIAL(polyFuncName, coeffGenerator, runTimeEval) \
0237   constexpr double polyFuncName(double x, unsigned order,              \
0238                                 unsigned derivative = 0) {             \
0239     switch (order) {                                                   \
0240       SETUP_POLY_ORDER(0u, derivative, coeffGenerator)                 \
0241       SETUP_POLY_ORDER(1u, derivative, coeffGenerator)                 \
0242       SETUP_POLY_ORDER(2u, derivative, coeffGenerator)                 \
0243       SETUP_POLY_ORDER(3u, derivative, coeffGenerator)                 \
0244       SETUP_POLY_ORDER(4u, derivative, coeffGenerator)                 \
0245       SETUP_POLY_ORDER(5u, derivative, coeffGenerator)                 \
0246       SETUP_POLY_ORDER(6u, derivative, coeffGenerator)                 \
0247       SETUP_POLY_ORDER(7u, derivative, coeffGenerator)                 \
0248       SETUP_POLY_ORDER(8u, derivative, coeffGenerator)                 \
0249       SETUP_POLY_ORDER(9u, derivative, coeffGenerator)                 \
0250       SETUP_POLY_ORDER(10u, derivative, coeffGenerator)                \
0251       default:                                                         \
0252         break;                                                         \
0253     }                                                                  \
0254     return runTimeEval(x, order, derivative);                          \
0255   }
0256
0257 EVALUATE_POLYNOMIAL(chebychevPolyTn, Chebychev::coeffientsFirstKind,
0258                     Chebychev::evalFirstKind);
0259 EVALUATE_POLYNOMIAL(chebychevPolyUn, Chebychev::coeffientsSecondKind,
0260                     Chebychev::evalSecondKind);
0261 EVALUATE_POLYNOMIAL(legendrePoly, Legendre::coefficients, Legendre::evaluate);
0262 #undef EVALUATE_POLYNOMIAL
0263 #undef SETUP_POLY_ORDER
0264
0265 }  // namespace Acts::detail