#include <gtest/gtest.h>
#include "modules/OrthPol.cpp"
#include "ExplicitExpressions.h"
#include "SumFormulas.cpp"

class OrthPolTest : public testing::Test {};

// Test that OrthPol returns correct values if x=0
TEST_F(OrthPolTest, ChebPolZeroCase) {
   double zero = 0.0;
   for (int i=0; i<=50; ++i) 
   {

      std::vector<double> res = OrthPol(1, i, 0);
      if (i % 2 == 1)
      {
        EXPECT_DOUBLE_EQ(res[0], zero) << "Chebyshev polinom of the first kind " << i << "th-order are not equal 0";
        EXPECT_DOUBLE_EQ(res[1], pow(-1.0, (i-1)/2)*i) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal "<< i << "*(-1)^" << (i-1)/2;
      }
      else
      {
        EXPECT_DOUBLE_EQ(res[0], pow(-1.0, i/2)) << "Chebyshev polinom of the first kind " << i << "th-order are not equal (-1)^" << i/2;
        EXPECT_DOUBLE_EQ(res[1], zero) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal 0";
      }
   }
   for (int i=0; i<=50; ++i) 
   {
      std::vector<double> res = OrthPol(2, i, 0);
      if (i % 2 == 1)
      {
        EXPECT_DOUBLE_EQ(res[0], zero) << "Chebyshev polinom of the second kind " << i << "th-order are not equal 0";
        EXPECT_DOUBLE_EQ(res[1], pow(-1.0, (i-1)/2)*(i+1)) << "Derivative of Chebyshev polinom of the second kind " << i << "th-order are not equal "<< (i+1) << "*(-1)^" << (i-1)/2;
      }
      else
      {
        EXPECT_DOUBLE_EQ(res[0], pow(-1.0, i/2)) << "Chebyshev polinom of the second kind " << i << "th-order are not equal (-1)^" << i/2;
        EXPECT_DOUBLE_EQ(res[1], zero) << "Derivative of Chebyshev polinom of the second kind " << i << "th-order are not equal 0";
      }
   }
}

// Test that OrthPol returns correct values if x=+-1
TEST_F(OrthPolTest, ChebPolOneCase) {
   for (int i=0; i<=50; ++i) 
   {
    std::vector<double> res = OrthPol(1, i, 1);
    EXPECT_DOUBLE_EQ(res[0], 1) << "Chebyshev polinom of the first kind " << i << "th-order are not equal " << 1;
    EXPECT_DOUBLE_EQ(res[1], pow(i, 2)) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal "<< pow(i, 2);

    res = OrthPol(2, i, 1);
    EXPECT_DOUBLE_EQ(res[0], i+1) << "Chebyshev polinom of the second kind " << i << "th-order are not equal " << i+1;
    EXPECT_DOUBLE_EQ(res[1], i*(i+1)*(i+2)/3) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal "<< i*(i+1)*(i+2)/3;

    res = OrthPol(1, i, -1);
    EXPECT_DOUBLE_EQ(res[0], pow(-1, i)) << "Chebyshev polinom of the first kind " << i << "th-order are not equal " << pow(-1, i);
    EXPECT_DOUBLE_EQ(res[1], pow(-1, i-1) * pow(i, 2)) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal "<< pow(-1, i-1) * pow(i, 2);

    res = OrthPol(2, i, -1);
    EXPECT_DOUBLE_EQ(res[0], pow(-1, i)*(i+1)) << "Chebyshev polinom of the first kind " << i << "th-order are not equal " << pow(-1, i)*(i+1);
    EXPECT_DOUBLE_EQ(res[1], pow(-1, i+1) * (i*(i+1)*(i+2)/3)) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal "<< pow(-1, i+1) * (i*(i+1)*(i+2)/3);
   }
}

TEST_F(OrthPolTest, ExplicitExprChebyshevCase){
    double x = -1; // the left boundary
    double length = 2; //the length of the segmment
    double step = 0.2;
    std::vector<double> res_rec;
    double res_expl;
    for (int i=0; i <= static_cast<int>(length / step); ++i)
    {
        for (int j=0; j <= 20; ++j)
        {
            res_rec = OrthPol(1, j, x);
            res_expl = ExplicitExpressionCheb(1, j, x);
            if (res_rec[0] != 0.0)
            {
                EXPECT_NEAR(1, res_expl / res_rec[0], EPS_CHEB);
            }
            else
            {
                EXPECT_NEAR(res_rec[0], res_expl, EPS_CHEB);
            }

            res_rec = OrthPol(2, j, x);
            res_expl = ExplicitExpressionCheb(2, j, x);
            if (res_rec[0] != 0.0)
            {
                EXPECT_NEAR(1, res_expl / res_rec[0], EPS_CHEB);
            }
            else
            {
                EXPECT_NEAR(res_rec[0], res_expl, EPS_CHEB);
            }

        }
        x += step;
    }
}

