优化我！（C，性能）--bit-twiddling 问题的后续行动答案

【问题标题】：Optimize me! (C, performance) -- followup to bit-twiddling question优化我！（C，性能）--bit-twiddling 问题的后续行动
【发布时间】：2013-08-03 20:23:40
【问题描述】：

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

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

概述

此函数将至少被调用 10¹³ 次，并且可能多达 10¹⁵ 次。也就是说，这段代码很可能会运行 几个月，所以任何性能提示都会有所帮助。

此函数占程序时间的 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

如果数字能被 3 整除，则设置 smallprimes&1，如果数字能被 5 整除，则设置 smallprimes&2，如果数字能被 7 整除，则设置 smallprimes&4，如果数字能被 11 整除，则设置 smallprimes&8，依此类推。到表示 313 的最高有效位。可被素数的平方整除的数与仅可被该数整除的数的表示方式并无不同。（实际上，可以丢弃平方倍数；在另一个函数的预处理阶段，素数

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

以这种方式恢复所有因素后，使用代码最能解释的数学公式计算碱基数 1 strong pseudoprime。

目前的改进

提前退出测试。这显然可以节省工作量，因此我进行了更改。
适当的函数已经内联，所以__attribute__ ((inline)) 什么都不做。奇怪的是，标记主要功能 bases 和一些带有 __attribute ((hot)) 的助手会降低近 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;
}

【问题讨论】：

你能解释一下代码实际上试图计算什么吗？我在数论方面很糟糕，但在数论方面做得更好的人可能会看到你可以做出的一些算法改进。
有代表的人真的需要创建一个成熟的优化标签，为这种情况保留......
@slacker：是的……但不是。我提到计算1 << tmp 在这个实现中保证不会溢出，但是我编写了代码来处理这种情况。因为现在不适用该代码已被注释掉，但我将其留在原处以防我稍后对其进行概括。不过不错！
谢谢大家！我现在已经将程序运行到了内置格式的限制。我正在进行部分重写以使其更高。
只是一个侧节点，目前有一个用于代码审查的 beta 堆栈交换站点，您可以查看codereview.stackexchange.com

标签： c performance optimization math bit-manipulation

【解决方案1】：

您似乎在按因子进行除法浪费了很多时间。用除数的倒数乘以代替除法要快得多（除法：~15-80（！）个周期，取决于除数，乘法：~4个周期），如果当然你可以预先计算倒数。

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

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

编辑：

在咨询了Hacker's Delight（一本优秀的书，顺便说一句）之后，您似乎可以通过利用代码中的所有除法都是精确的事实（即余数是零）。

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

如何检查给定的 x 是否是 d 的精确倍数 - 即 x mod d = 0 ？简单的！ x mod d = 0 iff 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 替换为任何其他无符号类型 - 只需始终使用相同的类型即可。

证明、为什么和如何都可以在书中找到。强烈推荐阅读:-)。

【讨论】：

谢谢！现在正在调查。

【解决方案2】：

如果你的编译器支持 GCC 函数属性，你可以用这个属性标记你的纯函数：

ulong star(ulong n) __attribute__ ((const));

此属性向编译器表明函数的结果仅取决于其参数。优化器可以使用此信息。

您打开编码 vals() 而不是使用 __builtin_ctz() 有什么原因吗？

【讨论】：

@caf: 不是还有pure 属性吗？
@Jens Gustedt：是的，但是pure 没有const 严格（pure 函数可以访问全局变量，const 不允许）。
好点。我想我也会输入__attribute__ ((hot))。
我认为所有这些属性只对线外功能很重要。你真正想要的是always_inline
小心always_inline 往往会适得其反。优化器通常比程序员更了解 I$ 未命中和调用开销之间的权衡。我有一位一直在使用的同事，我尝试了他的代码选项-fno-inline，它甚至覆盖了always_inline，结果程序速度提高了 10%（并且在 SPARC 上使用昂贵的窗口寄存器调用）。

【解决方案3】：

目前还不清楚，您在搜索什么。很多时候，数论问题通过推导解决方案必须满足的数学属性来实现巨大的加速。

