在Python中抽样n个大于p的随机质数的最快方法是什么？

Question

我想创建一个包含n个大于p ~ 2⁵⁰的质数数据集。我希望这些质数不是连续的，而是有一些间隔，使得第i个和第(i+1)个质数之间的差距不仅仅是几个比特。

我在循环中使用Sympy的randprime(low, hi)函数。

p = [start]
for i in range(n):
    curr = int(randprime(start, 2 * start + 1))
    p.append(curr)
    start = curr

对于 n=10,000，这个过程会显著变慢。有没有更好（更快）的方法来实现我想要的质数抽样？

- scribe

---

并不是你的问题，但你能分享一下想要让素数“有些间隔”的原因吗？显然，这个问题可以通过反复调用“randprime(a,b)”来解决。 - Captain Trojan

就实现而言，我没有看到可以改进的地方（有点惊讶randprime需要转换为整数）。如果这对你来说太慢了，你可能需要完全改变方法。话虽如此，你所尝试做的事情本质上很慢，因为这不是一项容易的任务。 - Ma0

@CaptainTrojan 我正在进行一些研究，尝试学习算法。如果两个不同的质数只在一两个位上有所不同，它们实际上并没有教会算法什么。此外，附近的质数大约具有相同数量的0位，因此在分解时预测这些数字很容易。为了真正测试收敛解决方案，我需要那些在位模式上看起来真正不同的质数。 - scribe

基本上只有两种方法可以生成质数：记忆化（基本上是埃拉托斯特尼筛法）或不记忆化（即独立测试每个数字）。第一种方法需要内存，第二种方法需要时间。顺便说一句，Sympy的素性测试非常快，所以......基本上就是Roman Pavelka所说的。 - Amadan

我建议使用厄拉多塞筛法构建质数，然后从中抽样。 - oskros

所有的大质数都可以表示为6n±1的形式，这可能是加速搜索的一种方法。 - rossum

Answer 1

预先计算（或只需下载）一次质数列表，然后在运行时从列表中取样。

您可以使用此算法为小于2 ^ 50（~10 ^ 15）的上限生成列表，该上限比30倍小：https://primes.utm.edu/nthprime/algorithm.php

我不知道如何在合理的硬件设置下进一步进行。

- Roman Pavelka

---

没有各种大小的质数列表。2 ** 50大约比一万亿大1000倍，我怀疑是否存在如此大的质数列表。问题在于如何计算这样的列表。 - President James K. Polk

@PresidentJamesK.Polk，你说得很好，我已经编辑了我的回答。谢谢你指出这一点。 - Roman Pavelka

说句实话，前一万亿个质数可以在这里找到：http://compoasso.free.fr/primelistweb/page/prime/liste_online_en.php - David Kaftan

Answer 2

当您使用最后一个值（curr）来确定下一个范围时，您正在指数级地增加范围。平均而言，随机质数应该落在范围的中间，并且范围从X到2X。这将使每次迭代将范围向前移动大约1.5X的因素。通过10,000次迭代，您的范围将增加高达1.5 ^ 10000（2 ^ 5850）的因素，这很快将使即使sympy也很难产生质数。

如果您的目标仅是在第i个和（i + 1）个质数之间具有足够数量的不同位，那么您可以保持相同数量级，并在先前的质数上过滤出最少数量的不同位（而不是增加随机范围的数量级）。

例如：

def oneBits(N): return N%2 + oneBits(N//2) if N else 0

minDiff  = 25 # minimum number of differing bits from ith to (i+1)th
minValue = 2**50
maxValue = minValue*2-1 
primes   = [randprime(minValue,maxValue)]
count    = 10000
for _ in range(count-1):
    while True:
       p = randprime(minValue,maxValue)
       if oneBits(primes[-1]^p)>=minDiff: break
    primes.append(p)

请注意，我没有安装sympy，因此我使用纯随机数进行了不同的测试，以获取下一个质数：

内存高效的质数生成器：

def genPrimes(toN):
    skips   = dict()
    maxSkip = int(toN**0.5)
    if toN>=2: yield 2
    for p in range(3,toN+1,2):
        if p not in skips:
            yield p
            if p <= maxSkip: skips[p*p] = 2*p
        else:
            stride = skips.pop(p)
            multiple = p + stride
            while multiple in skips: multiple += stride
            skips[multiple] = stride

获取从N开始的下一个质数的函数：

primes = list(genPrimes(2**26)) # for 2^50 max base prime is √(2^51) ~ 2^26
def nextPrime(N):
    sieve     = [1]*N.bit_length()*20
    maxPrime  = int((N+len(sieve))**0.5)
    for p in primes:
        if p>maxPrime: break
        offset = (p-N%p)%p
        if offset>len(sieve): continue
        sieve[offset::p] = [0]*len(range(offset,len(sieve),p))
    for p,isPrime in enumerate(sieve,N):
        if isPrime: return p
    return nextPrime(N+len(sieve))

最小间隔50位质数（生成器）：

import random
def randomPrimes(count,bits=50):
    minVal    = 2**bits
    maxVal    = minVal*2-1
    minDiff   = bits//2
    prevPrime = 0
    
    def oneBits(N): return N%2 + oneBits(N//2) if N else 0
    for _ in range(count):
        while True:
            n = random.randint(minVal,maxVal)            
            if prevPrime and oneBits(prevPrime^n)<minDiff: continue
            n = nextPrime(n)
            if not prevPrime: break
            if oneBits(prevPrime^n)>=minDiff: break
        yield n

输出：

for p in randomPrimes(10000):
    print(p,f"{p:b}")
            
2071968049418461 111010111000111000110100111100100100101000011011101
1399795190350597 100111110010001101100110111000101000100001100000101
1530818178259709 101011100000100010101100001101110110001101011111101
1140103670657957 100000011001110101100010010010010111111001110100101
1911908333932859 110110010101101111011011001000101100101000100111011
1236033889571977 100011001000010101010010000111010110001010010001001
1752989992684999 110001110100101010111001001110011110000110111000111
1849449158362859 110100100100001000001110000000111010100011011101011
1349704431776567 100110010111000110010001101001101010011001100110111
1142235147712271 100000011101101101101011000001110101011011100001111
...

这大约需要0.5秒钟来完成每一次迭代（无论迭代次数多少）。我相信sympy应该比我自己编写的函数要快得多。

- Alain T.