【问题标题】:Find position of prime number查找素数的位置
【发布时间】:2012-12-17 02:38:19
【问题描述】:

我需要找到第 N 个素数的相反操作,即给定一个素数,我需要找到它在中的位置

2, 3, 5, 7...

质数可以很大,按10^7 的顺序排列。而且,还有很多。

我有一个可以二进制搜索的预先计算的素数索引,但我也有 50k 的空间限制!可以过筛吗?还是有其他快速的方法?

编辑: 非常感谢所有精彩的答案,我没想到他们!我希望它们对寻求相同的其他人有用。

【问题讨论】:

  • 首先,您是否计算出可能的答案集中有多少个?有几种方法可以解决这个问题。空间/速度权衡、压缩和数学技巧(如孪生素数)。
  • 当我再次检查时,有 664570+ 个质数小于 10^7。预先计算的素数表上的二进制搜索不是这里的选项。
  • stackoverflow.com/questions/3918968/… 上查看我的答案,了解有关存储素数表和快速筛分的想法。您应该能够获得 one 每个字节的两个素数(请参阅en.wikipedia.org/wiki/Prime_gap)。二分查找是不必要的,您可以根据 P(n) ≈ n log n 和从那里线性查找来估计从哪里开始查找表。
  • 我认为筛分在这里不会很实用.. 你可以看看 Darrick Henry Lehmer 的素数计数方法:en.wikipedia.org/wiki/Prime-counting_function 他能够在 59 年的 IBM 701 上做 10^10
  • 非常感谢大家。 @ColonelPanic 50k 我实际上是指源文件限制,它是在线判断问题的子问题!

标签: c algorithm math primes sieve


【解决方案1】:

如果您先验地知道输入是素数,您可以使用近似值 pi(n) ≈ n / log n 和一个小的修正表来计算四舍五入结果不足以得到的素数正确的值 n。除了缓慢的蛮力方法之外,我认为这是您在大小限制内的最佳选择。

【讨论】:

  • 该表的大小非常关键。素数定理主要适用于大数;它在这里可能不太适用,但只有尝试才会知道。
  • 我已经有一段时间没有玩质数定理的数字了,但我认为值得一试。
  • 我认为他正在寻找一个确切的数字,并且有已知的计算这个数字的好方法(Derrick Henry Lehmer 的方法对我来说是最好的,Legendre 也有一个看起来可能在这些中可行的方法约束)
  • 显然他想要确切的数字。我的观点是,您可以使用基于反转素数定理估计并调整结果函数的近似值,然后针对错误的输入添加修正表。
  • 这行不通。素数的分布太糟糕了。素数定理告诉你 pi(n) 对 n/log(n) 是渐近的。这意味着它们之间的乘法差异趋于零。您需要一些在 pi(n) 的附加常数内的东西,以便修复方法有希望工作。 (并且 n/log(n) 不在 pi(n) 的加法常数内。)
【解决方案2】:

你的范围只有一千万,这对于这种东西来说是很小的。我有两个建议:

1) 以方便的间隔创建一个 pi(n) 表,然后使用分段的 Eratosthenes 筛来计算包含所需值的两个表条目之间的素数。区间的大小决定了所需表的大小和计算结果的速度。

2) 使用勒让德的 phi(x,a) 函数和 Lehmer 的素数计数公式直接计算结果。 phi 函数需要一些存储空间,我不确定具体需要多少。

考虑到您的问题规模,我可能会选择第一个替代方案。我的博客上提供了segmented Sieve of EratosthenesLehmer's 素数计数功能的实现。

编辑 1:

经过反思,我有第三种选择:

3) 使用对数积分来估计 pi(n)。它是单调递增的,并且在您需要的时间间隔内始终大于 pi(n)。但差异很小,永远不会超过 200。因此,您可以预先计算所有小于 1000 万的值的差异,制作 200 个变化点的表格,然后在需要时计算对数积分并在桌子。或者你可以用黎曼的 R 函数做类似的事情。

第三种选择占用的空间最少,但我怀疑第一种选择所需的空间不会太大,而且筛分可能比计算对数积分更快。所以我会坚持我原来的建议。 my blog 处有对数积分和黎曼 R 函数的实现。

