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.
n_lim, dozens or thousands or millions? BTW I guess the first parameter isinvLorLbut either way something is wrong there. \$\endgroup\$