使第一个素数达到极限的最快方法
Fastest way to get the first prime numbers up to a limit
就CPU而言，使第一个素数达到极限的最快方法是什么？
在我的机器上计算前 1,000,000,000 个数字需要 4.8 秒。它甚至比从文件中读取它还要快！
public static void main(String[] args) {
long startingTime = System.nanoTime();
int limit = 1_000_000_000;
BitSet bitSet = findPrimes(limit);
System.out.println(2);
System.out.println(3);
System.out.println(5);
System.out.println(7);
limit = limit / 2 + 1;
for (int i = 6; i < limit; i++) {
if (!bitSet.get(i)) {
int p = i * 2 - 1;
//System.out.println(p);
}
}
System.out.println("Done in " + (System.nanoTime() - startingTime) / 1_000_000 + " milliseconds");
}
public static BitSet findPrimes(int limit) {
BitSet bitSet = new BitSet(limit / 2 + 1);
int size = (int) Math.sqrt(limit);
size += size % 2;
for (int prime = 3; prime < size; prime += 2) {
if (bitSet.get(prime / 2 + 1)) {
continue;
}
for (int i = prime / 2 + 1; i < bitSet.size(); i += prime) {
bitSet.set(i);
}
}
return bitSet;
}
计算质数的最有效方法是使用 Sieve of Eratosthenes，它最早发现于千年前的古希腊。
然而，问题在于它需要大量内存。当我刚刚声明一个 1,000,000,000 的布尔数组时，我的 Java 应用程序崩溃了。
所以我做了两个内存优化来让它工作
- 由于 2 之后的所有质数都是奇数，我通过将索引映射到
int p = i * 2 - 1; 来将数组大小减半
- 然后我没有使用布尔值数组，而是使用了
BitSet 对长整数数组进行位运算。
通过使用埃拉托色尼筛法进行的这两项优化，您可以在 4.8 秒内计算出 btSet。至于打印出来，那就另当别论了;)
有多种方法可以快速计算素数 -- 'fastest' 最终将取决于许多因素（硬件、算法、数据格式等的组合）。这已经在这里得到回答： Fastest way to list all primes below N 。然而，那是针对 Python，当我们考虑 C（或其他低级语言）
时，问题变得更有趣
您可以看到：
渐近地，阿特金筛法具有更好的计算复杂度（O(N) 或更低*——见注释），但实际上埃拉托色尼筛法 (O(N log log N)) 非常好，并且它的优点是非常容易实施：
/* Calculate 'isprime[i]' to tell whether 'i' is prime, for all 0 <= i < N */
void primes(int N, bool* isprime) {
int i, j;
/* first, set odd numbers to prime, evens to not-prime */
for (i = 0; i < N; ++i) isprime[i] = i % 2 == 1;
/* special cases for i<3 */
isprime[0] = isprime[1] = false;
/* Compute square root via newton's algorithm
* (this can be replaced by 'i * i <= N' in the foor loop)
*/
int sqrtN = (N >> 10) + 1024;
for (i = 0; i < 32; ++i) sqrtN = (N / sqrtN + sqrtN) / 2;
/* iterate through all odd numbers */
isprime[2] = true;
for (i = 3; i <= sqrtN /* or: i*i <= N */; i += 2) {
/* check if this odd number is still prime (i.e. if it has not been marked off) */
if (isprime[i]) {
/* mark off multiples of the prime */
j = 2 * i;
isprime[j] = false;
for (j += i; j < N; j += 2 * i) {
isprime[j] = false;
}
}
}
}
该代码非常清晰简洁，但需要大约 6 秒来计算 < 1_000_000_000（10 亿）的素数，并使用 N * sizeof(bool) == N == 1000000000 字节的内存。不太好
我们可以使用 bit-fiddling 和 uint64_t 类型（使用 #include <stdint.h>）来生成以相同复杂度运行的代码，但希望速度更快并且使用更少的内存。我们可以立即将内存减少到 N/8（每个布尔值 1 位而不是 1 字节）。此外，我们只能计算奇数并削减更多的内存——2 倍，使我们的最终使用量达到 N/16 字节。
代码如下：
/* Calculate that the (i//2)%64'th bit of 'isprime[(i//2)/64]' to tell whether 'i' is prime, for all 0 <= i < N,
* with 'i % 2 == 0'
*
* Since the array 'isprime' can be seen as a long bistring:
*
* isprime:
* [0] [1] ...
* +---------+---------+
* |01234....|01234....| ...
* +---------+---------+
*
* Where each block is a 64 bit integer. Therefore, we can map the odd natural numbers to this with the following formula:
* i -> index=((i / 2) / 64), bit=((i / 2) % 64)
*
* So, '1' is located at 'primes[0][0]' (where the second index represnts the bit),
* '3' is at 'primes[0][1]', 5 at 'primes[0][2]', and finally, 129 is at 'primes[1][0]'
*
* And so forth.
*
* Then, we use that value as a boolean to indicate whether the 'i' that maps to it is prime
*
*/
void primes_bs(int N, uint64_t* isprime) {
uint64_t i, j, b/* block */, c/* bit */, v, ii;
/* length of the array is roundup(N / 128) */
int len = (N + 127) / 128;
/* set all to isprime */
for (i = 0; i < len; ++i) isprime[i] = ~0;
/* set i==1 to not prime */
isprime[0] &= ~1ULL;
/* Compute square root via newton's algorithm */
int sqrtN = (N >> 10) + 1024;
for (i = 0; i < 32; ++i) sqrtN = (N / sqrtN + sqrtN) / 2;
/* Iterate through every word/chunk and handle its bits */
uint64_t chunk;
for (b = 0; b <= (sqrtN + 127) / 128; ++b) {
chunk = isprime[b];
c = 0;
while (chunk && c < 64) {
if (chunk & 1ULL) {
/* hot bit, so is prime */
i = 128 * b + 2 * c + 1;
/* iterate 3i, 5i, 7i, etc...
* BUT, j is the index, so is basically 3i/2, 5i/2, 7i/2,
* that way we don't need as many divisions per iteration,
* and we can use 'j += i' (since the index only needs 'i'
* added to it, whereas the value needs '2i')
*/
for (j = 3 * i / 2; j < N / 2; j += i) {
isprime[j / 64] &= ~(1ULL << (j % 64));
}
}
/* chunk will not modify itself */
c++;
chunk = isprime[b] >> c;
}
}
/* set extra bits to 0 -- unused */
if (N % 128 != 0) isprime[len - 1] &= (1ULL << ((N / 2) % 64)) - 1;
}
此代码在我的机器上用 ~2.7 秒计算所有小于 1000000000 (1_000_000_000) 的数字的素数，这是一个很大的改进！而且您使用的内存更少（其他方法使用的内存的 1/16。
但是，如果您在其他代码中使用它，界面会比较棘手——您可以使用一个函数来提供帮助：
/* Calculate if 'x' is prime from a dense bitset */
bool calc_is_prime(uint64_t* isprime, uint64_t x) {
/**/ if (x == 2) return true; /* must be included, since 'isprime' only encodes odd numbers */
else if (x % 2 == 0) return false;
return (isprime[(x / 2) / 64] >> ((x / 2) % 64))) & 1ULL;
}
您可以走得更远，将这种映射索引的方法推广到 odd-only 值 (https://en.wikipedia.org/wiki/Wheel_factorization)，但生成的代码可能会更慢——因为内部循环和操作是不再是 2 的幂，您要么必须进行 table 查找，要么需要一些非常聪明的技巧。您想要使用轮子分解的唯一情况是如果您的内存非常有限，在这种情况下 primorial(3) == 30 的轮子分解将允许您每 7 位存储 30 个数字（而不是每 7 位存储 2 个数字） 1 位）。或者，假设 primorial(5) == 2310，您可以为每 341
存储 2310 个数字
然而，大轮子分解将再次将一些指令从移位和 power-of-two 操作更改为 table 查找，对于大轮子，这实际上最终可能会慢得多
备注：
- 从技术上讲，存在 O(N / (log log N)) 的阿特金筛法版本，但您不再使用相同的密集数组（如果是，复杂度必须至少为 O(N ))
就CPU而言，使第一个素数达到极限的最快方法是什么？
在我的机器上计算前 1,000,000,000 个数字需要 4.8 秒。它甚至比从文件中读取它还要快！
public static void main(String[] args) {
long startingTime = System.nanoTime();
int limit = 1_000_000_000;
BitSet bitSet = findPrimes(limit);
System.out.println(2);
System.out.println(3);
System.out.println(5);
System.out.println(7);
limit = limit / 2 + 1;
for (int i = 6; i < limit; i++) {
if (!bitSet.get(i)) {
int p = i * 2 - 1;
//System.out.println(p);
}
}
System.out.println("Done in " + (System.nanoTime() - startingTime) / 1_000_000 + " milliseconds");
}
public static BitSet findPrimes(int limit) {
BitSet bitSet = new BitSet(limit / 2 + 1);
int size = (int) Math.sqrt(limit);
size += size % 2;
for (int prime = 3; prime < size; prime += 2) {
if (bitSet.get(prime / 2 + 1)) {
continue;
}
for (int i = prime / 2 + 1; i < bitSet.size(); i += prime) {
bitSet.set(i);
}
}
return bitSet;
}
计算质数的最有效方法是使用 Sieve of Eratosthenes，它最早发现于千年前的古希腊。
然而，问题在于它需要大量内存。当我刚刚声明一个 1,000,000,000 的布尔数组时，我的 Java 应用程序崩溃了。
所以我做了两个内存优化来让它工作
- 由于 2 之后的所有质数都是奇数，我通过将索引映射到
int p = i * 2 - 1; 来将数组大小减半
- 然后我没有使用布尔值数组，而是使用了
BitSet对长整数数组进行位运算。
通过使用埃拉托色尼筛法进行的这两项优化，您可以在 4.8 秒内计算出 btSet。至于打印出来，那就另当别论了;)
有多种方法可以快速计算素数 -- 'fastest' 最终将取决于许多因素（硬件、算法、数据格式等的组合）。这已经在这里得到回答： Fastest way to list all primes below N 。然而，那是针对 Python，当我们考虑 C（或其他低级语言）
时，问题变得更有趣您可以看到：
渐近地，阿特金筛法具有更好的计算复杂度（O(N) 或更低*——见注释），但实际上埃拉托色尼筛法 (O(N log log N)) 非常好，并且它的优点是非常容易实施：
/* Calculate 'isprime[i]' to tell whether 'i' is prime, for all 0 <= i < N */
void primes(int N, bool* isprime) {
int i, j;
/* first, set odd numbers to prime, evens to not-prime */
for (i = 0; i < N; ++i) isprime[i] = i % 2 == 1;
/* special cases for i<3 */
isprime[0] = isprime[1] = false;
/* Compute square root via newton's algorithm
* (this can be replaced by 'i * i <= N' in the foor loop)
*/
int sqrtN = (N >> 10) + 1024;
for (i = 0; i < 32; ++i) sqrtN = (N / sqrtN + sqrtN) / 2;
/* iterate through all odd numbers */
isprime[2] = true;
for (i = 3; i <= sqrtN /* or: i*i <= N */; i += 2) {
/* check if this odd number is still prime (i.e. if it has not been marked off) */
if (isprime[i]) {
/* mark off multiples of the prime */
j = 2 * i;
isprime[j] = false;
for (j += i; j < N; j += 2 * i) {
isprime[j] = false;
}
}
}
}
该代码非常清晰简洁，但需要大约 6 秒来计算 < 1_000_000_000（10 亿）的素数，并使用 N * sizeof(bool) == N == 1000000000 字节的内存。不太好
我们可以使用 bit-fiddling 和 uint64_t 类型（使用 #include <stdint.h>）来生成以相同复杂度运行的代码，但希望速度更快并且使用更少的内存。我们可以立即将内存减少到 N/8（每个布尔值 1 位而不是 1 字节）。此外，我们只能计算奇数并削减更多的内存——2 倍，使我们的最终使用量达到 N/16 字节。
代码如下：
/* Calculate that the (i//2)%64'th bit of 'isprime[(i//2)/64]' to tell whether 'i' is prime, for all 0 <= i < N,
* with 'i % 2 == 0'
*
* Since the array 'isprime' can be seen as a long bistring:
*
* isprime:
* [0] [1] ...
* +---------+---------+
* |01234....|01234....| ...
* +---------+---------+
*
* Where each block is a 64 bit integer. Therefore, we can map the odd natural numbers to this with the following formula:
* i -> index=((i / 2) / 64), bit=((i / 2) % 64)
*
* So, '1' is located at 'primes[0][0]' (where the second index represnts the bit),
* '3' is at 'primes[0][1]', 5 at 'primes[0][2]', and finally, 129 is at 'primes[1][0]'
*
* And so forth.
*
* Then, we use that value as a boolean to indicate whether the 'i' that maps to it is prime
*
*/
void primes_bs(int N, uint64_t* isprime) {
uint64_t i, j, b/* block */, c/* bit */, v, ii;
/* length of the array is roundup(N / 128) */
int len = (N + 127) / 128;
/* set all to isprime */
for (i = 0; i < len; ++i) isprime[i] = ~0;
/* set i==1 to not prime */
isprime[0] &= ~1ULL;
/* Compute square root via newton's algorithm */
int sqrtN = (N >> 10) + 1024;
for (i = 0; i < 32; ++i) sqrtN = (N / sqrtN + sqrtN) / 2;
/* Iterate through every word/chunk and handle its bits */
uint64_t chunk;
for (b = 0; b <= (sqrtN + 127) / 128; ++b) {
chunk = isprime[b];
c = 0;
while (chunk && c < 64) {
if (chunk & 1ULL) {
/* hot bit, so is prime */
i = 128 * b + 2 * c + 1;
/* iterate 3i, 5i, 7i, etc...
* BUT, j is the index, so is basically 3i/2, 5i/2, 7i/2,
* that way we don't need as many divisions per iteration,
* and we can use 'j += i' (since the index only needs 'i'
* added to it, whereas the value needs '2i')
*/
for (j = 3 * i / 2; j < N / 2; j += i) {
isprime[j / 64] &= ~(1ULL << (j % 64));
}
}
/* chunk will not modify itself */
c++;
chunk = isprime[b] >> c;
}
}
/* set extra bits to 0 -- unused */
if (N % 128 != 0) isprime[len - 1] &= (1ULL << ((N / 2) % 64)) - 1;
}
此代码在我的机器上用 ~2.7 秒计算所有小于 1000000000 (1_000_000_000) 的数字的素数，这是一个很大的改进！而且您使用的内存更少（其他方法使用的内存的 1/16。
但是，如果您在其他代码中使用它，界面会比较棘手——您可以使用一个函数来提供帮助：
/* Calculate if 'x' is prime from a dense bitset */
bool calc_is_prime(uint64_t* isprime, uint64_t x) {
/**/ if (x == 2) return true; /* must be included, since 'isprime' only encodes odd numbers */
else if (x % 2 == 0) return false;
return (isprime[(x / 2) / 64] >> ((x / 2) % 64))) & 1ULL;
}
您可以走得更远，将这种映射索引的方法推广到 odd-only 值 (https://en.wikipedia.org/wiki/Wheel_factorization)，但生成的代码可能会更慢——因为内部循环和操作是不再是 2 的幂，您要么必须进行 table 查找，要么需要一些非常聪明的技巧。您想要使用轮子分解的唯一情况是如果您的内存非常有限，在这种情况下 primorial(3) == 30 的轮子分解将允许您每 7 位存储 30 个数字（而不是每 7 位存储 2 个数字） 1 位）。或者，假设 primorial(5) == 2310，您可以为每 341
然而，大轮子分解将再次将一些指令从移位和 power-of-two 操作更改为 table 查找，对于大轮子，这实际上最终可能会慢得多
备注：
- 从技术上讲，存在 O(N / (log log N)) 的阿特金筛法版本，但您不再使用相同的密集数组（如果是，复杂度必须至少为 O(N ))