如果您确实在搜索使 MR 测试的非见证人数量最大化的整数（即您提到的 oeis.org/classic/A141768），那么可以使用非见证人的数量不能大于 phi(n)/4 并且有这么多非见证的整数是两个素数的乘积，形式为

(k+1)*(2k+1)

或者它们是具有 3 个质因数的卡迈克尔数。我认为序列中的所有整数都具有这种形式，并且可以通过证明所有其他整数的见证人的上限来验证这一点。例如。具有 4 个或更多因数的整数总是最多有 phi(n)/8 个非见证。类似的结果可以从你的其他整数的基数公式中得出。

至于微优化：只要您知道一个整数可以被某个商整除，那么就可以将除法替换为与模 2^64 取反的商的乘法。并且测试 n % q == 0 可以替换为测试

n * inverse_q

其中 inverse_q = q^(-1) mod 2^64 和 max_q = 2^64 / q。显然 inverse_q 和 max_q 需要预先计算，以提高效率，但由于您使用的是筛子，我认为这不应该成为障碍。

【讨论】：

我会非常认真地调查此事。这几天我身体不太好；很抱歉似乎忽略了您深思熟虑的评论。
@Charles 我倾向于认为运行程序来查找序列 A141768 的值是浪费时间。但我也认为您的程序可用于查找通过 Miller-Rabin 测试的随机 k 位整数的平均碱基数。将结果与 Damgård、Landrock、Pomerance 的论文“强可能素数检验的平均案例误差估计”进行比较会很有趣。但是，当然这只是我个人的看法。
当然，搜索的部分灵感来自 DLP。但我不认为这种方法对于检查他们的数字有多大用处，因为感兴趣的区域可能是 300 到 3000 位，其中分解是困难的并且不可能进行筛选。

【解决方案4】：

小优化但是：

ulong f;
ulong nn;
ulong nu = 0xFFFF;  // "Infinity" for the purpose of minimum
ulong prod = 1;

if (!smallprimes & !q)
    return 0;

// no need to do this operations before because of the previous return
f = __builtin_popcountll(smallprimes) + (q > 1) + (r > 1) + (s > 1);
nn = star(n);

顺便说一句：您应该编辑您的帖子以添加 star() 和您使用定义的其他功能

【讨论】：

我在您输入功能时编辑了我的问题。确切的代码不应该太重要——它们已经被优化了。但也许我还是会发布它们；我只是不想用太长的代码块把任何人赶走。
哦...这似乎是一个小优化（错过它对我来说当然是愚蠢的！），但在我的计算生命周期中，我估计这将节省 200 万亿个周期。

【解决方案5】：

尝试替换此模式（也适用于 r 和 q）：

n /= p; 
while (n % p == 0) 
  n /= p;

有了这个：

ulong m;
  ...
m = n / p; 
do { 
  n = m; 
  m = n / p; 
} while ( m * p == n);

在我的有限测试中，我通过消除模数获得了小幅加速 (10%)。

此外，如果 p、q 或 r 为常数，编译器将用乘法替换除法。如果 p、q 或 r 的选择很少，或者某些选项更频繁，则可以通过专门针对这些值的函数来获得一些好处。

【讨论】：

我会试试的。不幸的是，p、q 或 r 中没有一个是恒定的，并且值可能变化很大（从 317 到几百万）。
同时拥有一个汇编 DIV 指令用于除法（商）和模（余数）会比除法和乘法更快。
我还建议避免使用 %，它与除法一样昂贵，而乘法通常很快。

【解决方案6】：

您是否尝试过使用配置文件引导优化？

使用-fprofile-generate 选项编译和链接程序，然后在代表性数据集上运行程序（例如，一天的计算量）。

然后重新编译并将其与-fprofile-use 选项链接。

【讨论】：

【解决方案7】：

1) 我会让编译器吐出它生成的程序集，并尝试推断它所做的是否是它所能做的最好的......如果发现问题，请更改代码以使程序集看起来更好。这样，您还可以确保您希望它内联的函数（如星号和 vals）是真正内联的。（您可能需要添加编译指示，甚至将它们转换为宏）

