# A probability distribution function, to be called repeatedly during numerical integration

I am trying to speed up as much a possible this function in C++. As I explained in another post, Implementing multidimensional integral for a custom function in C++, this function will be used inside multi-dimensional integrals.

Do you have any tips that could help me speeding up the computation? For instance, by improving the for-loop?

double PDFfunction(double invL, int t, double invtau, double x0, double x, int n_lim) {

const double c = (M_PI/2) * (M_PI/2) * (2 * t * invtau);

double res = 0;

for(int n = 1; n <= n_lim; ++n){

res += exp(-1 * (n * n) * c) * cos(n * M_PI * x * invL) * cos(n * M_PI * x0 * invL);

}

return invL + (2 * invL * res);
}

• What are some typical values? About how high is n_lim, dozens or thousands or millions? BTW I guess the first parameter is invL or L but either way something is wrong there. May 23, 2022 at 23:53
• Thanks, I added the corrections to the question. I used the wrong script. Anyway, n_lim is on the order of hundreds (usually few tens). May 23, 2022 at 23:58

This is a possible solution to harolds open point about exp in his answer.

Note that exp(x+y) = exp(x) * exp(y) holds for any real numbers x,y.

With this, your exponential term exp(-c * n^2) can be rewritten into exp(-c * (n-1)^2) * exp(c) * exp(-2c * n). Note that the first term with (n-1) in it is the exponential function evaluated at the last step, and the second term is just a constant which can be evaluated once.
This leaves us with the last term: exp(-2c * n) = exp(-2c * (n-1)) * exp(-2c). And again, the first term is known from the last step, and the second term is just a constant.

So for the exponential term you need exactly one evaluation of exp(c). Everything else can be done with saving the two mentioned terms from the last iteration and a few multiplications.

However I have no idea what will happen with the accuracy. But given that only multiplications are used I think it should not be too problematic.
To counteract possible inaccuracies you could think about doing something like calculating the exponential term "accurately" using std::exp every x iterations.

Regarding loop unrolling which was mentioned in harolds answer: The values now actually depend on each other. But the above can of course be rewritten so that the new value depends on the value from four iterations earlier instead on the value from the last iteration. Maybe the compiler is even smart enough to figure that out, but I wouldn't count on it and at least check if doing it manually makes it go faster.

• Excellent, very fast May 24, 2022 at 22:35

A typical way to reduce such cosines with an angle that is steadily counting up by the same increment, is to take a vector and rotate it step by step. That way, one cosine and one sine are calculated, but the main calculation involves just a couple of multiplications and additions:

double PDFfunctionFaster(double invL, int t, double invtau, double x0, float x, int n_lim) {

const double cc = (M_PI / 2) * (M_PI / 2) * (2 * t * invtau);

double res = 0;

double cT1 = cos(M_PI * x * invL);
double sT1 = sin(M_PI * x * invL);
double cT2 = cos(M_PI * x0 * invL);
double sT2 = sin(M_PI * x0 * invL);
double a = 1;
double b = 0;
double c = 1;
double d = 0;

for (int n = 1; n <= n_lim; ++n) {
double t;
// rotate vector (a, b)
t = cT1 * a - sT1 * b;
b = sT1 * a + cT1 * b;
a = t;
// rotate vector (c, d)
t = cT2 * c - sT2 * d;
d = sT2 * c + cT2 * d;
c = t;
// the first coordinate of each vector is the cosine of its angle
res += exp(-1 * (n * n) * cc) * a * c;
}

return invL + (2 * invL * res);
}


This has a mild effect on the accuracy, but with n_lim on the order of hundreds or dozens, it's still really good.

On my PC this had a significant effect, shaving off about 30% of the time. It would have been more if it hadn't been for the exp, which is still very expensive (approximately 90% of the time is spent there), but I did not see a way to eliminate it.

While this transformation has introduced loop-carried dependencies, the exp is slow enough that they're not the problem, and you can still employ parallelization and SIMD because these dependencies are just an artifact of how the computation was arranged, not an inherent part of the calculation. For example to split the calculation into N independent parts, you could calculate N successive vectors first (by applying rotations), and have each task apply a rotation by N times the original angle (the corresponding sine/cosine pair is easily calculated up-front), and perform n += N in its loop.

