C++高精度计算归一化勒让德多项式

#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>
#include<cmath>
using namespace boost::multiprecision;

// 计算阶乘
cpp_int factorial(int n) {
    cpp_int result = 1;

    for (int i = 1; i <= n; i++) {
        result *= i;
    }

    return result;
}

// 计算双阶乘
cpp_int factorial2(int n) {
    cpp_int result = 1;

    for (int i = n; i >= 1; i -= 2) {
        result *= i;
    }

    return result;
}

// 计算勒让德多项式
double legendre(int n, double x) {
    double p_prev = 1.0;
    double p_curr = x;
    double p_next;

    for (int i = 2; i <= n; i++) {
        p_next = ((2 * i - 1) * x * p_curr - (i - 1) * p_prev) / i;
        p_prev = p_curr;
        p_curr = p_next;
    }

    return p_curr;
}

// 计算归一化勒让德多项式
double normalized_legendre(int n, int m, double x) {
    double p_n = legendre(n, x);
    double p_m = legendre(m, x);
    double norm_const = sqrt((2.0 * n + 1) / (2.0 * acos(-1)) * static_cast<double>(factorial(n - m)) / static_cast<double>(factorial(n + m)));

    double factorial2_2n_1 = static_cast<double>(factorial2(2 * n - 1));
    double factorial2_2m_1 = static_cast<double>(factorial2(2 * m - 1));

    return norm_const * p_n / sqrt(1 - x * x) * std::pow(-1, m) * factorial2_2n_1 / factorial2_2m_1;
}

int main() {
    int n = 1000;
    int m = 0;
    double x = 0.5;
    double result = normalized_legendre(n, m, x);

    std::cout << 'P_' << n << '^' << m << '(' << x << ') = ' << result << std::endl;

    return 0;
}