Eratosthenes 筛法 - 按位优化问题

Sieve of Eratosthenes - bitwise optimization problem

想必大家都遇到过埃拉托色尼筛按位运算的优化代码吧。我正在努力解决这个问题,但我对这个实现中的一个操作有疑问。这是 GeeksforGeeks 的代码:

bool ifnotPrime(int prime[], int x) { 
    // checking whether the value of element 
    // is set or not. Using prime[x/64], we find 
    // the slot in prime array. To find the bit 
    // number, we divide x by 2 and take its mod 
    // with 32. 
    return (prime[x / 64] & (1 << ((x >> 1) & 31))); 
} 

// Marks x composite in prime[] 
bool makeComposite(int prime[], int x) { 
    // Set a bit corresponding to given element. 
    // Using prime[x/64], we find the slot in prime  
    // array. To find the bit number, we divide x 
    // by 2 and take its mod with 32. 
    prime[x / 64] |= (1 << ((x >> 1) & 31)); 
} 

// Prints all prime numbers smaller than n. 
void bitWiseSieve(int n) { 
    // Assuming that n takes 32 bits, we reduce 
    // size to n/64 from n/2. 
    int prime[n / 64]; 

    // Initializing values to 0 . 
    memset(prime, 0, sizeof(prime)); 

    // 2 is the only even prime so we can ignore that 
    // loop starts from 3 as we have used in sieve of 
    // Eratosthenes . 
    for (int i = 3; i * i <= n; i += 2) { 

        // If i is prime, mark all its multiples as 
        // composite 
        if (!ifnotPrime(prime, i)) 
            for (int j = i * i, k = i << 1; j < n; j += k) 
                makeComposite(prime, j); 
    } 

    // writing 2 separately 
    printf("2 "); 

    // Printing other primes 
    for (int i = 3; i <= n; i += 2) 
        if (!ifnotPrime(prime, i)) 
            printf("%d ", i); 
} 

// Driver code 
int main() { 
    int n = 30; 
    bitWiseSieve(n); 
    return 0; 
} 

所以我的问题是:

  1. (prime[x/64] & (1 << ((x >> 1) & 31)) 的含义是什么,更具体地说 (1 << ((x >> 1) & 31));
  2. prime[x/64] 中,当我们使用 32 位整数时,为什么我们除以 64 而不是 32
  3. 如果n < 64int prime[n/64]是否正确?

1)x%32等价于x&31:逻辑与在最低5位。所以基本上 ((x>>1)&31) 意味着 ((x/2)%32)1<<x 表示 2^x 所以你要问的是 2^((x/2)%32).

2)实现中的一个优化是,它完全跳过了所有偶数。

3)n可以小于64

代码中存在多个问题:

  • makeComposite() 应该有 return 类型 void.
  • (1 << ((x >> 1) & 31)) 如果 x == 63 有未定义的行为,因为 1 << 31 溢出类型 int 的范围。您应该使用 1U 或最好使用 1UL 以确保 32 位。
  • 移动 table 条目并屏蔽最后一位会更简单:return (prime[x / 64] >> ((x >> 1) & 31)) & 1;
  • 不要假设类型 int 有 32 位,您应该使用 uint32_t 作为位数组的类型。
  • 如果 n 不是 64 的倍数,则数组 int prime[n / 64]; 太短。请改用 uint32_t prime[n / 64 + 1];。对于您的示例 n = 30,这是一个问题,因此创建的数组长度为 0
  • ifnotPrime(n) 只有 return 是奇数的有效结果。更改此函数并将其命名为 isOddPrime().
  • 可能会更好且更具可读性

关于您的问题:

what is the meaning of (prime[x/64] & (1 << ((x >> 1) & 31)) and more specifically (1 << ((x >> 1) & 31))?

x先除以2(右移一位),因为数组中只有奇数才有位,然后用31屏蔽结果,保留低5位作为位号在这个词中。对于任何无符号值 xx & 31 等价于 x % 32x / 64是测试这个位的字号。

如上所述,1 是一个 int,因此不应向左移动 31 个位置。使用 1UL 确保类型至少有 32 位并且可以移动 31 个位置。

in prime[x/64] why do we divide by 64 and not with 32, when we work with 32-bit integer?

数组中的位对应奇数,所以一个32位的字包含64个数的素数信息:32个已知是合数的偶数和32个奇数,如果数是复合的。

Is int prime[n/64] correct if n < 64?

不对,如果n不是64的倍数是不正确的:大小表达式应该是(n + 63) / 64,或者更好int prime[n/64 + 1]


这是修改后的版本,您可以在其中传递命令行参数:

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>

bool isOddPrime(const uint32_t prime[], unsigned x) {
    // checking whether the value of element
    // is set or not. Using prime[x/64], we find
    // the slot in prime array. To find the bit
    // number, we divide x by 2 and take its mod
    // with 32.
    return 1 ^ ((prime[x / 64] >> ((x >> 1) & 31)) & 1);
}

// Marks x composite in prime[]
void makeComposite(uint32_t prime[], unsigned x) {
    // Set a bit corresponding to given element.
    // Using prime[x/64], we find the slot in prime
    // array. To find the bit number, we divide x
    // by 2 and take its mod with 32.
    prime[x / 64] |= (1UL << ((x >> 1) & 31));
}

// Prints all prime numbers smaller than n.
void bitWiseSieve(unsigned n) {
    // Assuming that n takes 32 bits, we reduce
    // size to n/64 from n/2.
    uint32_t prime[n / 64 + 1];

    // Initializing values to 0 .
    memset(prime, 0, sizeof(prime));

    // 2 is the only even prime so we can ignore that
    // loop starts from 3 as we have used in sieve of
    // Eratosthenes .
    for (unsigned i = 3; i * i <= n; i += 2) {
        // If i is prime, mark all its multiples as composite
        if (isOddPrime(prime, i)) {
            for (unsigned j = i * i, k = i << 1; j < n; j += k)
                makeComposite(prime, j);
        }
    }

    // writing 2 separately
    if (n >= 2)
        printf("2\n");

    // Printing other primes
    for (unsigned i = 3; i <= n; i += 2) {
        if (isOddPrime(prime, i))
            printf("%u\n", i);
    }
}

// Driver code
int main(int argc, char *argv[]) {
    unsigned n = argc > 1 ? strtol(argv[1], NULL, 0) : 1000;
    bitWiseSieve(n);
    return 0;
}