If the exp can be rewritten as well (into what?) then it may pay to unroll the loop by a small factor (by reusing the same trick: calculate more initial vectors, and then rotate by a larger amount) to work around the loop-carried dependencies.

If the function is repeatedly called with the same t and invtau then the exponential part could be precomputed once and reused for each of those calls. I couldn't tell from your other post whether that is the case, but if it is, that could save a lot of time, since that's now the most expensive part by a large margin.

Using BlameTheBits' trick, which if I did it right could look something like this:

double PDFfunctionFaster(double invL, int t, double invtau, double x0, float x, int n_lim) {

const double cc = (M_PI / 2) * (M_PI / 2) * (2 * t * invtau);

double res = 0;

double cT1 = cos(M_PI * x * invL);
double sT1 = sin(M_PI * x * invL);
double cT2 = cos(M_PI * x0 * invL);
double sT2 = sin(M_PI * x0 * invL);
double a = 1;
double b = 0;
double c = 1;
double d = 0;
double ec = exp(cc);
double e = 1.0 / ec; // exp(-c * n^2)
double e2c = e * e;
double e2n = e2c; // exp(-2c * n)

for (int n = 1; n <= n_lim; ++n) {
double t;
// rotate vector (a, b)
t = cT1 * a - sT1 * b;
b = sT1 * a + cT1 * b;
a = t;
// rotate vector (c, d)
t = cT2 * c - sT2 * d;
d = sT2 * c + cT2 * d;
c = t;
// the first coordinate of each vector is the cosine
res += e * a * c;
// update exponential term
e2n *= e2c;
e *= ec * e2n;
}

return invL + (2 * invL * res);
}


In my tests this gave good results (not significantly inaccurate), and it's quite fast now, more than ten times as fast as the original on my PC. You could still "split" the loop-carried dependencies and add SIMD for extra speed.

• The PDFfunction is used to integrate out the variable t in a range between 0 and Infinity, so t is definitely not constant. However, invtau remains constant, so indeed I could give as argument of the PDFfunction a variable called combo, for instance, which I calculate a priori as combo = (M_PI / 2) * (M_PI / 2) * (2 * invtau). May 24, 2022 at 3:07
• I tried thePDFfunctionFaster inside a 4-dimensional integral and it does not really improve code. They basically have the same time run. May 24, 2022 at 19:34
• Sorry for my ignorance, but when you say SMID you mean parallelisation, like OpenMP? May 24, 2022 at 21:44
• I may have found a solution to the exp problem. See my answer below. (It got too complex and too hard to describe it understandably in a small comment so I wrote some kind of "extension answer".) May 24, 2022 at 22:10
• @CafféSospeso parallelism, but without threading. You can still add threading on top. Anyway there's an even more useful trick now thanks to BlameTheBits May 24, 2022 at 22:37

The best way to improve a for loop is not to use one at all.

Let’s start at the top:

const double c = (M_PI/2) * (M_PI/2) * (2 * t * invtau);


M_PI is not portable; it’s not a standard constant, despite what most people think. Normally the correct thing to do is define it yourself, usually with a guard in case it’s predefined. Starting in C++20, there is a new standard set of predefined constants:

// Need this:
#include <numbers>

const double c = std::pow(std::numbers::pi / 2, 2) * (2 * t * invtau);


Now the meat of the function is, of course, this loop, and that’s what you really want to optimize:

    double res = 0;

for(int n = 1; n <= n_lim; ++n){

res += exp(-1 * (n * n) * c) * cos(n * M_PI * x * invL) * cos(n * M_PI * x0 * invL);

}


As a general rule, you should not write naked for loops in modern C++. You should look at the loop, and consider what it’s really doing, and then use an algorithm—usually a standard algorithm—that does the job. In this case, this is clearly a reduction operation. That gives us a couple of options, but the safe option to start with is std::accumulate().

    auto const indices = std::ranges::iota_view{1, n_lim + 1};
auto const res = std::accumulate(
indices.begin(), indices.end(),
0.0,
[invL, t, x, x0, c](auto res, auto n)
{
auto const coeff = n * std::numbers::pi * invL;

return res + std::exp(-1 * std::pow(n, 2) * c) * std::cos(coeff * x) * std::cos(coeff * x0);
});


When I insert that into the function and benchmark it on my machine, it’s 1.1× faster than the original… but that’s probably bullshit; it’s so close, it’s probably just noise. In reality, there’s probably literally no difference between this and the original.

