优化我！ (C，表演)-纠缠不清的问题的跟进 [英] Optimize me! (C, performance) -- followup to bit-twiddling question
问题描述
感谢一些非常有用的stackOverflow用户，网址为>位旋转:设置了哪一位?，我已经构造了函数(在问题末尾张贴).
任何建议-甚至是很小的建议-都将受到赞赏.希望它可以使我的代码更好，但是至少它应该教给我一些东西. :)
概述
此函数将至少被调用10 13 次，并且可能多达10 15 .也就是说，此代码很可能会运行几个月，因此任何性能提示都将有所帮助.
此功能基于配置文件以及在不同配置中进行的十几次运行(优化与此处无关的某些参数)，占程序时间的72-77％.
此功能目前平均运行50个时钟.我不确定可以改善多少，但我很高兴看到它能在30天内运行.
重点观察
如果在计算的某个时刻您可以确定将返回的值很小(确切的值可以协商，例如，小于一百万)，您可以提早中止.我只对大价值感兴趣.
这是我希望节省最多时间的方式，而不是通过进一步的微优化来实现(尽管当然也欢迎这样做！).
性能信息
- smallprimes是一个位数组(64位)；平均将设置约8位，但可能少至0或多至12.
- q通常为非零. (请注意，如果q和smallprimes为零，该函数会提前退出.) r和s通常为0.如果q为零，则r和s也为0.如果r为零，则s也是如此.
- 正如结尾处的评论所述，nu通常在结尾处为1，因此我有一个有效的特殊情况.
- 特殊情况下的计算可能会出现溢出风险，但是通过适当的建模，我已经证明，对于我的输入，这种情况不会发生-因此，不必担心这种情况.
- 此处未定义的功能(ugcd，minuu，star等)已经过优化；没有一个需要很长时间才能运行. pr是一个小数组(全部在L1中).此外，此处调用的所有函数都是纯函数.
- 但是如果您真的很在意... ugcd是 gcd ，则minuu是最小值，vals是尾随二进制数0的数目，__builtin_ffs是最左边的二进制数1的位置，星号是(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；如果将数字除以7，则设置smallprimes& 8.可以被11整除，直到代表313的最高有效位.被质数的平方整除的数字与仅被该数整除的数字的表示方式不同. (实际上，可以丢弃平方的倍数;在预处理阶段，在另一个函数中，素数<= lim的平方的平方具有小素数，并且q设置为0，因此它们将被丢弃，其中lim的最佳值由实验确定)
q，r和s代表该数字的较大因素.通过将因子从n中除，可以找到任何剩余因子(可能大于数字的平方根，或者如果s不为零，甚至可能更小). 一旦以这种方式恢复了所有因子，则碱基数为1≤b≤1. n，其中n是强伪素，使用由代码最佳解释的数学公式进行计数到目前为止的改进
- 推动提前退出测试.显然可以节省工作，因此我进行了更改.
- 适当的功能已经内联，因此
__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;
}
您似乎在浪费很多时间进行因数分解.用除数(除数:〜15-80(！)周期的倒数，取决于除数，乘数:〜4个周期)，如果，您当然可以预先计算倒数.
虽然 q ， r ， s 似乎不太可能-由于这些变量的范围，这非常容易与 p 有关，它总是来自小型静态 pr [] 数组.预计算这些素数的倒数，并将它们存储在另一个数组中.然后，将其乘以第二个数组的倒数，而不是除以 p . (或制作一个单一的结构数组.)
现在，通过这种方法获得精确的除法结果需要一些技巧来补偿舍入误差.您可以在本文档(第138页)中找到该技术的详细内容.
在咨询了关于该主题的 Hacker's Delight (一本很棒的书，BTW)之后，您似乎可以利用代码中所有划分都是精确的事实(例如，余数是零).
似乎对于每个奇数且基数 B = 2 word_size 的除数 d ，都有一个唯一的乘积逆满足条件d⃰ < B和d·d⃰ ≡ 1 (mod B)的d⃰.对于每个 d 的精确倍数的 x ，这意味着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替换为任何其他无符号类型-只需始终使用相同的类型即可.
证明，为什么和如何在本书中均可用.衷心推荐阅读:-).
Thanks to some very helpful stackOverflow users at Bit twiddling: which bit is set?, I have constructed my function (posted at the end of the question).
Any suggestions -- even small suggestions -- would be appreciated. Hopefully it will make my code better, but at the least it should teach me something. :)
Overview
This function will be called at least 1013 times, and possibly as often as 1015. That is, this code will run for months in all likelihood, so any performance tips would be helpful.
This function accounts for 72-77% of the program's time, based on profiling and about a dozen runs in different configurations (optimizing certain parameters not relevant here).
At the moment the function runs in an average of 50 clocks. I'm not sure how much this can be improved, but I'd be thrilled to see it run in 30.
Key Observation
If at some point in the calculation you can tell that the value that will be returned will be small (exact value negotiable -- say, below a million) you can abort early. I'm only interested in large values.
This is how I hope to save the most time, rather than by further micro-optimizations (though these are of course welcome as well!).
Performance Information
- smallprimes is a bit array (64 bits); on average about 8 bits will be set, but it could be as few as 0 or as many as 12.
- q will usually be nonzero. (Notice that the function exits early if q and smallprimes are zero.)
- r and s will often be 0. If q is zero, r and s will be too; if r is zero, s will be too.
- As the comment at the end says, nu is usually 1 by the end, so I have an efficient special case for it.
- The calculations below the special case may appear to risk overflow, but through appropriate modeling I have proved that, for my input, this will not occur -- so don't worry about that case.
- Functions not defined here (ugcd, minuu, star, etc.) have already been optimized; none take long to run. pr is a small array (all in L1). Also, all functions called here are pure functions.
- But if you really care... ugcd is the gcd, minuu is the minimum, vals is the number of trailing binary 0s, __builtin_ffs is the location of the leftmost binary 1, star is (n-1) >> vals(n-1), pr is an array of the primes from 2 to 313.
- The calculations are currently being done on a Phenom II 920 x4, though optimizations for i7 or Woodcrest are still of interest (if I get compute time on other nodes).
- I would be happy to answer any questions you have about the function or its constituents.
What it actually does
Added in response to a request. You don't need to read this part.
The input is an odd number n with 1 < n < 4282250400097. The other inputs provide the factorization of the number in this particular sense:
smallprimes&1 is set if the number is divisible by 3, smallprimes&2 is set if the number is divisible by 5, smallprimes&4 is set if the number is divisible by 7, smallprimes&8 is set if the number is divisible by 11, etc. up to the most significant bit which represents 313. A number divisible by the square of a prime is not represented differently from a number divisible by just that number. (In fact, multiples of squares can be discarded; in the preprocessing stage in another function multiples of squares of primes <= lim have smallprimes and q set to 0 so they will be dropped, where the optimal value of lim is determined by experimentation.)
q, r, and s represent larger factors of the number. Any remaining factor (which may be greater than the square root of the number, or if s is nonzero may even be less) can be found by dividing factors out from n.
Once all the factors are recovered in this way, the number of bases, 1 <= b < n, to which n is a strong pseudoprime are counted using a mathematical formula best explained by the code.
Improvements so far
- Pushed the early exit test up. This clearly saves work so I made the change.
- The appropriate functions are already inline, so
__attribute__ ((inline))does nothing. Oddly, marking the main functionbasesand some of the helpers with__attribute ((hot))hurt performance by almost 2% and I can't figure out why (but it's reproducible with over 20 tests). So I didn't make that change. Likewise,__attribute__ ((const)), at best, did not help. I was more than slightly surprised by this.
Code
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;
}
You seem to be wasting much time doing divisions by the factors. It is much faster to replace a division with a multiplication by the reciprocal of divisor (division: ~15-80(!) cycles, depending on the divisor, multiplication: ~4 cycles), IF of course you can precompute the reciprocals.
While this seems unlikely to be possible with q, r, s - due to the range of those vars, it is very easy to do with p, which always comes from the small, static pr[] array. Precompute the reciprocals of those primes and store them in another array. Then, instead of dividing by p, multiply by the reciprocal taken from the second array. (Or make a single array of structs.)
Now, obtaining exact division result by this method requires some trickery to compensate for rounding errors. You will find the gory details of this technique in this document, on page 138.
EDIT:
After consulting Hacker's Delight (an excellent book, BTW) on the subject, it seems that you can make it even faster by exploiting the fact that all divisions in your code are exact (i.e. remainder is zero).
It seems that for every divisor d which is odd and base B = 2word_size, there exists a unique multiplicative inverse d⃰ which satisfies the conditions: d⃰ < B and d·d⃰ ≡ 1 (mod B). For every x which is an exact multiple of d, this implies x/d ≡ x·d⃰ (mod B). Which means you can simply replace a division with a multiplication, no added corrections, checks, rounding problems, whatever. (The proofs of these theorems can be found in the book.) Note that this multiplicative inverse need not be equal to the reciprocal as defined by the previous method!
How to check whether a given x is an exact multiple of d - i.e. x mod d = 0 ? Easy! x mod d = 0 iff x·d⃰ mod B ≤ ⌊(B-1)/d⌋. Note that this upper limit can be precomputed.
So, in code:
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
}
Assuming the pr[] array becomes an array of struct prime:
struct prime {
ulong p;
ulong inv_p; //equal to mulinv(p)
ulong limit; //equal to (ulong)-1 / p
}
the while(smallprimes) loop in your code becomes:
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;
}
And for the mulinv() function:
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;
}
}
Note you can replace ulong with any other unsigned type - just use the same type consistently.
The proofs, whys and hows are all available in the book. A heartily recommended read :-).
这篇关于优化我！ (C，表演)-纠缠不清的问题的跟进的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