2) 在多核机器上尝试这个很棒，但是这个循环是单线程的。我猜有一个伞式函数可以将负载分散到几个线程上，以便使用更多的内核？

3) 如果实际函数试图计算的内容不清楚，则很难建议加快速度。通常，最令人印象深刻的加速不是通过比特旋转来实现的，而是通过算法的改变来实现的。所以一点 cmets 可能会有所帮助;^)

4) 如果您真的想要 10* 或更高的速度，请查看 CUDA 或 openCL，它允许您在图形硬件上运行 C 程序。它闪耀着这样的功能！

5）您正在做大量的模数运算，然后彼此分开。在 C 中，这是 2 个单独的命令（首先是“/”，然后是“%”）。但是在汇编中，这是 1 个命令：'DIV' 或 'IDIV' 一次性返回余数和商：

B.4.75 IDIV: Signed Integer Divide

IDIV r/m8                     ; F6 /7                [8086]
IDIV r/m16                    ; o16 F7 /7            [8086]
IDIV r/m32                    ; o32 F7 /7            [386]

IDIV performs signed integer division. The explicit operand provided is the divisor; the dividend and destination operands are implicit, in the following way:

For IDIV r/m8, AX is divided by the given operand; the quotient is stored in AL and the remainder in AH.

For IDIV r/m16, DX:AX is divided by the given operand; the quotient is stored in AX and the remainder in DX.

For IDIV r/m32, EDX:EAX is divided by the given operand; the quotient is stored in EAX and the remainder in EDX.

因此它需要一些内联汇编，但我猜这会显着加快速度，因为您的代码中有几个地方可以从中受益。

【讨论】：

这些都是很好的建议 (+1)，但我真的不认为它们对我有多大帮助。对于#4，CUDA 在这项任务上会非常出色，但是为它重写程序并非易事，因为进入 GPU 的信息通过一个小管道传输（因此为了有效地完成它，我必须移动过去很多我的程序，这里没有看到，到 CUDA）。而且，这需要很多我没有的技能！对于#3，我花了很长时间为这个特定任务设计算法；我比程序员更像数学家。您在 #2 中的猜测是准确的——我将运行该程序的四个副本。
#1 是个好主意，当然值得。我的组装有点生锈了，但如果我卡住了，我有一本很好的手册。如果我在这方面有任何进展，我会通知您。
@charles 不要低估自己。 Cuda 运行 C 程序，因此您已经拥有很多技能。确实，它需要对数据流进行一些重新架构。这当然不是微不足道的，但回报会非常甜蜜。（10* 加速是典型的......两倍的速度并不少见）
我将查看生成的程序集。我认为 gcc 足够聪明，在这种情况下只发出一个 IDIV，但也许不是。
Um boyz，div 在 C 的标准库中，不需要重新实现。我想 GNU-C 能够以这种方式将其编译为内在的。 gnu.org/software/libtool/manual/libc/Integer-Division.html

【解决方案8】：

确保您的函数被内联。如果它们不合规，则开销可能会增加，尤其是在第一个 while 循环中。最好的确定方法是检查程序集。
您是否尝试过预计算 star( pr[__builtin_ffsll(bit)] ) 和 vals( pr[__builtin_ffsll(bit)] - 1) ？这将用一些简单的工作换取数组查找，但如果表足够小，这可能是值得的。
在你真正需要它之前不要计算f（接近尾声，在你提前退出之后）。您可以将 BASES_END 周围的代码替换为

BASES_END:
ulong addToF = 0;
if (n > 1) {
    nu = minuu(nu, vals(n - 1));
    prod *= ugcd(nn, star(n));
    addToF = 1;
}
// ... early out if nu == 1...
// ... compute f ...
f += addToF;

希望对您有所帮助。

【讨论】：

非常有用的建议！是的，这些肯定值得预先计算。
smallprimes 循环破坏了 smallprimes 中的值，所以我仍然需要在顶部执行 popcountll(smallprimes)。但将其余部分向下移动应该可以为我节省 50 万亿个周期。
确保你测试得很好。为缓存未命中交易一些位移并不是一个好的交易。
或者你可以复制 smallprimes 的原始值。我不确定 popcountll 实际编译成什么。
如果你告诉编译器使用最新最好的指令集扩展（我假设 OP 正在这样做，考虑到他对性能的需求），那么popcountll() 将编译为单个 POPCNT 指令，它AMD 需要 2 个时钟周期，Intel 需要 3 个时钟周期。