Okay, but can we improve this?

Well, the first thing I notice is that there is no reason the loop has to iterate in order. In other words, I could theoretically take your original loop and do this:

auto const n_mid = n_lim / 2;   // assuming n_lim >= 2

// These two loops can be executed concurrently, on different threads:
auto res1 = 0.0;
for (auto n = 1; n < n_mid; ++n)
{
res1 += exp(-1 * (n * n) * c) * cos(n * M_PI * x * invL) * cos(n * M_PI * x0 * invL);
}

auto res2 = 0.0;
for (auto n = n_mid; n <= n_lim; ++n)
{
res2 += exp(-1 * (n * n) * c) * cos(n * M_PI * x * invL) * cos(n * M_PI * x0 * invL);
}

auto const res = res1 + res2;


Indeed, in theory, every single iteration can be done concurrently, and in any order.

Now std::accumulate() does a reduction, but strictly in order. std::reduce() also does a reduction… but out of order. That means it can theoretically be faster.

So, in theory, we can get a speedup simply be replacing std::accumulate() with std::reduce():

    auto const indices = std::ranges::iota_view{1, n_lim + 1};
auto const res = std::reduce(
indices.begin(), indices.end(),
0.0,
[invL, t, x, x0, c](auto res, auto n)
{
auto const coeff = n * std::numbers::pi * invL;

return res + std::exp(-1 * std::pow(n, 2) * c) * std::cos(coeff * x) * std::cos(coeff * x0);
});


And we can go even a step further. We can explicitly tell the compiler that not only can each iteration be done in any order, they can be done concurrently. To that, we simply add execution policies, either std::execution::par_unseq or (starting in C++20) std::execution::unseq.

Now, realistically, right now, in mid-2022, no compiler is sophisticated enough to really take advantage of std::execution::par_unseq or std::execution::unseq… or even the out-of-order benefits of std::reduce() over std::accumulate(). So, I wouldn’t expect any improvement today. (And, in fact, in my tests, the final version is less than 0.02% faster than the std::accumulate() version, which, again, is just noise.)

So—again, realistically—on any contemporary compiler, there’s going to be no real difference between your original code, and my final code:

double PDFfunction(double invL, int t, double invtau, double x0, double x, int n_lim)
{
auto const c = std::pow(std::numbers::pi / 2, 2) * (2 * t * invtau);

auto const indices = std::ranges::iota_view{1, n_lim + 1};
auto const res = std::reduce(
std::execution::par_unseq,
indices.begin(), indices.end(),
0.0,
[invL, t, x, x0, c](auto res, auto n)
{
auto const coeff = n * std::numbers::pi * invL;

return res + std::exp(-1 * std::pow(n, 2) * c) * std::cos(coeff * x) * std::cos(coeff * x0);
});

return invL + (2 * invL * res);
}


And, in fact, Quick Bench shows no real difference. (Note that I had to switch the compiler from Clang to GCC, because Clang has a bug that breaks std::ranges::iota_view.)

However….

Using the algorithm creates the potential for future compiler advances to optimize the code. With the naked for loop, only a really, really smart compiler can maybe detect that the loop body can be parallelized and vectorized… but even then it probably shouldn’t, and it shouldn’t assume it can do the loop out of order, because that might change the behaviour (not sure about this, would have to think about it). But with the algorithm, it’s explicit; you’re not relying on the good graces of the compiler or its willingness to do things that might surprise you, you’re literally telling it: this loop can be done out-of-order, in parallel, concurrently, and even vectorized. If the compiler is sophisticated enough to understand that, and to actually generate code to do that, then you’ll see performance gains.

You might even want to consider offering even more control, and allowing users to choose the execution policy:

template <typename ExecutionPolicy>
requires std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>
auto PDFfunction(ExecutionPolicy&& policy, double invL, int t, double invtau, double x0, double x, int n_lim)
{
auto const c = std::pow(std::numbers::pi / 2, 2) * (2 * t * invtau);

auto const indices = std::ranges::iota_view{1, n_lim + 1};
auto const res = std::reduce(
std::forward<ExecutionPolicy>(policy),
indices.begin(), indices.end(),
0.0,
[invL, t, x, x0, c](auto res, auto n)
{
auto const coeff = n * std::numbers::pi * invL;

return res + std::exp(-1 * std::pow(n, 2) * c) * std::cos(coeff * x) * std::cos(coeff * x0);
});

return invL + (2 * invL * res);
}

