3

我一直在尝试在 OpenCL 中设计一个快速的二进制求幂实现。我当前的实现与本书中关于 pi的实现非常相似。

// Returns 16^n mod ak
inline double expm (long n, double ak)
{
    double r = 16.0;
    long nt;

    if (ak == 1) return 0.;
    if (n == 0) return 1;
    if (n == 1) return fmod(16.0, ak);

    for (nt=1; nt <= n; nt <<=1);

    nt >>= 2;

    do
    {
        r = fmod(r*r, ak);
        if ((n & nt) != 0)
            r = fmod(16.0*r, ak);
        nt >>= 1;
    } while (nt != 0);
    return r;
}

有没有改进的余地？现在我的程序大部分时间都花在这个函数上。

4

2 回答 2

2

我的第一个想法是将其矢量化，以实现约 1.6 倍的潜在加速。与原始中的 2 次乘法相比，每个循环使用 5 次乘法，但对于足够大的 N，循环数大约为四分之一。将所有doubles 转换为s，并将slong换成s 可能会提供一些加速，具体取决于使用的确切 GPU 等等。fmod%

inline double expm(long n, double ak) {

    double4 r = (1.0, 1.0, 1.0, 1.0);
    long4 ns = n & (0x1111111111111111, 0x2222222222222222, 0x4444444444444444,
            0x8888888888888888);
    long nt;

    if(ak == 1) return 0.;

    for(nt=15; nt<n; nt<<=4); //This can probably be vectorized somehow as well.

    do {
        double4 tmp = r*r;
        tmp = tmp*tmp;
        tmp = tmp*tmp;
        r = fmod(tmp*tmp, ak); //Raise it to the 16th power, 
                                       //same as multiplying the exponent 
                                       //(of the result) by 16, same as
                                       //bitshifting the exponent to the right 4 bits.

        r = select(fmod(r*(16.0,256.0,65536.0, 4294967296.0), ak), r, (ns & nt) - 1);
        nt >>= 4;
    } while(nt != 0); //Process n four bits at a time.

    return fmod(r.x*r.y*r.z*r.w, ak); //And then combine all of them.
}

编辑：我很确定它现在可以工作。

于 2014-01-21T13:03:37.870 回答
0
  • 要提取的循环nt = log2(n);可以
    if (n & 1) ...; n >>= 1;
    在 do-while 循环中替换为。
  • 鉴于最初 r = 16;， fmod(r*r, ak) 与 fmod(16*r,ak) 可以很容易地延迟以仅在每 N 次迭代左右计算模数 - 循环展开？
  • 还有为什么要使用 fmod？
于 2014-01-21T05:10:06.623 回答