c - 优化 maC，性能)——跟进 bit-twiddling 问题

Question

感谢 Bit twiddling: which bit is set? 上的一些非常有帮助的 stackOverflow 用户，我已经构建了我的函数(发布在问题的末尾)。

任何建议——即使是很小的建议——将不胜感激。希望它能让我的代码变得更好，但至少它应该教会我一些东西。 :)

概述

此函数将被调用至少 1013 次，可能多达 1015 次。也就是说，此代码很可能会运行数月，因此任何性能提示都会有所帮助。

这个函数占程序时间的 72-77%，基于分析和大约十几个不同配置的运行(优化此处不相关的某些参数)。

目前该函数平均运行 50 个时钟。我不确定这可以改进多少，但我很高兴看到它在 30 中运行。

重点观察

如果在计算中的某个时刻，您可以看出将返回的值很小(确切的值可以协商——例如，低于一百万)，您可以提前中止。我只对大值感兴趣。

这是我希望节省最多时间的方式，而不是通过进一步的微优化(尽管这些当然也受欢迎!)。

性能信息

smallprimes 是一个位数组(64 位)；平均将设置大约 8 位，但可以少至 0 或多至 12。

q 通常为非零值。 (请注意，如果 q 和 smallprimes 为零，函数会提前退出。)

r 和 s 通常为 0。如果 q 为零，则 r 和 s 也将为；如果 r 为零，则 s 也为零。

正如最后的评论所说， nu 通常在最后是 1，所以我有一个有效的特例。

特殊情况下的计算可能看起来有溢出的风险，但通过适当的建模，我已经证明，对于我的输入，这不会发生——所以不要担心这种情况。

此处未定义的功能(ugcd、minuu、star等)已经优化；没有人需要很长时间才能运行。 pr 是一个小数组(都在 L1 中)。此外，这里调用的所有函数都是 pure functions .

但如果你真的在乎... ugcd 是 gcd , minuu 是最小值，vals 是尾随二进制 0 的个数，__builtin_ffs 是最左边的二进制 1 的位置，star 是 (n-1) >> vals(n-1)，pr 是一个由 2 中的素数组成的数组到 313。

目前正在 Phenom II 920 x4 上进行计算，但仍然对 i7 或 Woodcrest 的优化感兴趣(如果我在其他节点上获得计算时间)。

我很乐意回答您对该功能或其组成部分的任何问题。

它实际上做了什么

添加以响应请求。你不需要阅读这部分。

输入是一个奇数 n，其中 1 < n < 4282250400097。其他输入提供了这个特定意义上的数字因式分解:

如果数字可以被 3 整除，则设置 smallprimes&1，如果数字可以被 5 整除，则设置 smallprimes&2，如果数字可以被 7 整除，则设置 smallprimes&4，如果数字可以被 11 整除，则设置 smallprimes&8，等等。表示 313 的有效位。可被素数平方整除的数与仅可被该数整除的数的表示方式不同。 (实际上，平方倍数可以丢弃；在另一个函数的预处理阶段，素数平方倍数 <= lim 有 smallprimes 和 q 设置为 0，因此它们将被丢弃，其中 lim 的最佳值由实验确定。 )

q、r 和 s 表示数字的较大因数。任何剩余的因子(可能大于数字的平方根，或者如果 s 非零甚至可能更小)可以通过从 n 中除以因子来找到。

一旦以这种方式恢复所有因子，碱基数，1 <= b < n，其中 n 为 strong pseudoprime使用代码最好解释的数学公式计算。

到目前为止的改进

插入提前退出测试。这显然节省了工作，所以我进行了更改。

适当的函数已经内联，所以 __attribute__ ((inline))什么也没做。奇怪的是，标记主函数bases以及一些 __attribute ((hot)) 的 helper 性能下降了近 2%，我不知道为什么(但它可以通过 20 多次测试重现)。所以我没有做那个改变。同样，__attribute__ ((const)) ，充其量也没有帮助。我对此感到非常惊讶。

代码