inline auto PDFfunction(double invL, int t, double invtau, double x0, double x, int n_lim)
{
return PDFfunction(std::execution::unseq, invL, t, invtau, x0, x, n_lim);
}


Today, there’s really no way to make your own execution policies, but in the future, you will be able to, for example, create a thread pool, and then pass a policy that will tell std::reduce() that it can split its work across the thread pool. Depending on the number of threads, that could be a speed improvement of several times.

Other future benefits would include being able to mark everything constexpr, and a range version of std::reduce(), so you can create the iota_view right in place, instead of needing to do it separately.

Anywho, the bottom line is this: the most dramatic way you can speed up your function is if you can get that loop to be done concurrently or vectorized. There are non-standard (non C++ standard, that is) ways to do that today, like OpenMP. And there are standard ways to do it—such as using the right algorithms (std::reduce()) and using execution policies—that are technically available today… but… while it possible to request concurrency and/or vectorization, I don’t think compilers are quite sophisticated enough to really use these requests. So, in practice, you won’t see massive speedups today… but you will, in the future.

• So, do you think that using a OpenMP could be a good alternative to speed up the code? (As explained for instance here: stackoverflow.com/questions/54525473/…). Also, I was thinking that perhaps you may not see time differences in these examples, while you may see them when increasing the number of calculations, for instance when I integrate PDFfunction over time between 0 and Infinity. May 24, 2022 at 17:41
• I don’t know OpenMP well enough to answer; I’ve never really used it because it’s so inflexible. What you need to know is: if the OpenMP pragma creates N “things” and spreads the loop over those “things”, then… what, exactly, are the “things”? Kernel threads? Userspace threads? Fibers? Coroutines? GPU threads? When exactly are they created? At the start of the program or (every time the function is called) at the start of the loop? You’d have to profile, but unfortunately, Quick Bench doesn’t support OpenMP, so you’ll have to check manually.
– indi
May 25, 2022 at 20:32
• I should make clear that OpenMP isn’t magical; there is nothing there you can’t do with current standard C++, and much more efficiently. And with libraries like SYCL, you could, today, write your function to, say, execute massively parallel on the GPU (which would be the way I’d go if I really wanted this function done fast). OpenMP just makes it trivially easy… but that comes at a cost. What makes the future executor standards so cool is that it would be just as easy as OpenMP… but much, much more powerful and flexible.
– indi
May 25, 2022 at 20:33

One final thing to consider is the eventual precision of your answer: you probably don't need to consider terms beyond $$\e^{-36}\$$, which is around machine epsilon for double precision.

