legendre_polynomials.cc

Go to the documentation of this file.

#include "math/legendre_polynomials.h"
 
#include <cmath>
#include <iomanip>
#include <limits>
#include <numeric>
#include <sstream>
#include <stdexcept>
 
#include "utils/constants.h"
#include "utils/misc.h"
 
namespace hummingbird {
double LegendrePolynomial(const int n, const double x) {
  switch (n) {
    case 0:
      return 1.0;
 
    case 1:
      return x;
 
    default:
      double previous = 1.0;
      double current = x;
      for (double j = 1; j < n; j++) {
        double next =
            ((2.0 * j + 1.0) * x * current - j * previous) / (j + 1.0);
        previous = current;
        current = next;
      }
      return current;
  }
}
 
double LegendrePolynomialPrime(const int n, const double x) {
  if (DoubleEqual(x, -1.0))
    return std::pow(-1, n - 1) * n * (n + 1) / 2.0;
  else if (DoubleEqual(x, 1.0))
    return n * (n + 1) / 2.0;
 
  switch (n) {
    case 0:
      return 0;
 
    case 1:
      return 1.0;
 
    default:
      return (static_cast<double>(n) * LegendrePolynomial(n - 1, x) -
              static_cast<double>(n) * x * LegendrePolynomial(n, x)) /
             (1.0 - x * x);
  }
}
 
double LegendrePolynomialPrimePrime(const int n, const double x) {
  return (2 * x * LegendrePolynomialPrime(n, x) -
          n * (n + 1) * LegendrePolynomial(n, x)) /
         (1.0 - x * x);
}
 
std::vector<double> AllLegendreRoots(const int n) {
  unsigned int n_even =
      static_cast<unsigned int>(std::floor(static_cast<double>(n) / 2.0));
  std::vector<double> computed_roots(n);
 
  for (auto i = 0; i < n_even; i++) {
    computed_roots[n - 1 - i] = LegendreRoot(n, i);
    computed_roots[i] = -computed_roots[n - 1 - i];
  }
  return computed_roots;
}
 
double LegendreRoot(const int n, const int k) {
  double x_old = ApproximateLegendreRoot(n, k);
  unsigned int safety = 10000;
  for (auto i = 0; i < safety; i++) {
    double x_new = x_old - LegendrePolynomial(n, x_old) /
                               LegendrePolynomialPrime(n, x_old);
    double relative_error = std::abs((x_new - x_old) / std::max(1.0, x_new));
    double backward_error = std::abs(LegendrePolynomial(n, x_new));
 
    if (backward_error < 1e-8 && relative_error < TOLERANCE) return x_new;
    x_old = x_new;
  }
  throw std::runtime_error("LegendreRoot did not converge.");
}
 
double ApproximateLegendreRoot(const int n, const int k) {
  if (k > n) throw std::invalid_argument("k must be less than n.");
 
  return std::cos((4.0 * static_cast<double>(k) + 3.0) /
                  (4.0 * static_cast<double>(n) + 2.0) * M_PI);
}
 
std::vector<double> AllLegendrePrimeRoots(const int n) {
  unsigned int n_even =
      static_cast<unsigned int>(std::floor(static_cast<double>(n - 1) / 2.0));
  std::vector<double> computed_roots(n - 1);
 
  for (auto i = 0; i < n_even; i++) {
    computed_roots[i] = LegendrePrimeRoot(n, i);
    computed_roots[n - 2 - i] = -computed_roots[i];
  }
  return computed_roots;
}
 
double LegendrePrimeRoot(const int n, const int k) {
  double x_old = ApproximateLegendrePrimeRoot(n, k);
  unsigned int safety = 100;
 
  double relative_error = std::numeric_limits<double>::infinity();
  double backward_error = std::numeric_limits<double>::infinity();
  for (auto i = 0; i < safety; i++) {
    double x_new = x_old - LegendrePolynomialPrime(n, x_old) /
                               LegendrePolynomialPrimePrime(n, x_old);
    relative_error = std::abs((x_new - x_old) / std::max(1.0, x_new));
    backward_error = std::abs(LegendrePolynomialPrime(n, x_new));
 
    if (backward_error < 1e-8 && relative_error < TOLERANCE) return x_new;
    x_old = x_new;
  }
  std::ostringstream error_message;
  error_message << std::scientific << std::setprecision(6)
                << "LegendrePrimeRoot did not converge. "
                << "Backward Error: " << backward_error
                << ", Relative Error: " << relative_error;
 
  throw std::runtime_error(error_message.str());
}
 
double ApproximateLegendrePrimeRoot(const int n, const int k) {
  if (k >= n - 1) throw std::invalid_argument("k must be less than n-1.");
  return -(ApproximateLegendreRoot(n, k) + ApproximateLegendreRoot(n, k + 1)) /
         2.0;
}
 
}  // namespace hummingbird