OpenCL中的快速实现二进制求幂实现答案

【问题标题】：Fast implementation binary exponentiation implementation in OpenCLOpenCL中的快速实现二进制求幂实现
【发布时间】：2014-01-20 23:24:59
【问题描述】：

我一直在尝试在 OpenCL 中设计一个快速的二进制求幂实现。我当前的实现与this book about 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;
}

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

【问题讨论】：

知道输入的一般范围是多少吗？
n 可以达到一百万左右。所以，输入范围相当大。
ak 实际上是双精度数还是整数？它的范围是多少？
抱歉耽搁了 - 我现在回来了。 ak 是（据我所知）一个整数。出于某种原因，此实现将其视为双精度（也许是为了避免在某些时候进行强制转换？）。我实际上并没有用 ak 作为 int 而不是 double 来测试整个程序，所以我不知道进行这样的切换是否会降低其准确性或影响后面部分的性能。

标签： c opencl pow exponentiation modular-arithmetic

【解决方案1】：

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

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.
}

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

【讨论】：

这看起来棒极了！不幸的是，我今天没有太多时间来弄清楚你的代码在做什么。我明天可能会去看看。
谢谢！我想我现在可以工作了。它的要点是它使用向量类型一次处理输入的 4 位。
我今天会测试一下。谢谢！
嗯...它似乎无法正常工作。例如，在第一个工作项的第一个迭代中，它必须计算 expm(6, 9)（换句话说，16^6 % 9）。 expm() 在答案为 1 时返回 9。
看起来我的 select 设置是向后的，但在 return 语句中忘记了 fmod，尽管我不知道这两者是如何导致该错误的。是否还有其他要报告的破案？

【解决方案2】：

提取nt = log2(n);的循环可以替换为
if (n & 1) ...; n >>= 1;
在 do-while 循环中。
考虑到最初 r = 16;，fmod(r*r, ak) 与 fmod(16*r,ak) 可以很容易地延迟计算模数，只在每第 N 次迭代左右 -循环展开？
为什么还要使用 fmod？

【讨论】：