C 代码的半随机缓慢

>我有一个更长的代码，在保持问题的同时尽可能减少。我的代码针对不同的参数值运行 MCMC 计算。对于某些值组合，代码运行时间要长得多，比典型情况慢约 100 倍。但是，它不应该，因为操作数不依赖于参数值。

我正在AMD64 Linux盒子上运行它，其中glibc-2.17在Gentoo上用GCC 4.8.1编译。编译标志在出现时并不重要。我还在另一个带有较旧AMD64处理器的Gentoo盒子上测试了它，结果是一样的。

我做了一堆测试：

• 我尝试使用 Valgrind 进行调试，它没有发现内存问题或其他令人讨厌的事情。
• 其次，我尝试在修复有问题的参数值的情况下运行代码，但没有遇到缓慢的问题。
• 我尝试在迭代之间放置sleep(4)，但没有任何变化。

该问题显示当迭代命中k = 1、i = j = 0时，转化为mu[0] = -0.05、mu[1] = -0.05和mu[2] = 0.05。正如我所说，对所有迭代使用此固定值消除了我所看到的问题。

以下是消除问题的方法：

• 更改限制。
• 固定mu[]系数。
• 从计算中删除dW3。
• 删除rand() .
• 删除q的计算。
• 删除s[j] 的更新。

我读过一些关于slowpow的内容，因此试图通过编写我自己的版本来消除exp。这解决了我在这个MWE上遇到的问题，但当重新实现的exp放在生产代码中时，就解决了

问题：是什么导致了半随机缓慢？

以下是 MWE 的代码。有关如何进行的所有帮助和建议将不胜感激。

注意：此代码是用g++编译的，尽管它本质上是C的。更改编译器不会更改任何内容。

关于分支预测：使用 删除其中一个if语句

q = exp(dW);        
q = q / (1.0 + q); 

无论dW的价值是多少，都不会改变代码的行为;如果这确实是由于分支预测，那一定是由于第二个if。

#include <cstdio>
#include <cstdlib>
#include <cmath>
inline int index(int const i, int const j, int const n)
{
    return (i + n) % n + ((j + n) % n) * n;
}
void get_sample(int* s, int n, double* mu)
{
    for (int i = 0; i < 10 * n * n; i++)
    {
        int j = i % (n * n); 
        int x = j % n;
        int y = (j - x) / n;
        double dW1 = mu[0] * (s[index(x - 1, y, n)] + s[index(x + 1, y, n)] + s[index(x, y - 1, n)] + s[index(x, y + 1, n)]);
        double dW2 = mu[1] * (s[index(x - 1, y - 1, n)] + s[index(x + 1, y - 1, n)] + s[index(x + 1, y + 1, n)] + s[index(x - 1, y + 1, n)]);
        double dW3 = mu[2] * (s[index(x - 1, y, n)] * s[index(x - 1, y - 1, n)] * s[index(x, y - 1, n)] + s[index(x - 1, y, n)] * s[index(x - 1, y + 1, n)] * s[index(x, y + 1, n)]
                                        + s[index(x, y + 1, n)] * s[index(x + 1, y + 1, n)] * s[index(x + 1, y, n)] + s[index(x + 1, y, n)] * s[index(x + 1, y - 1, n)] * s[index(x, y - 1, n)]);
        double dW = 2.0 * (dW1 + dW2 + dW3);
        double q;
        if (dW < 0.0)
        {   
            q = exp(dW);
            q = q / (1.0 + q); 
        }
        else
        {
            q = exp(-dW);
            q = 1.0 / (1.0 + q);
        } 
        double p = ((double) rand()) / ((double) RAND_MAX);
        if (p < q)
        {
            s[j] = 1;
        }
        else
        {
            s[j] = -1;
        }
    }
}
int main(int argc, char** argv)
{
    double mu[3];
    double limits[6] = {-0.05, 0.8, -0.05, 0.45, -0.45, 0.05};
    int s[16];
    for (int i = 0; i < 16; i++)
    {
        s[i] = -1;
    }
    for (int k = 0; k < 2; k++)
    {
        for (int j = 0; j < 2; j++)
        {
            for (int i = 0; i < 2; i++)
            {               
                mu[0] = limits[0] + ((limits[1] - limits[0]) * i);
                mu[1] = limits[2] + ((limits[3] - limits[2]) * j);
                mu[2] = limits[4] + ((limits[5] - limits[4]) * k);
                printf(" Computing (% .6lf, % .6lf, % .6lf)...n", mu[0], mu[1], mu[2]);
                for (int sample = 0; sample < 1000; sample++)
                {
                    get_sample(s, 4, mu);
                }                           
            }
        }
    }               
    return 0;
}