(I'd always be tempted to precompute what doesn't change in a loop/reduce/whatever.

static auto const _PI_SQUARE_2 = std::pow(std::numbers::pi, 2) / -2;
⋮
auto const
_c = _PI_SQUARE_2 * t * invtau,
PI_L = std::numbers::pi * invL,
x_PI_L = x * PI_L,
x0_PI_L = x0 * PI_L;
⋮
exp(n * n * _c) * cos(n * x_PI_L) * cos(n * x0_PI_L);


- not really expecting speedup, but easier reading.
(oh well, reading PI_L as π/L (π*invL) doesn't carry over to x_PI_L - can't find xPI_L (or x_L_PI?!) much of an improvement there.)
(I wouldn't go quite as far as substituting multiplication of the x_PI_Ls by n by accumulation for fear of avoidably accumulating numerical inaccuracy.)

)
• Just to be sure I understood. What you are basically suggesting is to set as many variable as possible as constants (if possible). Right? May 24, 2022 at 17:21
• I suggest to look for things that don't change between iterations, to at least try to handle them outside the loop. If that happens to be a sub-expression used in the loop, try to assign an illuminating name to it. If that fails, or the use of such a "constant" (it does not really need language support here) impairs readability, turn into a comment about its invariance a revert use in the loop. Code motion has been around for ages and compilers are quite good at it. May 24, 2022 at 19:00
• Ok, it makes sense. May 24, 2022 at 19:02

There are three problems to solve:

• faster cosines
• faster exponentiation
• improved precision

Cosines

In a vein similar to @Harold's answer, I would attempt to use Chebyshev Polynomials of the First Kind. They work well for cos(N*x) calculations. Using the recurrence relation for the Chebyshev polynomials in your loop will require (after removing 2x from the loop) only two multiplications and two subtractions, instead of 8 multiplications, 2 additions and 2 subtractions for the cosine rotations.

For the math part, see https://www.mathsjournal.com/pdf/2016/vol1issue1/PartA/1-1-18-475.pdf

Applying the double-angle formula, you get:

  cos 0y = 1
cos 1y = cos y
cos 2y = 2cos²y - 1
cos 3y = 4cos³y - 3 cos y
etc.


Replacing cos y as x and using the convention Tₘ(cos y) = cos(m*y), you get:

  T₀(x) = 1
T₁(x) = x
T₂(x) = 2x² - 1
T₃(x) = 4x³ - 3x
T₄(x) = 8x⁴ - 8x² + 1



Then the recurrence relation is:

  Tₘ₊₁(x) = 2xTₘ(x) - Tₘ₋₁


You can use this recurrence relation to compute cos((n+1)*x) from cos(n*x). Since 2x is repeated, pull it out of the loop.

Exponentiation

You can use a different recurrence relation here (no Chebyshev needed!)

Outside the loop, compute:

 e2c = exp(-2*c)
e1c = exp(-c)


We will maintain two running products, f and g. Here are the initial values:

  f(1) = exp(-c) = e1c
g(1) = exp(-2*c) = e2c


The recurrence relation is simple:

  f(n+1) = f(n) * g(n) * exp(-c) = f(n) * g(n) * e1c
g(n+1) = g(n) * exp(-2*c) = g(n) * e2c


Your product to use in multiplication in your formula is f(n). As you can see, there are only three multiplications per loop iteration. This should be much faster than the exponentiation call.

Precision

To deal with loss of precision, I would also try employing Kahan summation (also called compensated summation) in your accumulator. This will double your apparent accuracy if you have cancellation problems.

Here is Rosetta code example for C:

https://rosettacode.org/wiki/Kahan_summation#C

This summation has a small cost, but is worth it.

• Both of these strength reductions are inherently sequential. It is likely that a modern vectorizing compiler, such as ICX with AVX2 enabled, or an OpenMP parallel loop, could generate faster code on current hardware. But, clever! May 27, 2022 at 4:02

After some tests with the Compiler Explorer, my first suggestion is to compile this code with ICX (if you are running it on an Intel or AMD CPU). With the flags -std=c++20 -O3 -mavx512f -mavx512cd, either Clang 14.0.0 or ICX 2022 are capable of transforming your code into vectorized code using the 512-bit zmm registers, but GCC 12.1 does not. With the -std=c++20 -march=x86-64-v3 -O3 flags, ICX can compile the code to use AVX, but GCC and Clang do not.

One alternative approach that might improve the code on some compilers is to transform the loop into explicit vector operations using std::valarray. (Warning: untested.)

#include <numbers> // For pi_v
#include <numeric> // For iota
#include <stdexcept>
#include <valarray>

constexpr double pi = std::numbers::pi_v<double>;

static std::valarray<double> iota_valarray( size_t i0, size_t n )
{
if (n < i0) {
throw std::domain_error("Upper bound greater than lower bound.");
}

std::valarray<double> is(n-i0+1);
std::iota( std::begin(is), std::end(is), i0 );
return is;
}

double PDFfunction( const double invL,
const int t,
const double invtau,
const double x0,
const double x,
const int n_lim)
{
if (n_lim < 1) {
throw std::domain_error("Upper bound less than 1.");
}

const double c = (pi/2) * (pi/2) * (2 * t * invtau);

const std::valarray<double> is = iota_valarray( 1, (unsigned)n_lim );
const double res = (std::exp(-c*(is*is)) *
std::cos((pi*x*invL)*is) *
std::cos((pi*x0*invL)*is)
).sum();

return invL + (2 * invL * res);
}


Specifically, I found that this enables 256-bit vectorization on some targets with Clang 14 that the original code did not. It also enables the code to be written as static single assignments, and you might consider it clearer to read. In benchmarks with GCC 11.2, your code, this code, and indi’s all perform similarly.

Since the SIMD operations need the index, this code first creates an iota vector holding the indices. In practice, GCC, Clang and ICX all appear to compile this code using expression templates, avoiding the need to create any temporary valarray objects in memory.

It would also be possible to pass a constexpr function to std::valarray<T>::apply().