在Python中采样大于p的n个随机素数的最快方法？答案

【问题标题】：Fastest way to sample n random primes greater than p in Python?在Python中采样大于p的n个随机素数的最快方法？
【发布时间】：2021-06-23 08:46:23
【问题描述】：

我想创建一个包含大于 p ~ 2⁵⁰ 的 n 个素数的日期集。我希望这些素数不是连续的，但它们之间有一些空间，所以 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，这会变得非常缓慢。有没有更好（更快）的方法来完成我想要的主要采样？

【问题讨论】：

不是你的问题，但是你能不能通过将想要你的素数的原因分享给have some space in between来满足我的好奇心，这显然是通过反复调用randprime(a, b)来解决的？
在实现方面，我没有看到任何可以改进的地方（对于必须将 randprime 强制转换为 int 有点惊讶）。如果这对您来说太慢，您可能需要一起更改方法。话虽如此，您尝试做的事情本身就很慢，因为这不是一件容易的事。
@CaptainTrojan 我正在做一些研究，我正在尝试学习算法，如果两个不同的素数仅在一位或两位上有所不同，那么它们并没有真正教给算法任何东西。此外，附近的素数具有大约相同数量的 0 位，因此在因式分解时预测这些素数很容易。为了真正测试一个收敛的解决方案，我想要在位模式中看起来不同的素数。
基本上只有两种生成素数的方法：memoise（基本上是埃拉托色尼筛法）或 not memoise（即单独测试每个数字）。第一个需要记忆，第二个需要时间。 FWIW，Sympy 的素性测试非常快，所以...基本上，正如 Roman Pavelka 所说。
我建议使用 Eratosthenes 筛构建素数，然后从中取样

标签： python primes sampling

【解决方案1】：

预计算（或只是download）一次素数列表，然后在运行时对列表进行采样。

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

我不知道如何通过合理的硬件设置走得更远。

【讨论】：

没有所有大小的素数列表。 2**50 大约是一万亿的 1000 倍，我怀疑那里有这种大小的素数的列表。问题是如何计算这样一个列表。
@PresidentJamesK.Polk，这是一个很好的观点，我已经编辑了我的答案。感谢您指出这一点。
不管怎样，前万亿素数都可以在这里找到：compoasso.free.fr/primelistweb/page/prime/liste_online_en.php

【解决方案2】：

当您使用最后一个值 (curr) 来确定下一个范围时，您正在以指数方式增加范围。平均而言，随机素数应该落在范围的中间，范围从 X 到 2X。这将在每次迭代时将范围向前移动大约 1.5 倍。通过 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 应该比我的自制函数快得多。

【讨论】：