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.

Advertisements

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 = [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
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: