Implementing numerical integration I have this C++ code that implements a rectangular numerical integration.  It evaluates the \$K\$-dimensional integral 

$$\int_{u_K = 0}^{\gamma}\int_{u_{K-1} = 0}^{\gamma-u_K}\cdots\int_{u_2}^{\gamma-u_K-\cdots-u_3}F_U(\gamma-\sum_{k=2}^Ku_k)f_U(u_2)\cdots f_U(u_K)\,du_2\cdots du_K$$
where \$F_U\$ is the cdf and \$f_U\$ is the pdf.

#include &lt;iostream&gt;
#include &lt;cmath&gt;
using namespace std;

float pdf(float u){
     return (1/(pow(1+u, 2)));
}

float cdf(float u){
      return (1 - 1/(u+1));
}
// The main function that implements the numerical integration, 
//and it is a recursive function
float integ(float h, int k, float du){
      float res = 0;
      if (k == 1){
         res =  cdf(h);
      }else{
         float u = 0;
      while (u &lt; h){
         res += integ(h - u, k - 1, du)*pdf(u)*du;
         u += du;
     }
}
     return res;
}

int main(){
    float du = 0.0001;
    int K = 3;
    float gamma[4] = {0.31622777, 0.79432823, 
                1.99526231, 5.01187234};
    int G = 50;
    int Q = 2;
    for (int i = 0; i &lt; 4; i++){
        if ((G-Q*(K-1)) &gt; 0){
            float gammath = (gamma[i]/Q)*(G-Q*(K-1));
            cout&lt;&lt;1-integ(gammath, K, du)&lt;&lt; endl;
    }

 }

    return 0;
}


I am facing a speed problem, although I switched to C++ from Python and MATLAB, because C++ is faster. The problem is that I need a small step size du to get an  accurate evaluation of the integration. 

Basically, I want to evaluate the integral at 4 different points defined by gammath, which is a function of other defined parameters. 

Is there anyway I can speed up this program? I already have 25x+ speed factor over the same code in Python, but still the code takes too long (I ran it all night, and it wasn't finished in the morning). And this is only for K=3, and G=50. In other cases I want to test K = 10, and G = 100 or 300.
I see a few things that may help you improve your program.

Avoid pow for float

The use of pow converts the argument to a double and returns a double.  If you're casting the result to a float anyway, the use of std::powf would be more appropriate.  But see the next suggestion for an even better approach.

Prefer multiplication to pow with small, fixed integral powers

The content of the pdf function could be rewritten like this instead:

 return 1.0f/((u+1)*(u+1));


Usealittlemorewhitespace

With everything crowded together on the line, it's hard to read and understand some lines like this:

float gammath = (gamma[i]/Q)*(G-Q*(K-1));


I'd write it like this instead:

float gammath{gamma[i] / Q) * (G - Q * (K - 1)};


This also uses the C++11 uniform initialization syntax, which I prefer.

Use constexpr where appropriate

Any time you can make a numerical constant or expression constexpr is one more opportunity for compiler optimization.  In this case, virtually everything in main, including du, K, gamma, G, and Q can be made constexpr simply by adding that keyword in front of the declaration.  With that, the compiler is free to evaluate expressions like this:

G - Q * (K - 1)


once at compile-time rather than at runtime.  It can save quite a bit of time if you use it.

Move invariants out of the loop

The value of the expression listed above does not change within the loop (or within the program at all, as presented).  Smart compilers can usually optimize that by themselves, but doing it manually may give you insight.  For example, we can make a constant of that expression outside the loop:

constexpr auto m{G - Q * (K - 1)};


Now the inner if boils down to if (m &gt; 0) but since it's invariant within the loop, move it outside.  Now we have the encapsulating if outside and the for loop inside which makes a lot more sense.

Use C++11 threading features

The calculations for each value of gamma are completely independent.  For that reason, they could easily be calculated in parallel.  Here's a way to do that:

int main() {
    constexpr float du = 0.0001;
    constexpr int K = 3;
    constexpr float gamma[4] = {0.31622777, 0.79432823, 1.99526231, 5.01187234};
    constexpr int G = 50;
    constexpr int Q = 2;
    constexpr auto m{G - Q * (K - 1)};
    std::vector&lt;std::future&lt;float&gt;&gt; chunks;
    if (m &gt; 0) {
        for (int i = 0; i &lt; 4; i++) {
            float gammath{gamma[i] / Q * m};
            chunks.push_back(std::async(integ, gammath, K, du));
        }
    }
    std::mutex cout_mutex;
    for (int i=0; i &lt; 4; ++i) {
        std::lock_guard&lt;std::mutex&gt; lock(cout_mutex);
        std::cout &lt;&lt; gamma[i] &lt;&lt; '\t' &lt;&lt; 1 - chunks[i].get() &lt;&lt; '\n';
    }
}


This requires &lt;future&gt;, &lt;mutex&gt;, and &lt;vector&gt; as #includes.

Parallelize the recursive integrals

If there are 100,000 steps in the integral, why do that with a single core?  If you have a multicore machine like I do, you could divide that work among multiple cores by using the same kind of std::async technique shown above.  This would dramatically speed the results (and can be divided among as many cores as you have).  I did not implement that, but the code above should give a hint as to how to do such a thing generically.

Use static for functions which can have file scope

If the pdf, cdf and integ functions are never used outside the file in which they're defined, the compiler can get quite clever about inlining code and other optimizations.  For that reason, I'd recommend making those functions static if at all possible.

Results

I didn't have the patience to run this at the du value in the program (0.0001) so I ran it at 50 times that at 0.005 to speed up testing.  Obviously, that's an option if you can accept the lower resolution, but I just did it to speed up timing tests.

The original (with the modified value of du) took 1.06 seconds, while the modified version took 0.723s.  At full resolution (du = 0.0001) the results took 30 minutes and 3 seconds on the same 64-bit 8-core Linux machine.