优化我!(C,性能)-后续的位旋转问题

你似乎在浪费很多时间按因素划分。用除数倒数的乘法(除法:~15-80)代替除法要快得多( ! )周期,取决于除数,乘法:~4个周期),

而这似乎不太可能与 q , , -由于这些变量的范围,这是非常容易做到的 它总是来自于小的,静态的 pr[] 数组。预先计算这些素数的倒数,并将它们存储在另一个数组中。然后,不是除以 ,乘以第二个数组的倒数。(或创建单个结构数组。)

现在,用这种方法得到精确的除法结果需要一些技巧来补偿舍入误差。你会发现这种技术的血淋淋的细节 this document ,第138页。

咨询后 高效程序的奥秘 (一本好书,顺便说一句)在这个问题上,通过利用代码中所有除法都是精确的(即余数为零)这一事实,似乎可以使它变得更快。

这是奇怪的和基本的 = 2 ^字号 存在唯一的乘法逆 满足以下条件: d⃰ < B 和 d·d⃰ ≡ 1 (mod B) . 每 十 这是 ,这意味着 x/d ≡ x·d⃰ (mod B) . 这意味着你可以简单地用乘法代替除法,不需要额外的修正,检查,舍入问题等等。(这些定理的证明可以在书中找到。) 注意 这个乘法逆 等于前面方法定义的倒数!

如何检查给定的 是 -即。 x mod d = 0 ? 容易的! x模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 与任何其他无符号类型-只要使用相同的类型一致。

证据, 为什么? s和 怎样 这本书里都有。衷心推荐阅读:-)。

如果编译器支持GCC函数属性,则可以使用此属性标记纯函数:

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

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

你打开代码有什么原因吗 vals() 而不是使用 __builtin_ctz() ?

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

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

至于微优化:只要你知道一个整数可以被某个商整除,那么就可以用2^64的商的倒数来代替乘法的除法。测试n%q==0可以替换为一个测试

n*反向\u q<最大\u q,

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() 以及其他你使用的函数定义

尝试替换此模式(也适用于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的选项很少,或者某些选项更为频繁,那么您可以通过专门化这些值的函数来获得一些好处。

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

-fprofile-generate 选项,然后在具有代表性的数据集上运行程序(例如,一天的计算值)。

-fprofile-use 而不是选择。

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

2) 在多核机器上尝试这个很好,但是这个循环是单线程的。我猜有一个伞形函数可以将负载分散到几个线程上,从而使用更多的内核?

5) 你们正在做大量的模运算和除法运算。在C中,这是两个独立的命令(首先是“/”,然后是“%”)。但是在汇编中,这是一个命令:“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. while 循环。最好的方法是检查组件。

2. star( pr[__builtin_ffsll(bit)] ) 和 vals( pr[__builtin_ffsll(bit)] - 1) ? 这将用一些简单的工作来换取数组查找,但是如果表足够小的话,这可能是值得的。

3. f 直到你真的需要它(接近尾声,在你早出晚归之后)。你可以用类似这样的代码替换你的代码

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;

希望有帮助。

首先是吹毛求疵;-)你应该对你使用的类型更加小心。在某些地方,您似乎假设ulong是64位宽的,请使用 uint64_t

uint32_t )也许比大的更有效率。尤其是 第32页 有一条汇编指令可以一次性进行除法和模运算,叫做 divl

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

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

您是否尝试了第一个while循环的查表版本?你可以分开 smallprimes 在4个16位值中,查找它们的贡献并合并它们。但也许你需要副作用。

你有没有试着传入一个素数数组而不是把它们分开 smallprimes , q , r 和 s ? 因为我不知道外部代码是做什么的,我可能错了,但是有可能你也有一个函数可以把一些素数转换成 位图,在这个函数中,您可以有效地将位图转换回素数数组。此外,您似乎对 小素数 , q r s . 它应该可以为你节省一点点的处理每一个电话。

而且,你似乎知道传递的素数是分开的 n . 你对除n以外的每个素数的幂有足够的了解吗?如果可以通过将该信息传递给此函数来消除模运算,则可以节省大量时间。换句话说,如果 是 pow(p_0,e_0)*pow(p_1,e_1)*...*pow(p_k,e_k)*n_leftover e_i 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);

你要传递n的完全分解,所以你要分解连续的整数,然后在这里使用分解的结果。在我看来,你可能会受益于这样做的时候,找到一些因素。

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

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

#include 
#define SIZE (1024*1024) //must be 2^n
#define MASK (SIZE-1)
typedef struct {
    int p;
    int next;
} p_type;

p_type primes[SIZE];
int sieve[SIZE];

void init_sieve()
{
    int i,n;
    int count = 1;
    primes[1].p = 3;
    sieve[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;
                break;
            }
        }
        if (flag==0)
        {
            count++;
            primes[count].p = n;
            sieve[n>>1] = count;
        }
    }
}

int main()
{
    int ptr,n;
    init_sieve();
    printf("init_done\n");

    // factor odd numbers starting with 3
    for (n=1;1000000000>n;n++)
    {
        ptr = sieve[n&MASK];
        if (ptr == 0) //prime
        {
//            printf("%d is prime",n*2+1);
        }
        else //composite
        {
//            printf ("%d has divisors:",n*2+1);
            while(ptr!=0)
            {
//                printf ("%d ",primes[ptr].p);
                sieve[n&MASK]=primes[ptr].next;
                //move the prime to the next number it divides
                primes[ptr].next = sieve[(n+primes[ptr].p)&MASK];
                sieve[(n+primes[ptr].p)&MASK] = ptr;
                ptr = sieve[n&MASK];
            }
        }
//        printf("\n");
    }
    return 0;
}

其思想是为筛选中的每个条目维护一个链表。数字是通过简单地从链表中取出它们的系数来计算的。当它们被拉出时,它们被插入到下一个数字的列表中,下一个数字将被那个素数整除。这也是非常友好的缓存。筛子大小必须大于因子基中的最大素数。这样,这个筛子可以在7小时内达到2**40,这似乎是你的目标(除了需要64位的n)。

希望有帮助。

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

如果你要马上离开 (!smallprimes&!q) 为什么不在调用函数之前进行测试,这样可以节省函数调用的开销呢?

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

在没有分支和goto的情况下,将它们实际创建为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);

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

如果您使用的除数在编译时是未知的,但在运行时经常使用(多次除以相同的数),那么我建议使用 libdivide 库,它基本上在运行时实现编译器对编译时常量所做的优化(使用移位掩码等)。这可以带来巨大的好处。同时避免使用x%y==0表示z=x/y,z*y==x作为 ergosys 上述建议也应该有明显的改善。

你上面的代码是优化版吗?如果是的话,仍然有太多的除法操作大大消耗了CPU周期。

这个代码天生就有点执行过度

if (!smallprimes & !q)
    return 0;

更改为逻辑与;

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

会使它短路更快而不影响q值

以及以下代码

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

用于查找smallprimes的最后一组位。你为什么不用更简单的方法呢

ulong p = pr[__builtin_ctz(smallprimes)];

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