【解决方案9】：

首先进行一些吹毛求疵 ;-) 你应该更加小心你正在使用的类型。在某些地方，您似乎认为 ulong 是 64 位宽，请在此处使用 uint64_t。对于所有其他类型，请仔细重新考虑您对它们的期望并使用适当的类型。

我可以看到的优化是整数除法。你的代码做了很多，这可能是你正在做的最昂贵的事情。小整数 (uint32_t) 的除法可能比大整数更有效。特别是对于uint32_t，有一条汇编指令可以一次性完成除法和取模操作，称为divl。

如果您使用适当的类型，您的编译器可能会为您完成这一切。但是你最好检查一下汇编器（选项-S 到 gcc），正如有人已经说过的那样。否则很容易在这里和那里包含一些小的汇编程序片段。我在我的一些代码中发现了类似的东西：

register uint32_t a asm("eax") = 0;
register uint32_t ret asm("edx") = 0;
asm("divl %4"
    : "=a" (a), "=d" (ret)
    : "0" (a), "1" (ret), "rm" (divisor));

如您所见，它使用了特殊寄存器 eax 和 edx 以及类似的东西......

【讨论】：

在这个应用程序中，有 #ifdef 语句保证除了 64 位的情况之外永远不会看到这段代码（我认为 ulong 是一个 typedef）。
@Charles：好的，但这只是问题的一部分。在不需要全宽的地方使用较小的类型。首先，这可能会释放优化器的压力，因为寄存器仍然是稀缺资源。然后，对于我的示例中的ldiv，可能有更快的汇编器操作可用于较小的类型。这样做时，更容易阅读 uint32_t 与 uint64_t 然后 ulong ;-)
我可以保证 q、r 和 s 可以放入 uint32_t，所以在这种情况下，这可能会节省一些时间。（我会检查实际性能差异并回复您。）在ulong上，我只是按照其他200,000 LOC的编码风格......

【解决方案10】：

您是否尝试过第一个 while 循环的表查找版本？您可以将smallprimes 划分为 4 个 16 位值，查找它们的贡献并将它们合并。但也许你需要副作用。

【讨论】：

没有副作用本身，从这里调用的所有函数都是纯函数。因此，也许可以使用多个表进行查找……不确定。另外：我把这个澄清编辑到我的问题中，因为我最初并不清楚。

【解决方案11】：

您是否尝试过传入素数数组而不是将它们拆分为smallprimes、q、r 和s？由于我不知道外部代码的作用，我可能是错的，但是您也有可能将一些素数转换为smallprimes 位图，在此函数内部，您将位图转换回一组素数，有效。此外，您似乎对smallprimes、q、r 和s 的元素进行了相同的处理。它应该为您节省每次调用的少量处理。

此外，您似乎知道传入的素数除以n。你对除 n 的每个素数的幂有足够的了解吗？如果您可以通过将该信息传递给此函数来消除模运算，则可以节省大量时间。换句话说，如果n 是pow(p_0,e_0)*pow(p_1,e_1)*...*pow(p_k,e_k)*n_leftover，并且如果您对这些e_is 和n_leftover 有更多的了解，那么将它们传入将意味着您在此函数中不必做很多事情。

可能有一种方法可以通过更少的模运算来发现n_leftover（n 的未分解部分），但这只是一种预感，因此您可能需要稍微尝试一下。这个想法是使用gcd 反复从 n 中删除已知因子，直到你摆脱所有已知的素因子。让我给出一些几乎是c的代码：

factors=p_0*p_1*...*p_k*q*r*s;
n_leftover=n/factors;
do {
    factors=gcd(n_leftover, factors);
    n_leftover = n_leftover/factors;
} while (factors != 1);

