#include "UtilitiesP.h"
#include "Simpson.h"

int main(void)
{
    double agm(1/functionsP::agm(1, 2));
    double integral( (2/constantsP::kPi) * Simpson(0, constantsP::kPi*0.5, constantsP::kEps, functionsP::eliptic));
    assert(fabs(agm-integral)<constantsP::kEps);

    double firstVal(constantsP::kPi*0.5);
    integral = Simpson(constantsP::kAlmostZero, 10000, constantsP::kEps, functionsP::Dirichlet);

    std::cout << firstVal << std::endl;
    std::cout << integral << std::endl;

    double secondVal(sqrt(constantsP::kPi)*0.5);
    integral = Simpson(constantsP::kAlmostZero, 12, constantsP::kEps, functionsP::Puasson);

    std::cout << secondVal << std::endl;
    std::cout << integral << std::endl;

    double thirdVal(constantsP::kPi/sin(constantsP::kNumberFromEuler*constantsP::kPi));
    integral = Simpson(constantsP::kAlmostZero, 13, constantsP::kEps, functionsP::Euler);

    std::cout << thirdVal << std::endl;
    std::cout << integral << std::endl;


    return 0;
}

#include <stdexcept>
#include <iostream>

//#define NDEBUG
#include <cassert>

namespace constantsP
{
    const double kPi(3.1415926535897932);
    const double kEps(1.e-12);
    const double kAlmostZero(1.e-7);
    const double kBigNumber(30);
    const double kNumberFromEuler(0.4);
};

namespace functionsP
{
    double agm(const double&, const double&);
    double eliptic(const double&);
    double Dirichlet(const double&);
    double Puasson(const double&);
    double Euler(const double&);
};

#include "UtilitiesP.h"

double functionsP::agm(const double& a, const double& b)
{
    if(a <= 0 || b <= 0)
    {throw std::invalid_argument("Arithmetic-geometric mean is defined only for 0 < a < b.");}

    double aPrev(a), bPrev(b), aCurrent(a), bCurrent(b);
    do
    {
        aPrev = aCurrent;
        bPrev = bCurrent;
        aCurrent = sqrt(aPrev*bPrev);
        bCurrent = (aPrev+bPrev)*0.5;
    }while((aPrev < aCurrent) && (bCurrent < bPrev) && (aPrev < bPrev));

    double res((aPrev <= aCurrent)?aPrev:bPrev);

    return res;
}

double functionsP::eliptic(const double& x)
{
    return 1/sqrt(pow(1*sin(x), 2) + pow(2*cos(x), 2));
}

double functionsP::Dirichlet(const double& x)
{
    return sin(2*x)/x;
}

double functionsP::Puasson(const double& x)
{
    return exp(-(x*x));
}

double functionsP::Euler(const double& x)
{
    return pow(x, constantsP::kNumberFromEuler-1)/(1+x);
}

#include <cmath>

double Simpson (const double& a, const double& b, const double& eps, double(*const f)(const double&));

#include "Simpson.h"

double Simpson (const double& a, const double& b, const double& eps, double(*const f)(const double&))
{
    double step((b-a)*0.5);
    double s1(step*(f(a)+f(b))), s2(0), s4(4*step*f(a+step));
    double previousSum(0), currentSum(s1+s4);
    int n(2);
    do
    {
        previousSum = currentSum;
        n += n;
        step *= 0.5;

        s1 *= 0.5;
        s2 = 0.5*s2+0.25*s4;

        s4 = 0;
        for(int i(1); i < n; i+=2)
        {s4 += f(a+i*step);}
        s4 = 4*step*s4;

        currentSum = s1+s2+s4;
    }
    while(eps<fabs(previousSum-currentSum));

    return currentSum/3;
}