Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cmath>
       7             : #include <numeric>
       8             : #include <utility>
       9             : #include <vector>
      10             : 
      11             : #include "DataStructures/DataVector.hpp"
      12             : #include "Utilities/ErrorHandling/Error.hpp"
      13             : #include "Utilities/MakeWithValue.hpp"
      14             : 
      15             : namespace Spectral {
      16             : /// Clenshaw's algorithm for evaluating linear combinations of
      17             : /// basis functions obeying a given three term recurrence relation.
      18             : /// See numerical recipes 3rd edition 5.4.2 for details.
      19             : /// `alpha` and `beta` define the recursion
      20             : /// \f$ f_n(x) = \alpha f_{n-1}(x) + \beta f_{n-2}(x) \f$
      21             : /// `f0_of_x` and `f1_of_x` are \f$ f_0(x)\f$ and \f$f_1(x)\f$ respectively
      22             : template <typename CoefficientDataType, typename DataType>
      23           1 : DataType evaluate_clenshaw(const std::vector<CoefficientDataType>& coefficients,
      24             :                            const DataType& alpha, const DataType& beta,
      25             :                            const DataType& f0_of_x, const DataType& f1_of_x) {
      26             :   if (coefficients.empty()) {
      27             :     return make_with_value<DataType>(f0_of_x, 0.0);
      28             :   }
      29             : 
      30             :   DataType y_upper{make_with_value<DataType>(f0_of_x, 0.0)};  // y_{N+2} = 0
      31             :   DataType new_y_lower;
      32             :   // Do the loop to calculate y_1 and y_2 (the accumulate function leaves
      33             :   // y_lower_final = y_1, and y_upper = y_2)
      34             :   const DataType y_lower_final =
      35             :       std::accumulate(coefficients.rbegin(), coefficients.rend() - 1,
      36             :                       make_with_value<DataType>(f0_of_x, 0.0),
      37             :                       [&y_upper, &alpha, &beta, &new_y_lower](
      38             :                           const auto& y_lower, const auto& c_k) {
      39             :                         new_y_lower = c_k + alpha * y_lower + beta * y_upper;
      40             :                         y_upper = y_lower;
      41             :                         return new_y_lower;
      42             :                       });
      43             :   return f0_of_x * coefficients[0] + f1_of_x * y_lower_final +
      44             :          beta * f0_of_x * y_upper;
      45             : }
      46             : 
      47             : }  // namespace Spectral