优化我！ （C，表现） – 跟随苦涩的问题

感谢bit Twiddling的一些非常有用的stackOverflow用户：设置了哪个位？ ，我已经构建了我的函数（在问题的最后发布）。

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

概观

此function将被调用至少10 ¹³次，并且可能经常被调用10 ^15次 。 也就是说，此代码很可能会运行数月 ，因此任何性能提示都会有所帮助。

此function占计划时间的72-77％，基于分析和不同配置中的大约十二次运行（优化此处不相关的某些参数）。

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

重点观察

如果在计算的某个时刻你可以告诉你将返回的值很小（准确值可协商 – 比如说，低于一百万） 你可以提前中止 。 我只对大价值感兴趣。

这就是我希望节省大部分时间的方式，而不是通过进一步的微观优化（尽管这些当然也是受欢迎的！）。

绩效信息

• smallprimes是一个位数组（64位）; 平均大约8位将被设置，但它可以少至0或多达12。
• q通常是非零的。 （请注意，如果q和smallprimes为零，则函数会提前退出。）
• r和s通常为0.如果q为零，r和s也将为0; 如果r为零，则s也是如此。
• 正如最后的评论所说，nu到底通常是1，所以我有一个有效的特殊情况。
• 特殊情况下的计算可能会出现溢出风险，但通过适当的建模我已经certificate，对于我的输入，这不会发生 – 所以不要担心这种情况。
• 此处未定义的function（ugcd，minuu，star等）已经过优化; 无需花很长时间才能运行。 pr是一个小数组（全部在L1中）。 此外，这里调用的所有函数都是纯函数 。
• 但是如果你真的在乎… ugcd是gcd ，minuu是最小值，vals是尾随二进制0的数量，__ builtin_ffs是最左边的二进制1的位置，star是（n-1）>> vals（n- 1），pr是从2到313的素数数组。
• 目前正在Phenom II 920 x4上进行计算，尽管对i7或Woodcrest的优化仍然很有意义（如果我在其他节点上获得计算时间）。
• 我很乐意回答您对该function或其成员的任何问题。

实际上是做什么的

添加以响应请求。 您无需阅读此部分。

输入是奇数n，其中1 <n <4282250400097.其他输入提供了此特定意义上的数字因子分解：

如果数字可被3整除，则设置smallprimes＆1;如果数字可被5整除，则设置smallprimes＆2;如果数字可被7整除，则设置smallprimes＆4;如果数字可被11整除，则设置smallprimes＆8，等等到最多表示313的有效位。可以被素数的平方整除的数字与仅被该数字整除的数字表示不同。 （实际上，可以丢弃多个正方形;在另一个函数的预处理阶段，素数的多个平方<= lim具有小的primes而q设置为0因此它们将被丢弃，其中lim的最佳值通过实验确定。 ）

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

一旦以这种方式恢复所有因子，使用由代码最佳解释的数学公式计算n为强伪荧光的碱基数1 <= b <n。

到目前为止的改进

• 推迟提前退出测试。 这显然可以节省工作，所以我做了改变。
• 相应的函数已经内联，因此__attribute__ ((inline))不执行任何操作。 奇怪的是，用__attribute ((hot))标记主要functionbases和一些助手会使性能下降近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个周期）， IF当然可以预先计算倒数。

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

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

编辑：

在咨询了Hacker’s Delight （一本优秀的书，BTW）之后，通过利用代码中所有分区都是精确的（即余数为零）这一事实，你似乎可以更快地完成它。

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

如何检查给定的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替换为任何其他无符号类型 - 只需始终使用相同的类型。

书中提供了证据， 原因和方法 。 衷心推荐阅读:-)。

Answer 2

如果您的编译器支持GCC函数属性，则可以使用以下属性标记纯函数：

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

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

您是否有理由使用开放式编码的vals()而不是使用__builtin_ctz() ？

Answer 3

它仍然有点不清楚，你在寻找什么。 数值理论问题经常通过导出解决方案必须满足的数学特性来实现巨大的加速。

如果您确实正在搜索最大化MR测试的非证人数量的整数（即您提到的oeis.org/classic/A141768）那么可能会使用非证人的数量不能更大比phi（n）/ 4并且有这么多非证人的整数要么是forms的两个素数的乘积

第（k + 1）*（2K + 1）

或者他们是Carmichael数字，有3个主要因素。 我认为在上面的某些限制中，序列中的所有整数都具有这种forms，并且可以通过certificate所有其他整数的证人的上限来validation这一点。 例如，具有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需要预先计算，才能有效，但由于你使用筛子，我认为这不应该是一个障碍。

Answer 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); 

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

Answer 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几乎没有选择，或者某些选项更频繁，您可以通过专门设置这些值的函数来获得一些东西。

Answer 6

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

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

3）如果实际函数试图计算的内容不清楚，则很难建议加速。 通常情况下，最令人印象深刻的加速并不是通过比特琐事来实现的，而是随着算法的变化而实现的。 所以一些评论可能会有所帮助; ^）

4）如果您真的想要加速10 *或更高，请查看CUDA或openCL，它允许您在图形硬件上运行C程序。 它充满了这些function！

