【发布时间】：2022-01-17 08:00:25
【问题描述】：
我正在尝试在 C 中实现一个高效的分段素数筛。它基本上是一个 Eratosthenes 筛，但每个段都被分割成可以很好地放入缓存中的大小。
在我的版本中，有一个标志位数组，其中每个位都是一个连续的奇数。当每个位是已知素数的倍数时，通过用AND 屏蔽来擦除每个位。
这部分代码消耗了大约 90% 的运行时间。每一个脏点的代码都有我在cmets中解释过的原因，但是整体操作非常简单。
- 获取质数。
- 计算它的平方和它比缓存块起始点所代表的数略大的倍数。
- 拿大一点的。
- 擦除该位，将基本素数添加到自身两次，然后重复直到缓存块的末尾。
就是这样。
有一个名为primesieve 的程序可以非常快速地执行此操作。它比我的版本快大约 3 倍。我阅读了它关于算法及其代码的文档，并在我的代码中应用了任何合理的东西。
由于有一个已知程序比我的要快得多，我将进一步调查他们在做什么以及我没有做什么，但在此之前，我发布了这个问题以获得额外的帮助，如果你能帮助我找出答案哪个部分没有高效运行。
再说一遍，这个单一的例程消耗了 90% 的运行时间，所以我真正专注于让这部分运行得更快。
这是旧版本，我在发帖后做了一些修改，那个在这个下面。 cmets 仍然适用。
#include <stdlib.h>
#include <stdio.h>
//size of the cache block (64K)
#define C 0x10000
static unsigned sq(unsigned a) {
return a * a;
}
//`f` is the array of flags. `st` is the starting point being a multiple of
//`C * 16`. Each byte can hold 16 odd numbers, thus having `C` multiplied by 16.
//`p` is the array of prime numbers from 67 up to `sqrt(st + C * 16)`. Primes
//below 67 are handled specially by other routines.
__attribute__((noinline)) //Don't inline for testing.
static void sieve_bit(unsigned *f, unsigned st, unsigned *p) {
//Doing table access than computing the mask in runtime was about 15% faster
//on my machine. I'm guessing that probably the bit-shifts on x86 only being
//able to use the CL register as a variable counter creates a dependency
//when there is a series of shifts, making it slower than multiple L1 cache
//accesses.
static const unsigned m[] = {
~(1u << 0), ~(1u << 1), ~(1u << 2), ~(1u << 3), ~(1u << 4),
~(1u << 5), ~(1u << 6), ~(1u << 7), ~(1u << 8), ~(1u << 9),
~(1u << 10), ~(1u << 11), ~(1u << 12), ~(1u << 13), ~(1u << 14),
~(1u << 15), ~(1u << 16), ~(1u << 17), ~(1u << 18), ~(1u << 19),
~(1u << 20), ~(1u << 21), ~(1u << 22), ~(1u << 23), ~(1u << 24),
~(1u << 25), ~(1u << 26), ~(1u << 27), ~(1u << 28), ~(1u << 29),
~(1u << 30), ~(1u << 31)
};
unsigned p2 = sq(*p);
do {
unsigned r = st % *p;
//This calculates the starting point. The result of `n` will be an odd
//number multiple of `*p` a bit larger than `st`.
unsigned n = st + *p - r + (-(r & 1) & *p);
//If the square of `*p` is larger than `n`, use it instead.
if (p2 > n) {
n = p2;
}
//The loop is unrolled 8 times, which gave the fastest result on my
//machine.
for (;; n += *p * 16) {
//Jumps to the next stage when there are no space left for 8
//operations to be done. This could simply be a `break` followed by
//a loop, but this was slightly (1~3%) faster than the simpler
//alternative.
int d = st + C * 16 - n;
//as mentioned in the comments, the value of `d` relies on
//implementation dependent behaviour for 2's complement machines.
//`n` can be larger than `st + C * 16`, in which case the wrapped
//unsigned value is converted to a negative `int`.
if (d <= (int)*p * 14) {
if (d <= (int)*p * 6) {
if (d <= (int)*p * 2) {
if (d <= (int)*p * 0) goto L0; else goto L1;
} else {
if (d <= (int)*p * 4) goto L2; else goto L3;
}
} else {
if (d <= (int)*p * 10) {
if (d <= (int)*p * 8) goto L4; else goto L5;
} else {
if (d <= (int)*p * 12) goto L6; else goto L7;
}
}
}
//The multiples of primes are erased.
#define ur(i) do {\
unsigned _n = (n - st) / 2 + *p * i;\
f[_n / 32] &= m[_n % 32];\
} while (0)
ur(0); ur(1); ur(2); ur(3); ur(4); ur(5); ur(6); ur(7);
}
//Erase with the leftovers.
L7: ur(6);
L6: ur(5);
L5: ur(4);
L4: ur(3);
L3: ur(2);
L2: ur(1);
L1: ur(0);
L0:
p2 = sq(*++p);
#undef ur
} while (p2 < st + C * 16);
}
//This could break if `TscInvariant` CPUID is false, which is probably rare on
//modern machines?
static inline unsigned long long rdtscp() {
unsigned _;
return __builtin_ia32_rdtscp(&_);
}
//This isn't used in actuall code. Just a simple one for a test. Fill with prime
//numbers enough to sieve all 32-bit primes.
static inline void fillPrimes(unsigned *p) {
for (unsigned n = 67, c = 2; n <= 65521; n += c ^= 6) {
for (unsigned d = 5, c = 4; d * d <= n; d += c ^= 6) {
if (!(n % d)) goto next;
}
*p++ = n;
next:;
}
}
int main() {
unsigned p[8000];
fillPrimes(p);
unsigned f[C / sizeof(unsigned)];
puts("start sieve");
unsigned long long c = rdtscp();
for (int i = 0; i < 2000; ++i) {
volatile unsigned *vf = f;
sieve_bit((unsigned *)vf, C * 16 * i, p);
}
c = rdtscp() - c;
printf("%llu\n", c);
return 0;
}
稍作休息后，我可以在循环中找到一些冗余计算，将它们取出来获得大约 6~7% 的速度。
static void sieve_bit(unsigned *f, unsigned st, unsigned *p) {
static const unsigned m[] = {
~(1u << 0), ~(1u << 1), ~(1u << 2), ~(1u << 3), ~(1u << 4),
~(1u << 5), ~(1u << 6), ~(1u << 7), ~(1u << 8), ~(1u << 9),
~(1u << 10), ~(1u << 11), ~(1u << 12), ~(1u << 13), ~(1u << 14),
~(1u << 15), ~(1u << 16), ~(1u << 17), ~(1u << 18), ~(1u << 19),
~(1u << 20), ~(1u << 21), ~(1u << 22), ~(1u << 23), ~(1u << 24),
~(1u << 25), ~(1u << 26), ~(1u << 27), ~(1u << 28), ~(1u << 29),
~(1u << 30), ~(1u << 31)
};
unsigned p2 = sq(*p);
do {
unsigned n =
(p2 > st + *p * 2 ? p2 - st : *p - st % *p + (-(st % *p & 1) & *p)) / 2;
for (;; n += *p * 8) {
int d = C * 8 - (int)n;
if (d <= (int)*p * 7) {
if (d <= (int)*p * 3) {
if (d <= (int)*p * 1) {
if (d <= (int)*p * 0) goto L0; else goto L1;
} else {
if (d <= (int)*p * 2) goto L2; else goto L3;
}
} else {
if (d <= (int)*p * 5) {
if (d <= (int)*p * 4) goto L4; else goto L5;
} else {
if (d <= (int)*p * 6) goto L6; else goto L7;
}
}
}
#define ur(i) f[(n + *p * i) / 32] &= m[(n + *p * i) % 32]
ur(0); ur(1); ur(2); ur(3); ur(4); ur(5); ur(6); ur(7);
}
L7: ur(6);
L6: ur(5);
L5: ur(4);
L4: ur(3);
L3: ur(2);
L2: ur(1);
L1: ur(0);
L0:
p2 = sq(*++p);
#undef ur
} while (p2 < st + C * 16);
}
【问题讨论】：
-
尽量用左/右移位替换乘法/除法来优化。比如
(a * 16)可以写成(a << 4)。 -
@kiner_shah 我厌倦了阅读汇编输出。编译器显然是这样做的。
-
可以在codereview问吗？
-
@kiner_shah 编译器无论如何都会进行这种优化，即使没有启用优化，甚至 TCC（不完全以优化闻名）也会这样做：godbolt.org/z/eqhxGcP38
-
@pmg，这个问题需要解决才能适合Code Review。您应该将提问者指向A guide to Code Review for Stack Overflow users，因为那里有些事情的处理方式不同 - 例如。我们需要很好地描述代码的用途以提供上下文，并且问题标题应该简单地说明代码做了什么（问题总是，“如何我可以改进吗？”）。代码正确运行很重要。如果可能，包括单元测试。
标签： c optimization primes