TEST_F(OrthPolTest, ExplicitExprLaguerreCase){
    double x = -10; // the left boundary
    double length = 20; //the length of the segmment
    double step = 0.5;
    std::vector<double> res_rec;
    double res_expl;
    for (int i=0; i <= static_cast<int>(length / step); ++i)
    {
        for (int j=0; j <= 20; ++j)
        {
            res_rec = OrthPol(3, j, x);
            res_expl = ExplicitExpressionLaguerre(j, x);
            if (res_rec[0] != 0.0)
            {
                EXPECT_NEAR(1, res_expl / res_rec[0], EPS_LAG);
            }
            else
            {
                EXPECT_NEAR(res_rec[0], res_expl, EPS_LAG);
            }
        }
        x += step;
    }
}

TEST_F(OrthPolTest, ExplicitExprHermiteCase){
    double x = -10; // the left boundary
    double length = 20; //the length of the segmment
    double step = 0.5;
    std::vector<double> res_rec;
    double res_expl;
    for (int i=0; i <= static_cast<int>(length / step); ++i)
    {
        for (int j=0; j <= 20; ++j)
        {
            res_rec = OrthPol(4, j, x);
            res_expl = ExplicitExpressionHermite(j, x);
            if (res_rec[0] != 0.0)
            {
                EXPECT_NEAR(1, res_expl / res_rec[0], EPS_HERM);
            }
            else
            {
                EXPECT_NEAR(res_rec[0], res_expl, EPS_HERM);
            }
        }
        x += step;
    }
}

TEST_F(OrthPolTest, SumFormulaChebCase){
    double x = -1; // the left boundary
    double length = 2; //the length of the segmment
    double step = 0.2;
    std::vector<double> U;
    std::vector<double> T;
    double left_part;
    double right_part;
    for (int i=0; i <= static_cast<int>(length / step); ++i)
    {
        for (int j=0; j <= 200; ++j)
        {
            left_part = SumFormulaCheb1(j, x);
            U = OrthPol(2, 2*j, x);
            right_part = 0.5 * (1 + U[0]);
            if (right_part != 0){
                EXPECT_NEAR(left_part / right_part, 1, pow(10, -9));
            }
            else{
                EXPECT_NEAR(left_part, right_part, pow(10, -9));
            }

            if (2*j > 1){
                left_part = SumFormulaCheb2(j, x);
                U = OrthPol(2, 2*j - 1, x);
                right_part = 0.5 * U[0];
                if (right_part != 0){
                    EXPECT_NEAR(left_part / right_part, 1, pow(10, -9));
                }
                else{
                    EXPECT_NEAR(left_part, right_part, pow(10, -9));
                }
            }

            if (x != 1 && x != -1){
                left_part = SumFormulaCheb3(j, x);
                T = OrthPol(1, 2*j + 2, x);
                right_part = 0.5 * (1 - T[0]) / (1 - pow(x, 2));
                if (right_part != 0){
                    EXPECT_NEAR(left_part / right_part, 1, pow(10, -9));
                }
                else{
                    EXPECT_NEAR(left_part, right_part, pow(10, -9));
                }
            }

            if (x != 1 && x != -1){
                 left_part = SumFormulaCheb4(j, x);
                 T = OrthPol(1, 2*j + 1, x);
                 right_part = 0.5 * (x - T[0]) / (1 - pow(x, 2));
                if (right_part != 0){
                    EXPECT_NEAR(left_part / right_part, 1, pow(10, -9));
                }
                else{
                    EXPECT_NEAR(left_part, right_part, pow(10, -9));
                }
            }

        }
        x += step;
    }
}

TEST_F(OrthPolTest, SumFormulaLagCase){
    double x = -1; // the left boundary
    double length = 2; //the length of the segmment
    double step = 0.1;
    std::vector<double> L;
    double left_part;
    double right_part;
    for (int i=0; i <= static_cast<int>(length / step); ++i)
    {
        for (int j=0; j <= 25; ++j)
        {
            for (double mu=0.1; mu <=0.9; mu += 0.1){
                left_part = SumFormulaLag1(j, mu, x);
                L = OrthPol(3, j, mu * x);
                right_part = L[0];
            if (right_part != 0){
                EXPECT_NEAR(left_part / right_part, 1, pow(10, -9));
            }
            else{
                EXPECT_NEAR(left_part, right_part, pow(10, -9));
            }            }
        }
        x += step;
    }
}

TEST_F(OrthPolTest, SumFormulaHermiteCase){
    double x = -1; // the left boundary
    double length = 2; //the length of the segmment
    double step = 0.1;
    std::vector<double> H;
    double left_part;
    double right_part;
    for (int i=0; i <= static_cast<int>(length / step); ++i)
    {
        for (int j=0; j <= 25; ++j)
        {
            right_part = SumFormulaHermite1(j, x);
            H = OrthPol(4, j, 2 * x);
            left_part = H[0];
            if (left_part != 0){
                EXPECT_NEAR(right_part / left_part, 1, pow(10, -9));
            }
            else{
                EXPECT_NEAR(left_part, right_part, pow(10, -9));
            }        }
        x += step;
    }
}

int main(int argc, char **argv) {
  testing::InitGoogleTest(&argc, argv);
  return RUN_ALL_TESTS();
}