我完全不确定这会比您拥有的代码更好，更不用说您可以在其他答案中找到的组合 mod/div 建议，但我认为值得一试。我觉得这将是一场胜利，尤其是对于具有大量小质因数的数字。

【讨论】：

这听起来对我来说也是一个胜利——与其传递位图，不如传递一个指向数组的指针，其中数组的每个元素都是相应素数的幂在n 的因式分解中。
是的，我真的做不到。我已经花费 2 GB 将数字存储为位图，如果我将它们存储为数组，我无法在内存中存储尽可能多的数字，所以我不得不花更多的时间进行筛选。这将是一笔巨大的成本。
很公平。如果内存是个问题，那么使用位图绝对是必要的。
我的答案中的筛子实际上产生了一个除以 n 的素数链表，而不是指数。它只消耗 4MB，所以它应该适合缓存 :-)
你似乎已经感觉到我在想什么——关于缓存，但我想我没有写下来。 :) 无论如何，鉴于您不知道指数，我添加了一个算法，可以使用 gcd 而不是模（编辑答案）发现未分解的部分。您可能想尝试一下。

【解决方案12】：

您传递的是 n 的完整分解，因此您要分解连续整数，然后在此处使用该分解的结果。在我看来，您可能会在找到这些因素时从这样做中受益。

顺便说一句，我有一些非常快速的代码可以在不进行任何除法的情况下找到您正在使用的因子。它有点像筛子，但会很快非常产生连续数字的因子。如果您认为有帮助，可以找到并发布。

edit 不得不在这里重新创建代码：

＃包括 #define SIZE (1024*1024) //必须是 2^n #define 面具 (SIZE-1) 类型定义结构{ 诠释 p; 下一个； } p_type; p_type 素数[大小]; int筛子[大小]; 无效的 init_sieve() { 诠释我，n; 整数计数 = 1；素数[1].p = 3; 筛子[1] = 1; 对于 (n=5;尺寸>n;n+=2) { 整数标志 = 0；对于 (i=1;count>=i;i++) { 如果 ((n%primes[i].p) == 0) { 标志 = 1; 休息; } } 如果（标志==0） { 计数++；素数[计数].p = n; 筛[n>>1] =计数； } } } 主函数（） { 整数点，n； init_sieve(); printf("init_done\n"); // 分解以 3 开头的奇数对于 (n=1;1000000000>n;n++) { ptr = 筛子[n&MASK]; if (ptr == 0) //素数 { // printf("%d 是素数",n*2+1); } 否则//复合 { // printf ("%d 有除数:",n*2+1); 而（指针！= 0） { // printf ("%d",primes[ptr].p); 筛子[n&MASK]=primes[ptr].next; //将素数移动到它所除的下一个数字 primes[ptr].next = 筛子[(n+primes[ptr].p)&MASK]; 筛子[(n+primes[ptr].p)&MASK] = ptr; ptr = 筛子[n&MASK]; } } // printf("\n"); } 返回0； }

init 函数创建一个因子基并初始化筛子。这在我的笔记本电脑上大约需要 13 秒。然后再过 25 秒，所有高达 10 亿的数字都被分解或确定为素数。小于 SIZE 的数字永远不会报告为素数，因为它们在因子基数中有 1 个因子，但可以更改。

这个想法是为筛子中的每个条目维护一个链表。通过简单地将它们的因子从链表中拉出来对数字进行因子分解。当它们被拉出时，它们被插入到下一个可以被那个素数整除的数字的列表中。这对缓存也非常友好。筛子大小必须大于因子基中的最大素数。事实上，这个筛子可以在大约 7 小时内运行到 2**40，这似乎是你的目标（除了 n 需要是 64 位）。

您的算法可以合并到其中，以便在识别时利用这些因素，而不是将位和大素数打包到变量中以传递给您的函数。或者您的函数可以更改为采用链表（您可以创建一个虚拟链接以传递因子基数之外的素数）。

希望对你有帮助。

顺便说一句，这是我第一次公开发布这个算法。

【讨论】：

我正在筛选生成因子，是的。（它实际上是连续的奇数——偶数与这个问题无关。）程序最终只花费了大约四分之一的时间在筛子上，其余的时间在这里。但是，如果您的代码很聪明，请随时发布。

