// Functions from this file are implemented using the book "Handbook of Mathematical Function With Formulas, Graphs and Mathematical Tables" 
// by Milton Abramowitz and Irene A. Stegun, tenth printing, 1972.
// Table 22.12 "Sum Formulas"

#pragma once
#include <iostream>
#include <cmath>
#include <modules/OrthPol.h>

double SumFormulaCheb1(const int& n, const double& x) // The left part of 22.12.2
{
    double res = 0;
    for (int m=0; m<=n; ++m)
    {
        std::vector<double> T = OrthPol(1, 2*m, x);
        res += T[0];
    }
    return res;
}

double SumFormulaCheb2(const int& n, const double& x) // The left part of 22.12.3
{
    double res = 0;
    for (int m=0; m<=n-1; ++m)
    {
        std::vector<double> T = OrthPol(1, 2*m + 1, x);
        res += T[0];
    }
    return res;
}

double SumFormulaCheb3(const int& n, const double& x) // The left part of 22.12.4
{
    double res = 0;
    for (int m=0; m<=n; ++m)
    {
        std::vector<double> U = OrthPol(2, 2*m, x);
        res += U[0];
    }
    return res;
}

double SumFormulaCheb4(const int& n, const double& x) // The left part of 22.12.5
{
    double res = 0;
    for (int m=0; m<=n-1; ++m)
    {
        std::vector<double> U = OrthPol(2, 2*m + 1, x);
        res += U[0];
    }
    return res;
}

int C_n_k(const int& n, int k)
{
    if (k > n) return 0;
    if (k * 2 > n) k = n-k;
    if (k == 0) return 1;

    int result = n;
    for( int i = 2; i <= k; ++i ) {
        result *= (n-i+1);
        result /= i;
    }
    return result;
}

double SumFormulaLag1(const int& n, const double& mu, const double& x) // The left part of 22.12.7
{
    double res = 0;
    for (int m=0; m<=n; ++m)
    {
        std::vector<double> L = OrthPol(3, n - m, x);
        int C_nm = C_n_k(n, m);
        res += C_nm * pow(mu, n - m) * pow(1 - mu, m) * L[0];
    }
    return res;
}

double SumFormulaHermite1(const int& n, const double& x) // The right part of 22.12.8, case x = y
{
    double res = 0;
    for (int k=0; k<=n; ++k)
    {   
        std::vector<double> H_1 = OrthPol(4, k, pow(2, 0.5) * x);
        std::vector<double> H_2 = OrthPol(4, n - k, pow(2, 0.5) * x);
        int C_nk = C_n_k(n, k);
        res += C_nk * H_1[0] * H_2[0];
    }
    double m = static_cast<double>(n) / 2;
    res = res / pow(2, m);
    return res;
}