ulong bases(ulong smallprimes, ulong n, ulong q, ulong r, ulong s)
{
    if (!smallprimes & !q)
        return 0;

    ulong f = __builtin_popcountll(smallprimes) + (q > 1) + (r > 1) + (s > 1);
    ulong nu = 0xFFFF;  // "Infinity" for the purpose of minimum
    ulong nn = star(n);
    ulong prod = 1;

    while (smallprimes) {
        ulong bit = smallprimes & (-smallprimes);
        ulong p = pr[__builtin_ffsll(bit)];
        nu = minuu(nu, vals(p - 1));
        prod *= ugcd(nn, star(p));
        n /= p;
        while (n % p == 0)
            n /= p;
        smallprimes ^= bit;
    }
    if (q) {
        nu = minuu(nu, vals(q - 1));
        prod *= ugcd(nn, star(q));
        n /= q;
        while (n % q == 0)
            n /= q;
    } else {
        goto BASES_END;
    }
    if (r) {
        nu = minuu(nu, vals(r - 1));
        prod *= ugcd(nn, star(r));
        n /= r;
        while (n % r == 0)
            n /= r;
    } else {
        goto BASES_END;
    }
    if (s) {
        nu = minuu(nu, vals(s - 1));
        prod *= ugcd(nn, star(s));
        n /= s;
        while (n % s == 0)
            n /= s;
    }

    BASES_END:
    if (n > 1) {
        nu = minuu(nu, vals(n - 1));
        prod *= ugcd(nn, star(n));
        f++;
    }

    // This happens ~88% of the time in my tests, so special-case it.
    if (nu == 1)
        return prod << 1;

    ulong tmp = f * nu;
    long fac = 1 << tmp;
    fac = (fac - 1) / ((1 << f) - 1) + 1;
    return fac * prod;
}

Answer 1

您似乎浪费了很多时间按因素进行除法。用除数的倒数(除数:~15-80( ! ) 个周期，取决于除数，乘法:~4 个周期)、 的倒数替换除法要快得多如果 当然，您可以预先计算倒数。

虽然这对于 q、r、s 似乎不太可能 - 由于这些变量的范围，但很容易使用 p，它总是来自小的静态 pr[] 数组。预先计算这些素数的倒数并将它们存储在另一个数组中。然后，不是除以 p，而是乘以取自第二个数组的倒数。 (或制作单个结构数组。)

现在，通过这种方法获得精确的除法结果需要一些技巧来补偿舍入误差。您可以在 this document 中找到该技术的详细信息。 ，第 138 页。

编辑:

在就该主题咨询 Hacker's Delight(一本优秀的书，顺便说一句)之后，似乎您可以通过利用代码中的所有除法都是精确的(即余数为零)这一事实来使其更快。

似乎对于每个为奇数且基数为 B = 2word_size 的除数 d，都存在唯一满足条件的乘法逆 d⃰:d⃰ < B和 d·d⃰ ≡ 1 (mod B) .对于每个是 d 的精确倍数的 x，这意味着 x/d ≡ x·d⃰ (mod B) .这意味着您可以简单地用乘法替换除法，无需添加更正、检查、舍入问题等。 (这些定理的证明可以在书中找到。)备注 这个乘法逆不需要等于前面方法定义的倒数!

如何检查给定的 x 是否是 d 的精确倍数 - 即 x mod d = 0 ?简单! x mod d = 0伊夫x·d⃰ mod B ≤ ⌊(B-1)/d⌋ .请注意，此上限可以预先计算。

所以，在代码中:

unsigned x, d;
unsigned inv_d = mulinv(d);          //precompute this!
unsigned limit = (unsigned)-1 / d;   //precompute this!

unsigned q = x*inv_d;
if(q <= limit)
{
   //x % d == 0
   //q == x/d
} else {
   //x % d != 0
   //q is garbage
}

假设 pr[]数组变成了 struct prime 的数组:

struct prime {
   ulong p;
   ulong inv_p;  //equal to mulinv(p)
   ulong limit;  //equal to (ulong)-1 / p
}

while(smallprimes)代码中的循环变为:

while (smallprimes) {
    ulong bit = smallprimes & (-smallprimes);
    int bit_ix = __builtin_ffsll(bit);
    ulong p = pr[bit_ix].p;
    ulong inv_p = pr[bit_ix].inv_p;
    ulong limit = pr[bit_ix].limit;
    nu = minuu(nu, vals(p - 1));
    prod *= ugcd(nn, star(p));
    n *= inv_p;
    for(;;) {
        ulong q = n * inv_p;
        if (q > limit)
            break;
        n = q;
    }
    smallprimes ^= bit;
}

而对于 mulinv()功能:

ulong mulinv(ulong d) //d needs to be odd
{
   ulong x = d;
   for(;;)
   {
      ulong tmp = d * x;
      if(tmp == 1)
         return x;
      x *= 2 - tmp;
   }
}

注意你可以替换 ulong与任何其他无符号类型 - 只需始终使用相同的类型。

书中提供了证明、原因和方法。强烈推荐阅读:-)。