使用#pragma omp 进行硬并行化以找到第 N 个素数

Question

这个问题的目标是能够获得第 2.000.000 个素数，并能够分辨出第 2.000.000 个素数是哪个。

我们从这段代码开始：

#include <stdlib.h>
#include <stdio.h>

#define N 2000000

int p[N];

main(int na,char* arg[])
{
int i;
int pp,num;

printf("Number of primes to find: %d\n",N);

p[0] = 2;
p[1] = 3;
pp = 2;
num = 5;

while (pp < N)
{
  for (i=1; p[i]*p[i] <= num ;i++)
    if (num % p[i] == 0) break;
  if (p[i]*p[i] > num) p[pp++]=num;
  num += 2;
}

printf("The %d prime is: %d\n",N,p[N-1]);
exit(0);
}

现在我们被要求通过 pragma omp 使这个进程线程化。这是我到目前为止所做的：

#include <stdlib.h>
#include <stdio.h>

#define N 2000000
#define D 1415

int p[N];
main(int na,char* arg[])
{
int i,j;
int pp,num;

printf("Number of primes to find: %d\n",N);

p[0] = 2;
p[1] = 3;

pp = 2;
num = 5;

while (pp < D)
{
    for (i=1; p[i]*p[i] <= num ;i++)
        if (num % p[i] == 0) break;
    if (p[i]*p[i] > num) p[pp++]=num;
    num += 2;
}

int success = 0;
int t_num;
int temp_num = num;
int total = pp;

#pragma omp parallel num_threads(4) private(j, t_num, num, success)
{
    t_num = omp_get_thread_num();
    num = temp_num + t_num*2;

    #pragma omp for ordered schedule(static,4)
    for(pp=D; pp<N; pp++) {
        success = 0;
        while(success==0) {
            for (i=1; p[i]*p[i] <= num;i++) {
                if (num % p[i] == 0) break;
            }
            if (p[i]*p[i] > num) {
                p[pp] = num;
                success=1;
            }
            num+=8;
        }

    }
}

//sort(p, 0, N);

printf("El %d primer es: %d\n",N,p[N-1]);

exit(0);
}

现在让我解释一下我的“部分”解决方案，以及我的问题。

第一个 D 素数是通过序列码获得的，所以现在我可以检查大量数字的可分性。

每个线程运行一个素数对角线，因此线程之间没有依赖关系，也不需要同步。但是，这种方法存在以下问题：

1. 一个线程可能比另一个线程生成更多的素数
2. 作为问题 1. 的直接结果，它将生成 N 个素数，但它们不会被排序，因此当素数计数器 'pp' 达到 'N' 时，最后一个素数不是第 2.000.000 个素数，而是更高级的素数。
3. 也可能是当它生成 2.000.000 个素数时，能够到达真正的第 2.000.000 个素数的线程可能没有足够的时间将其放在素数数组“p”上。

问题/困境是：

我如何知道第 2.000.000 个素数何时生成？

提示： 有人告诉我，我应该分批（比方说）10.000 个素数候选者。然后当我不知道的事情发生时，我会知道最后一批 10.000 个候选者包含第 2.000.000 个素数，我可以使用快速排序对其进行排序。

我希望我说清楚了，这真的很艰难，我只是不停地尝试了几天。

Answer 1

如果您只需要 2000000 个素数，您可以维护一个 ~4.1MB 大小的位数组，并为每个找到的素数翻转位。不需要排序。通过实施仅赔率表示方案将您的位数组大小减半。

在段中使用Sieve of Eratosthenes，其大小与sqrt(top_value_of_range) 成比例（或类似的东西 - 目标是在每个段上执行大致相同的工作量）。对于 n=2000000、n*(log n + log(log n)) == 34366806 和 prime[771]^2 == 34421689（从 0 开始），预计算前 771 个奇数素数。

每个工人也可以计数，因为它会翻转位，所以当它们全部完成时，您将知道每个范围的计数，并且只需要扫描包含第 2 百万个素数的一个范围，在结束，找到那个素数。或者让每个工作人员根据其范围维护自己的位数组 - 您只需保留一个，并且可以丢弃其他的。

计算埃拉托色尼筛的伪代码为：

Input: an integer n > 1

Let A be an array of bool values, indexed by integers 3, 5, ... upto n,
initially all set to true.

count := floor( (n-1)/2 )
for i = 3, 5, 7, ..., while i^2 ≤ n:
  if A[i] is true:
    for j = i^2, i^2 + 2i, i^2 + 4i, ..., while j ≤ n:
      if A[j] is true:
        A[j]  := false
        count := count - 1

Now all 'i's such that A[i] is true are prime,
and 'count' is the total count of odd primes found.

Answer 2

我可以想到两种方法。

1. 一旦您获得了第 2 百万个素数的候选，您的线程将继续计算低于您的候选的素数，直到您没有丢失任何素数。然后你可以对素数列表进行排序并从中取出百万分之一。

2. 如果您的线程正在生成连续素数块，它们应该单独维护这些块，然后可以随后将素数块重新组合到主列表中。一旦找到第 2 百万个素数，执行重组的线程就可以终止程序。