Fastest way to sample n random primes greater than p in Python?

user_12035739 0
Viewed 156

I want to create a dateset of n primes greater than p ~ 2⁵⁰. I want these primes to not be consecutive but have some space in between so the difference between iᵗʰ and (i+1)ᵗʰ prime is not just a few bits.

I am using Sympy's randprime(low, hi) in a loop,

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

This gets significantly slow for n=10,000. Is there a better (faster) way of accomplishing the prime sampling that I want?

sampling primes python
2 Answers

Precompute (or just download) a list of primes once and then sample the list during the runtime.

You can generate the list for 30x smaller upper bound than 2^50 (~10^15) with this algorithm: https://primes.utm.edu/nthprime/algorithm.php

I don't know how to get further with reasonable hardware setup.

user_121185460

When you use the last value (curr) to determine the next range, you are increasing the range exponentially. On average the random prime should fall about midway of the range and the range goes from X to 2X. This will move the range forward by a factor of roughly 1.5X at each iteration. With 10,000 iterations your range will increase by a factor of up to 1.5^10000 (2^5850) which will soon make it very hard for even sympy to produce primes.

If your objective is merely to have a sufficient number of differing bits between the ith and (i+1)th prime, you could stay in the same order of magnitude and filter on a minimal number of distinct bits with the previous prime (instead of increasing the magnitude of the random range).

for example:

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)

Note that I don't have sympy so I tested this a little differently using purely random numbers for which I get the next prime:

Memory efficient primes generator:

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

Function to get next prime from 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))

Minimally spaced 50-bit primes (generator):

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

output:

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
...

This steadily takes roughly 0.5 second per iteration (no matter how many iterations are done). I believe sympy should be much faster than my home made fucntions.

user_52375600
Related