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

Pages: 1 2

### One Response to “Generating Large Random Primes”

1. Paul said

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 = *SIEVE_SIZE
for p in prs:
offset = (-(base + p) // 2) % p
sieve[offset::p] =  * 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
```