编辑 2:

正如 cmets 所指出的那样,这并不是很好。请忽略我的第三条建议。

为了弥补我在提出一个不起作用的解决方案时的错误,我编写了一个程序,该程序使用 pi(n) 值表和 Eratosthenes 的分段筛来计算 n 的 pi(n) 值10000000。我将使用 Python,而不是原始海报要求的 C,因为 Python 更简单,更易于阅读。

我们从计算小于一千万平方根的筛选质数开始;这些素数将用于构建 pi(n) 的值表和执行计算最终答案的筛子。一千万的平方根是 3162.3。我们不想使用 2 作为筛选素数——我们将只筛选奇数,并将 2 视为一种特殊情况——但我们确实希望下一个素数大于平方根,因此列表筛选素数永远不会耗尽(这会导致错误)。所以我们使用这个非常简单的埃拉托色尼筛法来计算筛分素数:

def primes(n):
    b, p, ps = [True] * (n+1), 2, []
    for p in xrange(2, n+1):
        if b[p]:
            ps.append(p)
            for i in xrange(p, n+1, p):
                b[i] = False
    return ps

埃拉托色尼筛分两部分。首先,列出小于目标数的数字,从 2 开始。然后,从第一个未划线的数字开始,重复遍历该列表,并从列表中划掉该数字的所有倍数。最初,2 是第一个未划线的数字,因此划掉 4、6、8、10 等。然后 3 是下一个未划线的数字,因此划掉 6、9、12、15 等。然后 4 作为 2 的倍数被划掉,下一个未划掉的数字是 5,所以划掉 10、15、20、25,以此类推。继续,直到所有未交叉的数字都被考虑在内;未交叉的数字是素数。 p 上的循环依次考虑每个数,如果未交叉,则 i 上的循环将多个数划掉。

primes 函数返回一个包含 447 个素数的列表:2, 3, 5, 7, 11, 13, ..., 3121, 3137, 3163。我们从列表中删除 2 并将 446 个筛选素数存储在全局 ps 变量:

ps = primes(3163)[1:]

我们需要的主要函数计算一个范围内的素数。它使用一个筛子,我们将把它存储在一个全局数组中,这样它就可以被重用,而不是在每次调用 count 函数时重新分配:

sieve = [True] * 500

count 函数使用 Eratosthenes 的分段筛来计算从 lo 到 hi 范围内的素数(lo 和 hi 都包含在该范围内)。该函数有四个for 循环:第一个清除筛子,最后一个计算素数,另外两个执行筛分,其方式类似于上面显示的简单筛子:

