# Computing the improper integral of several functions

How can I optimize and make prettier the following code?

Main.cpp

#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;
}


Utilities.h

#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&);
};


Utilities.cpp

#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);
}


Simpson.h

#include <cmath>

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


Simpson.cpp

#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;
}


After running the code you will see that it works well for the first three integrals. But it works very bad for the last one. It takes very much time to compute it.

What can I change in order to make it work faster?

• So I gave this a try, and this is clearly not just a question of optimizing the code. You have something severely busted in there. In the 4th example, n baloons to tens of millions before we get even close to epsilon.
– user128454
Nov 28, 2017 at 4:23

You are assuming the functions in <cmath> (e.g. fabs(), sin()) etc. are in the global namespace already. However, it is unspecified in the standard whether or not these functions exist in the global namespace. The easiest way to ensure they are in the global namespace is to add using namespace std; to the appropriate source file, but that's bad practice. A better practice would be to add using declarations (preferably in the most local scope possible), or simply add std:: to the beginning of each function call.

There is little point in passing a primitive like double by const reference rather than by value. The performance is about the same and it saves you typing to pass by value.

Areas of your code are difficult to read due to lack of spacing. For example:

while(eps<fabs(previousSum-currentSum));


The same line is easier to read like this:

while (eps < fabs(previousSum - currentSum));


Similarly, this:

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


would be easier to read with the { and } on their own lines like this:

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


If your compiler supports C++11 or later you may want to initialize with curly braces ({, }) to avoid the most vexing parse.