优化我!(C,表现) - 跟随喋喋不休的问题

Question

感谢bit Twiddling的一些非常有用的stackOverflow用户:设置了哪个位？,我已经构建了我的函数(在问题的最后发布).

任何建议 - 即使是小建议 - 将不胜感激.希望它能使我的代码更好,但至少它应该教会我一些东西.:)

概观

此函数将被调用至少10 ¹³次,并且可能经常被调用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的位置,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具有小的primes而q设置为0因此它们将被丢弃,其中lim的最佳值通过实验确定. )

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

一旦以这种方式恢复所有因子,使用由代码最佳解释的数学公式计算n为强伪荧光的碱基数1 <= b <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;
}

Answer 1

你似乎浪费了很多时间来分析这些因素.用除数的倒数替换除法要快得多(除法:〜15-80(!)周期,取决于除数,乘法:~4个周期),IF当然可以预先计算倒数.

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

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

编辑:

在咨询了Hacker's Delight(一本优秀的书,BTW)之后,通过利用代码中所有分区都是精确的(即余数为零),你似乎可以更快地完成它.

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

如何检查给定的x是否是d的精确倍数- 即x mod d = 0？简单!x mod d = 0iff 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并且有这么多非证人的整数要么是形式的两个素数的乘积

第(k + 1)*(2K + 1)

或者他们是Carmichael数字,有3个主要因素.我认为在上面的某些限制中,序列中的所有整数都具有这种形式,并且可以通过证明所有其他整数的证人的上限来验证这一点.例如,具有4个或更多因子的整数总是具有至多phi(n)/ 8个非证人.对于其他整数的基数,可以从您的公式中得出类似的结果.

至于微优化:每当你知道一个整数可以被一个商整除时,就有可能用一个乘以模2 ^ 64的商的倒数来代替除法.测试n%q == 0可以用测试代替

n*inverse_q <max_q,

其中inverse_q = q ^( - 1)mod 2 ^ 64和max_q = 2 ^ 64/q.显然,inverse_q和max_q需要预先计算,才能有效,但由于你使用筛子,我认为这不应该是一个障碍.