【解决方案13】：

只是一个想法，但如果您还没有使用编译器优化选项，也许会有所帮助。另一个想法是，如果钱不是问题，您可以使用英特尔 C/C++ 编译器，假设您使用的是英特尔处理器。我还假设其他处理器制造商（AMD 等）会有类似的编译器

【讨论】：

我目前使用的关键优化是-O3 -fno-strict-aliasing -fomit-frame-pointer。为了不破坏程序的其余部分，进一步的编译器优化可能必须是函数属性。 icc 是不可能的，因为我使用的是 Phenom II 920（你可能在我上面的文字墙中错过了这个......）。
@Charles：你有没有告诉编译器专门针对这个 CPU 模型进行优化？一个小小的-march=amdfam10就能创造奇迹。
@Charles：请注意，如果您不使用 -march 指定最小架构，编译器将不会发出任何旧 CPU 不支持的指令 - 您只能使用 SSE2 x86-64。这意味着没有 POPCNT 或 LZCNT 等。
如果您不想弄乱程序的其余部分，您不能将函数本身放入库中并最大限度地优化库，然后正常编译程序的其余部分我也从 AMD 那里找到了这个，但我不知道它是如何工作的，因为我不使用 AMD 芯片。 developer.amd.com/cpu/open64/pages/default.aspx
我相信我正在指定架构，但我会检查一下。如果不是这样，这可以节省我一些时间！

【解决方案14】：

如果您要在(!smallprimes&!q) 上立即退出，为什么不在调用函数之前进行测试，并节省函数调用开销？

此外，您似乎实际上拥有 3 个不同的线性函数，除了 smallprimes 循环。 bases1(s,n,q)、bases2(s,n,q,r) 和 bases3(s,n,q,r,s)。

实际上将它们创建为 3 个没有分支和 goto 的独立函数并调用适当的函数可能是一种胜利：

if (!(smallprimes|q)) { r = 0;}
else if (s) { r = bases3(s,n,q,r,s);}
else if (r) { r = bases2(s,n,q,r); }
else        { r = bases1(s,n,q);

如果之前的处理已经为调用代码提供了一些关于要执行哪个函数的“知识”并且您不必对其进行测试，这将是最有效的。

【讨论】：

没有函数开销；该函数将被内联（-finline-functions-call-once 包含在 -O1 及更高版本中）。另一个建议会很好，只是调用函数不会有任何知识。

【解决方案15】：

如果您使用的除法数字在编译时未知，但在运行时经常使用（多次除以相同的数字），那么我建议使用 libdivide 库，它基本上在运行时实现编译器对编译时间常数的优化（使用移位掩码等）。这可以提供巨大的好处。此外，避免将 x % y == 0 用于 z = x/y, z * y == x as ergosys 以上建议也应该有可衡量的改进。

【讨论】：

这个想法+1，但这无济于事——它有很多数字，范围很广，每个数字都使用一次。
啊，真可惜。哦，好吧，libdivide 仍然是您了解的有用工具。作者也有一个很棒的博客：ridiculousfish.com.

【解决方案16】：

您顶帖上的代码是优化版本吗？如果是，仍然有太多的除法操作，极大地消耗了 CPU 周期。

这段代码有点不必要地过度执行

if (!smallprimes & !q)
    return 0;

改为逻辑和&&

if (!smallprimes && !q)
    return 0;

在不评估 q 的情况下使其更快地短路

还有下面的代码

ulong bit = smallprimes & (-smallprimes);
ulong p = pr[__builtin_ffsll(bit)];

用于查找小素数的最后一组。你为什么不使用更简单的方法

ulong p = pr[__builtin_ctz(smallprimes)];

另一个导致性能下降的罪魁祸首可能是程序分支过多。您可以考虑更改为其他一些更少分支或无分支的等价物

【讨论】：

!smallprimes && !q 中的附加分支实际上使其比!smallprimes & !q 更昂贵。至于另一个，我在代码后面需要bit 的值，如果我直接在smallprimes 上进行查询，则数组的长度为2^64，因此不适合内存。