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