Generating Large Random Primes
June 9, 2017
We begin by computing the sieving primes less than 65,000 (well, 65,536, but who’s counting?):
(define sievers (primes (expt 2 16)))
We randomly choose a range of size 100000, storing its lower bound in base. The two nested do loops perform the sieving; the inner do loop starts at the offset within the sieve of the smallest multiple of the current sieving prime, calculated by that strange-looking modulo expression, and increments by p. Then the inner-loop on i checks the numbers that survived the sieve for primality; if none are found, outer-loop chooses another base and restarts the calculation:
(define (rand-prime lo hi)
(define (rand-odd lo hi)
(let ((n (randint lo hi)))
(if (odd? n) n (+ n 1))))
(let outer-loop ((base (rand-odd lo (- hi 100000))))
(let ((sieve (make-vector 50000 #t)))
(do ((ps (cdr sievers) (cdr ps))) ((null? ps))
(let ((p (car ps)))
(do ((i (modulo (* -1/2 (+ base p)) p) (+ i p)))
((<= 50000 i)))))
(let inner-loop ((i 0))
(if (<= 50000 i) (outer-loop (rand-odd lo (- hi 100000)))
(if (and (vector-ref sieve i) (prime? (+ base i i))) (+ base i i)
(inner-loop (+ i 1))))))))
Here’s an example:
> (rand-prime #e1e15 #e1e16) 2382987027564103
We reused the random number generator, Sieve of Eratosthenes, and Miller-Rabin primality checker from the previous exercise. You can run the program at http://ideone.com/jc57fP. This algorithm will always be faster than the algorithm of the previous exercise, but there is no guarantee that the output is prime. Of course, it’s good enough for some purposes.
In Python.
SIEVE_SIZE = 50000 prs = list(primes(65000))[1:] def large_prime(lo, hi): while True: base = random.randrange(lo, hi - 2*SIEVE_SIZE) | 1 sieve = [1]*SIEVE_SIZE for p in prs: offset = (-(base + p) // 2) % p sieve[offset::p] = [0] * len(range(offset, SIEVE_SIZE, p)) for i, f in enumerate(sieve): if f: cand = base + 2 * i if is_prime(cand): # miller-rabin return cand