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 two integral. 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?

What comments about the style do you have?