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] = factorial(i + D, i + 1) * 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 + d <= l; k += 2u) {
0107     sum += pow(x, k - d) * coeff(l, k) *
0108            (d > 0u ? factorial(k - d, k - d + 1u) : 1u);
0109   }
0110   return sum;
0111 }
0112 }  // namespace Legendre
0113
0114 /// @brief Implementation of the chebychev polynomials of the first & second kind
0115 ///        (c.f. https://en.wikipedia.org/wiki/Chebyshev_polynomials)
0116 namespace Chebychev {
0117 //// @brief Calculates the k-th coefficient of the Chebychev polynomial of the first kind (T_{n})
0118 /// @param n: Order of the polynomial
0119 /// @param k: Coefficient mapped to the monomial n-2k
0120 constexpr double coeffTn(const unsigned n, const unsigned k) {
0121   assert(n >= 2 * k);
0122   if (n == 0) {
0123     return 1.;
0124   }
0125   const double sign = (k % 2 == 1 ? -1. : 1.);
0126   const double t_k = sign * static_cast<double>(factorial(n - k - 1)) /
0127                      static_cast<double>(factorial(k) * factorial(n - 2 * k)) *
0128                      static_cast<double>(n) *
0129                      pow(2., static_cast<int>(n - 2 * k - 1));
0130   return t_k;
0131 }
0132 /// @brief Collects the coefficients of the n-th Chebychev polynomial
0133 ///        of the first kind and returns them via an array
0134 /// @tparam Tn: Order of the Chebychev 1-st kind polynomial
0135 template <unsigned Tn>
0136 constexpr std::array<double, Tn + 1> coeffientsFirstKind() {
0137   std::array<double, Tn + 1> coeffs{filledArray<double, Tn + 1>(0.)};
0138   for (unsigned k = 0; 2 * k <= Tn; ++k) {
0139     coeffs[Tn - 2 * k] = coeffTn(Tn, k);
0140   }
0141   return coeffs;
0142 }
0143 /// @brief Evaluates the chebychev polynomial of the first kind
0144 /// @param x: Domain value to evaluate
0145 /// @param n: Order of the chebychev polynomial
0146 /// @param d: Order of the derivative
0147 constexpr double evalFirstKind(const double x, const unsigned n,
0148                                const unsigned d = 0u) {
0149   double result{0.};
0150   for (unsigned k = 0u; 2u * k + d <= n; ++k) {
0151     result += coeffTn(n, k) * pow(x, n - 2u * k - d) *
0152               (d > 0 ? factorial(n - 2u * k, n - 2u * k - d + 1u) : 1);
0153   }
0154   return result;
0155 }
0156 /// @brief Calculates the k-th coefficient of the Chebychev polynomial of the second kind (U_{n})
0157 /// @param n: Order of the polynomial
0158 /// @param k: Coefficient mapped to the monomial n-2k
0159 constexpr double coeffUn(const unsigned n, const unsigned k) {
0160   assert(n >= 2u * k);
0161   if (n == 0u) {
0162     return 1.;
0163   }
0164   const double sign = (k % 2 == 1u ? -1. : 1.);
0165   const double u_k = sign * binomial(n - k, k) * pow(2., n - 2u * k);
0166   return u_k;
0167 }
0168
0169 /// @brief Collects the coefficients of the n-th Chebychev polynomial
0170 ///        of the second kind and returns them via an array
0171 /// @tparam Un: Order of the Chebychev 2-nd kind polynomial
0172 template <unsigned Un>
0173 constexpr std::array<double, Un + 1> coeffientsSecondKind() {
0174   std::array<double, Un + 1> coeffs{filledArray<double, Un + 1>(0.)};
0175   for (unsigned k = 0u; 2u * k <= Un; ++k) {
0176     coeffs[Un - 2u * k] = coeffUn(Un, k);
0177   }
0178   return coeffs;
0179 }
0180 /// @brief Evaluates the chebychev polynomial of the second kind
0181 /// @param x: Domain value to evaluate
0182 /// @param n: Order of the chebychev polynomial
0183 /// @param d: Order of the derivative
0184 constexpr double evalSecondKind(const double x, const unsigned n,
0185                                 const unsigned d) {
0186   double result{0.};
0187   for (unsigned k = 0u; 2u * k + d <= n; ++k) {
0188     result += coeffUn(n, k) * pow(x, n - 2u * k - d) *
0189               (d > 0u ? factorial(n - 2u * k, n - 2u * k - d + 1u) : 1u);
0190   }
0191   return result;
0192 }
0193 }  // namespace Chebychev
0194
0195 /// @brief Helper macros to setup the evaluation of the n-th orthogonal
0196 ///        polynomial and of its derivatives. The  coefficients of the
0197 ///        n-th polynomial and of the first two derivatives
0198 ///        are precalculated at compile time. For higher order derivatives
0199 ///        a runtime evaluation function is used
0200 /// @param order: Explicit order of the orthogonal polynomial
0201 /// @param derivative: Requested derivative from the function call
0202 /// @param coeffGenerator: Templated function which returns the coefficients
0203 ///                        as an array
0204 #define SETUP_POLY_ORDER(order, derivative, coeffGenerator) \
0205   case order: {                                             \
0206     constexpr auto polyCoeffs = coeffGenerator<order>();    \
0207     switch (derivative) {                                   \
0208       case 0u:                                              \
0209         return derivativeSum<0>(x, polyCoeffs);             \
0210       case 1u:                                              \
0211         return derivativeSum<1>(x, polyCoeffs);             \
0212       case 2u:                                              \
0213         return derivativeSum<2>(x, polyCoeffs);             \
0214       default:                                              \
0215         break;                                              \
0216     }                                                       \
0217     break;                                                  \
0218   }
0219 /// @brief Helper macro to write a function to evaluate an orthogonal
0220 ///        polynomial. The first 10 polynomial terms together with their
0221 ///        first two derivatives are compiled as constexpr expressions
0222 ///        for the remainders the fall back runTime evaluation function is used
0223 /// @param polyFuncName: Name of the final function
0224 /// @param coeffGenerator: Name of the templated function over the polynomial order
0225 ///                        generating the coefficients of the n-th orthogonal
0226 ///                        polynomial
0227 /// @param runTimeEval: Evaluation function with the signature
0228 ///                       double runTimeEval(const double x,
0229 ///                                          const unsigned order,
0230 ///                                          onst unsigned derivative)
0231 ///                     which is used as fallback if higher order derivatives or
0232 ///                     higher order polynomials are requested from the user.
0233 #define EVALUATE_POLYNOMIAL(polyFuncName, coeffGenerator, runTimeEval) \
0234   constexpr double polyFuncName(double x, unsigned order,              \
0235                                 unsigned derivative = 0) {             \
0236     switch (order) {                                                   \
0237       SETUP_POLY_ORDER(0u, derivative, coeffGenerator)                 \
0238       SETUP_POLY_ORDER(1u, derivative, coeffGenerator)                 \
0239       SETUP_POLY_ORDER(2u, derivative, coeffGenerator)                 \
0240       SETUP_POLY_ORDER(3u, derivative, coeffGenerator)                 \
0241       SETUP_POLY_ORDER(4u, derivative, coeffGenerator)                 \
0242       SETUP_POLY_ORDER(5u, derivative, coeffGenerator)                 \
0243       SETUP_POLY_ORDER(6u, derivative, coeffGenerator)                 \
0244       SETUP_POLY_ORDER(7u, derivative, coeffGenerator)                 \
0245       SETUP_POLY_ORDER(8u, derivative, coeffGenerator)                 \
0246       SETUP_POLY_ORDER(9u, derivative, coeffGenerator)                 \
0247       SETUP_POLY_ORDER(10u, derivative, coeffGenerator)                \
0248       default:                                                         \
0249         break;                                                         \
0250     }                                                                  \
0251     return runTimeEval(x, order, derivative);                          \
0252   }
0253
0254 EVALUATE_POLYNOMIAL(chebychevPolyTn, Chebychev::coeffientsFirstKind,
0255                     Chebychev::evalFirstKind);
0256 EVALUATE_POLYNOMIAL(chebychevPolyUn, Chebychev::coeffientsSecondKind,
0257                     Chebychev::evalSecondKind);
0258 EVALUATE_POLYNOMIAL(legendrePoly, Legendre::coefficients, Legendre::evaluate);
0259 #undef EVALUATE_POLYNOMIAL
0260 #undef SETUP_POLY_ORDER
0261
0262 }  // namespace Acts::detail