5）你正在做大量的模数并相互分开。 在C中，这是两个单独的命令（第一个’/’然后’％’）。 但是在汇编中这是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. 

所以它需要一些内联汇编，但我猜测会有一个显着的加速，因为你的代码中有一些地方可以从中受益。

Answer 7

1. 确保您的函数内联。 如果它们不在线，则开销可能会增加，尤其是在第一个while循环中。 最好的方法是检查assembly。

2. 你有没有试过预计star( pr[__builtin_ffsll(bit)] )和vals( pr[__builtin_ffsll(bit)] - 1) ？ 这将为数组查找交换一些简单的工作，但如果表足够小，它可能是值得的。

3. 在你真正需要它之前不要计算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; 

希望有所帮助。

Answer 8

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

我能看到的优化是整数除法。 你的代码做了很多，这可能是你正在做的最昂贵的事情。 小整数划分（ uint32_t ）可能比大整数划分更有效。 特别是对于uint32_t有一个汇编指令，它一次执行除法和模divl ，称为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以及类似的东西……

Answer 9

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

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

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

Answer 10

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

Answer 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_i和n_leftover更多了解，那么将它们传入意味着你在这个function中不需要做很多事情。

---

可能有一种方法可以用较少的模运算来发现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建议，但我认为值得一试。 我觉得这将是一场胜利，特别是对于拥有大量小素因子的数字。

Answer 12

你传递的是n的完全因子分解，所以你要将连续的整数分解，然后在这里使用该分解的结果。 在我看来，在找到这些因素时，您可能会从中做到这一点。

顺便说一下，我有一些非常快速的代码来找到你正在使用的因素，而不进行任何划分。 它有点像筛子，但很快产生连续数字的因子。 如果您认为它可能有所帮助，可以找到并发布。

编辑必须在这里重新创建代码：

 ＃包括 
 #define SIZE（1024 * 1024）//必须为2 ^ n
 #define MASK（SIZE-1）
 typedef struct {
     int p;
     int next;
 } p_type;

 p_type primes [SIZE];
 int筛[SIZE];

 void init_sieve（）
 {
     int i，n;
     int count = 1;
    素数[1] .p = 3;
    筛[1] = 1;
     for（n = 5; SIZE> n; n + = 2）
     {
         int flag = 0;
         for（i = 1; count> = i; i ++）
         {
             if（（n％primes [i] .p）== 0）
             {
                 flag = 1;
                打破;
             }
         }
         if（flag == 0）
         {
            计数++;
             primes [count] .p = n;
            筛[n >> 1] =计数;
         }
     }
 }

 int main（）
 {
     int ptr，n;
     init_sieve（）;
    的printf（ “INIT_DONE \ n”）;

     //从3开始的奇数数
     for（n = 1; 1000000000> n; n ++）
     {
         ptr =筛[n＆MASK];
         if（ptr == 0）// prime
         {
 // printf（“％d是素数”，n * 2 + 1）;
         }
        否则//复合
         {
 // printf（“％d有除数：”，n * 2 + 1）;
            而（PTR！= 0）
             {
 // printf（“％d”，primes [ptr] .p）;
                筛[N MASK] =素数[PTR]的.next;
                 //将素数移动到它所分割的下一个数字
                 primes [ptr] .next = sieve [（n + primes [ptr] .p）＆MASK];
                筛[（n + primes [ptr] .p）＆MASK] = ptr;
                 ptr =筛[n＆MASK];
             }
         }
 // printf（“\ n”）;
     }
    返回0;
 }

init函数创建一个因子库并初始化筛子。 我的笔记本电脑大约需要13秒钟。 然后，所有高达10亿的数字都会在另外25秒内被考虑或确定为素数。 小于SIZE的数字从未报告为素数，因为它们在因子基数中有1个因子，但可以改变。

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

您的算法可以合并到此中，以便在识别它们时使用因子，而不是将位和大质数打包成变量以传递给您的函数。 或者可以更改您的函数以获取链接列表（您可以创建一个虚拟链接以传入因子库之外的素数）。

希望能帮助到你。

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

Answer 13

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

Answer 14

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

此外，您似乎有效地拥有3个不同的函数，除了smallprimes循环外，它们都是线性的。 bases1(s,n,q) ， bases2(s,n,q,r)和bases3(s,n,q,r,s) 。

实际创建那些没有分支和gotos的3个单独函数可能是一个胜利，并调用适当的函数：

 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); 

如果先前的处理已经给调用代码一些“知识”执行哪个函数并且您不必测试它，那么这将是最有效的。

Answer 15

如果您正在使用的分区具有在编译时未知但在运行时经常使用的数字（除以相同的数字多次），那么我建议使用libdivide库，它基本上在运行时实现编译器为编译时常量做的优化（使用移位掩码等）。 这可以提供巨大的好处。 同样避免使用x％y == 0来表示像z = x / y，z * y == x这样的ergosys也应该有一个可衡量的改进。

Answer 16

您post中的代码是优化版本吗？ 如果是的话，仍有太多的除法操作会大大超出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)]; 

降低性能的另一个罪魁祸首可能是程序分支过多。 您可以考虑更改为其他一些较少分支或无分支的等效项