def count(lo, hi):
    for i in xrange(500):
        sieve[i] = True
    for p in ps:
        if p*p > hi: break
        q = (lo + p + 1) / -2 % p
        if lo+q+q+1 < p*p: q += p
        for j in xrange(q, 500, p):
            sieve[j] = False
    k = 0
    for i in xrange((hi - lo) // 2):
        if sieve[i]: k += 1
    return k

函数的核心是循环for p in ps,它执行筛选,依次获取每个筛选质数 p。当筛选素数的平方大于范围的限制时,循环终止,因为所有素数都将在该点被识别(我们需要下一个大于平方根的素数的原因是为了存在筛选素数停止循环)。神秘变量 q 是 p 在 lo 到 hi 范围内的最小倍数入筛的偏移量(注意不是 p 在范围内的最小倍数,而是 p 在范围内的最小倍数的偏移量的索引范围,这可能会造成混淆)。 if 语句在引用一个完全平方数时增加 q。然后 j 上的循环从筛子中击出 p 的倍数。

我们以两种方式使用count 函数。第一次使用建立一个 pi(n) 值的表,该表是 1000 的倍数;第二种用途在表内插值。我们将表存储在一个全局变量 piTable 中:

piTable = [0] * 10000

我们根据原始请求选择参数 1000 和 10000 以将内存使用量控制在 50 KB 以内。 (是的,我知道最初的发布者放宽了这个要求。但我们无论如何都可以兑现它。)一万个 32 位整数将占用 40,000 字节的存储空间,从 lo 到 hi 的 1000 范围内筛选只需要 500 个字节存储空间,速度非常快。您可能想尝试其他参数以查看它们如何影响程序的空间和时间使用。通过调用count函数一万次来构建piTable

for i in xrange(1, 10000):
    piTable[i] = piTable[i-1] + \
        count(1000 * (i-1), 1000 * i)

到目前为止,所有计算都可以在编译时而不是运行时完成。当我在ideone.com 进行这些计算时,它们花费了大约 5 秒,但那段时间不算在内,因为当程序员第一次编写代码时,它可以永远完成一次。作为一般规则,您应该寻找机会将代码从运行时移动到编译时,以使您的程序运行得非常快。

剩下的就是写一个实际计算小于等于n的素数个数的函数:

def pi(n):
    if type(n) != int and type(n) != long:
        raise TypeError('must be integer')
    if n < 2: return 0
    if n == 2: return 1
    i = n // 1000
    return piTable[i] + count(1000 * i, n+1)

第一个if 语句进行类型检查。第二个if 语句返回对荒谬输入的正确响应。第三条if 语句专门处理2;我们的筛分使 1 成为素数,2 成为合数,两者都不正确,因此我们在此处进行修复。然后将 i 计算为小于请求 n 的 piTable 的最大索引,return 语句将 piTable 中的值与表值和请求值之间的素数计数相加; hi 限制是 n+1,否则在 n 是素数的情况下,它不会被计算在内。例如,说:

print pi(6543223)

将导致数字 447519 显示在终端上。

pi 函数非常快。在ideone.com,大约半秒内计算了对 pi(n) 的一千次随机调用,因此每个调用大约半毫秒;这包括生成素数和求和结果的时间,因此实际计算 pi 函数的时间甚至不到半毫秒。这对我们在建表方面的投资来说是一个相当不错的回报。

如果您对使用素数进行编程感兴趣,我在blog 上做了很多工作。请前来参观。

【讨论】:

  • 重新“制作一个包含 200 个变化点的表格”——即使误差很小,它也不是单调的,因此可以在 n 高达 10^7 的情况下更改超过 200 次。你数过变化的次数吗?
  • 不,我没有。感谢您的指正。无论如何,更改的数量很少,因此只需要一点空间。不过,我仍然认为筛分更可取。
  • 在 3 和 1000 之间,floor(Li(x)) - pi(x) 的值改变了 293 次,而 floor(Li(x) + .5) - pi(x) 改变了值260 次。这是行不通的。
  • 这远远超出我的想象。我同意这行不通。那么,让我们继续使用分段筛。
  • 我喜欢#1,它让人想起了解决不同问题的巨型逐步算法。像该算法一样,我希望您可以获得大约 sqrt(N) 的存储和工作。但是,有一个问题:您如何找到包含该值的两个表条目?我认为您可以使用素数定理进行猜测,然后进行更正(也许只是线性探测以找到括号条目)。
【解决方案3】:

我建议在这里使用启发式混合模型。存储每个第 n 个素数,然后通过素数测试进行线性搜索。为了加快速度,您可以使用快速而简单的素数测试(例如 a==2 的 Fermat 测试)并预先计算误报。根据输入的最大大小和存储限制进行一些微调应该很容易解决。

【讨论】:

  • 如果按照你说的实现的话,两个数之间会有大约 60 个素数。
  • 还有大约 800 个通过测试的小于 10^8 的非素数。计算 2^n (mod n) 是 O(log n);在实际水平上,这意味着要测试约 1000 个素数来获得计数。这种方法更适合一次性测试而不是批量测试,但仍然应该很快。
  • 这是个好主意。您可以通过使用 Miller-Rabin 的确定性变体使其变得更好并避免错误和 Carmichael 数;检查证人 2、7 和 61 适用于任何少于几十亿的输入。
  • 这是一个空间与时间的权衡;确定性 Fermat 测试非常快(并且易于实现)以至于难以被击败,并且存储 2064 个数字的成本非常小。 (另外,事实证明 2064 是小于 10^8 的正确数字;850 小于 10^7。如果我们测试 2 和 3,那么它会下降到 492)。
  • Miller-Rabin 一样快,之后您不需要“备份”表查找或试用除法。当您可以使用 Miller-Rabin 代替时,没有理由使用直接 Fermat 检验。
【解决方案4】:

您的建议是最好的。预先计算(或download)小于 10^7 的素数列表,然后对它们进行二进制搜索。

只有 664579 个小于 10^7 的素数,因此该列表将占用约 2.7 MB 的空间。解决每个实例的二进制搜索将非常快速 - 只需约 20 次操作。

【讨论】:

  • 但他说他的空间限制为 50k
  • OP writes "50k 我实际上是指源文件限制"
【解决方案5】:

这里有一些有效的代码。您应该使用适用于您的输入范围的确定性Miller-Rabin 测试替换基于试验划分的素数测试。在适当的小范围内筛选素数会比试除法更好,但这是朝着错误方向迈出的一步。

#include <stdio.h>
#include <bitset>
using namespace std;

short smallprimes[549]; // about 1100 bytes
char in[19531]; // almost 20k

// Replace me with Miller-Rabin using 2, 7, and 61.
int isprime(int j) {
 if (j<3) return j==2;
 for (int i = 0; i < 549; i++) {
  int p = smallprimes[i];
  if (p*p > j) break;
  if (!(j%p)) return 0;
 }
 return 1;
}

void init() {
 bitset<4000> siv;
 for (int i = 2; i < 64; i++) if (!siv[i])
  for (int j = i+i; j < 4000; j+=i) siv[j] = 1;
 int k = 0;
 for (int i = 3; i < 4000; i+=2) if (!siv[i]) {
  smallprimes[k++] = i;
 }

 for (int a0 = 0; a0 < 10000000; a0 += 512) {
  in[a0/512] = !a0;
  for (int j = a0+1; j < a0+512; j+=2)
   in[a0/512] += isprime(j);
 }
}

int whichprime(int k) {
 if (k==2) return 1;
 int a = k/512;
 int ans = 1 + !a;
 for (int i = 0; i < a; i++) ans += in[i];
 for (int i = a*512+1; i<k; i+=2) ans += isprime(i);
 return ans;
}

int main() {
 int k;
 init();
 while (1 == scanf("%i", &k)) printf("%i\n", whichprime(k));
}

【讨论】:

  • 我喜欢这个——定期存储一个小于 x 的素数的表格比我存储素数本身的方法更有效,因为查找是微不足道的。
【解决方案6】:

以下听起来像是您正在寻找的内容。 http://www.geekviewpoint.com/java/numbers/index_of_prime。在那里你会找到代码和单元测试。由于您的列表相对较小(即10^7),它应该处理它。

基本上你会找到2n 之间的所有素数,然后计算所有小于n 的素数以找到索引。此外,如果n 不是素数,则函数返回-1

【讨论】:

    【解决方案7】:

    我就这样做过一次。写了一段代码,即给定n,可以很快找到第n个素数,最多n=203542528,所以约2e8。或者,它可以倒退,对于任何数字 n,可以知道有多少质数小于 n。

    使用数据库。我将所有素数存储到某个点(我的上限的 sqrt)。在您的情况下,这意味着您将所有素数存储到 sqrt(1e7)。其中有 446 个,您可以以压缩形式存储该列表,因为到该点的最大差异仅为 34。超出该点,存储每个第 k 个素数(对于某个 k 值)。然后快速筛子就足够了在很短的时间间隔内生成所有素数。

    所以在 MATLAB 中,要找到第 1e7 个素数:

    nthprime(1e7)
    ans =
       179424673
    

    或者,它可以找到小于 1e7 的素数个数:

    nthprime(1e7,1)
    ans =
          664579
    

    关键是,这样的数据库很容易构建和搜索。如果你的数据库不能超过 50k,应该没有问题。

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 2015-01-31
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2018-11-26
      • 2012-10-13
      • 1970-01-01
      相关资源
      最近更新 更多