一种在C中找到最接近无符号长整数（32位宽）的素数的方法？

Question

我正在寻找一种方法来找到最接近的素数。大于或小于，没关系，只是最接近的（最好没有溢出。）至于速度，如果它可以在 1GHz 机器上大约 50 毫秒内计算出来（在软件中，在 Linux 中运行），我会欣喜若狂。

Answer 1

在 (2^32 - 1) 范围内的 largest prime gap 是 (335)。在一次性设置之后，有 (6542) 个小于 (2^16) 的素数可以制成表格并用于筛选连续的奇数值。显然，只有素数

或者：deterministic variant of the Miller-Rabin test，带有 SPRP 基数：{2, 7, 61} 足以证明 32 位值的素数。由于测试的复杂性（需要取幂等），我怀疑对于这么小的候选人来说它是否会一样快。

编辑：实际上，如果乘法/减法可以保持为 32 位求幂（可能需要 64 位支持），则 M-R 测试可能会更好。主要间隙通常会小得多，从而使筛子设置成本过高。如果没有大型查找表等，您也可能会从更好的缓存局部性中获得提升。

此外：素数 {2, 3, 5, 7, 11, 13, 17, 19, 23} = (223092870) 的乘积。显式测试 [2, 23] 中的任何候选人。计算最大公约数：g = gcd(u, 223092870UL)。如果(g != 1)，候选人是复合的。如果(g == 1 && u < (29 * 29))，候选人（u > 23）肯定是素数。否则，请继续进行更昂贵的测试。使用 32 位算术的单个 gcd 测试非常便宜，根据 Mertens 的 (?) 定理，这将检测到 ~ 68.4% 的所有奇合数。

Answer 2

更新 2：修复（以严厉的方式）一些导致小 n 错误答案的错误。感谢 Brett Hale 的注意！还添加了一些断言来记录一些假设。

更新：我对此进行了编码，它似乎足够快以满足您的要求（在

它是用 C++ 编写的，因为那是我的筛子代码（我改编自）所在的，但转换为 C 应该很简单。内存使用量也（相对）小，您可以通过检查看到。

您可以看到，由于调用函数的方式，返回的数字是适合 32 位的最接近的素数，但实际上这是一回事，因为 2^32 周围的素数是 4294967291 和 4294967311。

我试图确保不会因为整数溢出而出现任何错误（因为我们正在处理高达 UINT_MAX 的数字）；希望我没有在那里犯错。如果您想使用 64 位类型（或者您知道您的数字将小于 2^32-256），则可以简化代码，因为您不必担心在循环条件中回绕。只要您愿意将小素数计算/存储到所需的限制，这个想法也适用于更大的数字。

我还应该注意，对于这些数字，小素数筛运行得非常快（粗略测量为 4-5 毫秒），所以如果你特别缺乏内存，每次都运行它而不是存储小素数是可行（在这种情况下，您可能希望使 mark[] 数组更节省空间）

#include <iostream>
#include <cmath>
#include <climits>
#include <cassert>

using namespace std;

typedef unsigned int UI;

const UI MAX_SM_PRIME = 1 << 16;
const UI MAX_N_SM_PRIMES = 7000;
const UI WINDOW = 256;

void getSMPrimes(UI primes[]) {
  UI pos = 0;
  primes[pos++] = 2;

  bool mark[MAX_SM_PRIME / 2] = {false};
  UI V_SM_LIM = UI(sqrt(MAX_SM_PRIME / 2));
  for (UI i = 0, p = 3; i < MAX_SM_PRIME / 2; ++i, p += 2)
    if (!mark[i]) {
      primes[pos++] = p;
      if (i < V_SM_LIM)
        for (UI j = p*i + p + i; j < MAX_SM_PRIME/2; j += p)
          mark[j] = true;
      }
  }

UI primeNear(UI n, UI min, UI max, const UI primes[]) {
  bool mark[2*WINDOW + 1] = {false};

  if (min == 0) mark[0] = true;
  if (min <= 1) mark[1-min] = true;

  assert(min <= n);
  assert(n <= max);
  assert(max-min <= 2*WINDOW);

  UI maxP = UI(sqrt(max));
  for (int i = 0; primes[i] <= maxP; ++i) {
    UI p = primes[i], k = min / p;
    if (k < p) k = p;
    UI mult = p*k;
    if (min <= mult)
      mark[mult-min] = true;
    while (mult <= max-p) {
      mult += p;
      mark[mult-min] = true;
      }
    }

  for (UI s = 0; (s <= n-min) || (s <= max-n); ++s)
    if ((s <= n-min) && !mark[n-s-min])
      return n-s;
    else if ((s <= max-n) && !mark[n+s-min])
      return n+s;

  return 0;
  }

int main() {
  UI primes[MAX_N_SM_PRIMES];
  getSMPrimes(primes);

  UI n;
  while (cin >> n) {
    UI win_min = (n >= WINDOW) ? (n-WINDOW) : 0;
    UI win_max = (n <= UINT_MAX-WINDOW) ? (n+WINDOW) : UINT_MAX;

    if (!win_min)
      win_max = 2*WINDOW;
    else if (win_max == UINT_MAX)
      win_min = win_max-2*WINDOW;

    UI p = primeNear(n, win_min, win_max, primes);
    cout << "found nearby prime " << p << " from window " << win_min << ' ' << win_max << '\n';
    }
  }

---

如果您知道最高 2^16 的素数，则可以在该范围内筛选区间（只有 6542

基本上，做一个常规的 Eratosthenes 筛选来获得“小”素数（比如前 7000 个）。显然你只需要在程序开始时做一次，但应该很快。

然后，假设您的“目标”数字是“a”，请考虑某个 n 值的区间 [a-n/2, a+n/2)。可能 n = 128 是一个合理的起点；如果第一个数字中的数字都是复合的，您可能需要尝试相邻的间隔。

对于每个“小”质数 p，在范围内划掉它的倍数，使用除法来查找从哪里开始。一种优化是您只需要从 p*p 开始交叉多个倍数（这意味着一旦 p*p 高于区间，您就可以停止考虑素数）。

除了前几个素数之外，大多数素数在区间内都有一个或零倍数；要利用这一点，您可以预先忽略前几个素数的倍数。最简单的方法是忽略所有偶数，但忽略 2、3 和 5 的倍数并不少见；这留下了与 1、7、11、13、17、19、23 和 29 mod 30 一致的整数（有 8 个，在筛选大范围时很好地映射到一个字节的位）。

...有点离题了；无论如何，一旦您处理了所有小素数（直到 p*p > a+n/2），您只需在区间中查找您没有划掉的数字；因为您希望最接近起点的人在那里寻找并在两个方向上向外搜索。