Generating Random Large Primes

March 31, 2017

Here is our version of Pocklington’s test, which takes p, k and a and returns N if N is prime or #f if N is composite:

(define (pocklington p k a)
  (and (odd? p) (prime? p)
       (< 1 k (* 2 (+ p 1)))
       (positive? (modulo k p))
       (let ((n (+ (* 2 k p) 1)))
         (and (< 1 a n)
              (= (expm a (/ (- n 1) 2) n) (- n 1))
              (= (gcd (+ (expm a k n) 1) n) 1)
              n)))) ; return n if found, else #f

> (pocklington 2333 2001 2)
9336667
> (pocklington 2333 2017 2)
#f

Next we have function ratchet that repeatedly tests combinations of k and a until it finds a prime. Instead of using the pocklington function given above, we use a smaller internal function, because the full function does too much work:

(define (ratchet p lo hi verbose?)
  (define (pock p k a n)
    (and (= (expm a (/ (- n 1) 2) n) (- n 1))
         (= (gcd (+ (expm a k n) 1) n) 1)))
  (let loop ((k (randint lo hi)))
    (let ((n (+ (* 2 k p) 1)))
      (cond ((pock p k 2 n) (when verbose? (display (list p k 2)) (newline)) n)
            ((pock p k 3 n) (when verbose? (display (list p k 3)) (newline)) n)
            (else (loop (randint lo hi)))))))

Now all that’s left is to work the ratchet. Note the computation of the range of k, which may not be 1 ≤ k < 2 × p on the last ratcheting:

(define (rand-prime lo hi . args)
  (let ((verbose? (if (pair? args) (car args) #f)))
    (if (< hi 10000) (fortune (drop-while (lambda (x) (< x lo)) (primes hi)))
      (let loop ((p (fortune (cdr (primes 10000)))))
        (if (<= lo p) p
          (let ((k-lo (ceiling (/ (- lo 1) 2 p)))
                (k-hi (floor (/ (- hi 1) 2 p))))
            (loop (ratchet p (if (< (* 2 p) k-lo) 1 k-lo)
                           (min (* 2 p) k-hi) verbose?))))))))

We are careful to handle small ranges properly, and we choose the initial prime randomly from the primes less than ten thousand. Here’s an example computation of a 50-digit prime:

> (rand-prime #e1e49 #e1e50 #t)
(607 1022 2)
(1240709 147645 2)
(366368960611 227126727978 2)
(166424366512554387349117 236422002865808466771968 3)
(78692764113142983989384966585012763035290304513 618 2)
97264256443844728210879818699075775111618816378069

You can run the program at http://ideone.com/eVR8RA, where you will also see contributions from the Standard Prelude, a Sieve of Eratosthenes, a Miller-Rabin prime number checker, and the 147 random number generator from the previous exercise.

The function given above can fail if the range for k on the last ratchet is small and yields no primes. At an interactive prompt, simply interrupt the computation and try again. For a non-interactive program, it would be necessary to timeout and restart a runaway computation, but we won’t do that because there is no standard timer in Scheme.

By the way, Eric Weisstein’s entry on Pocklington’s Criterion at Wolfram’s MathWorld is wrong; it omits some of the conditions on a. I sent a message to notify MathWorld of the error.

Posted by programmingpraxis

Filed in Exercises

18 Comments »

18 Responses to “Generating Random Large Primes”

Milbrae said

April 2, 2017 at 5:53 PM

Simplistic Python v3.x

from __future__ import division
from random import randrange

def gcd(m, n):
    while n:
        r = m % n
        m = n
        n = r
    return m

def simplePrimeTrial(n):
    if n < 5:
        return (n == 2) or (n == 3)

    if (n % 2 == 0) or (n % 3 == 0) or (n % 7 == 0):
        return False

    d = 7
    while d * d <= n:
        if n % d == 0:
            return False
        d += 2

    return True

#
if __name__ == "__main__":

    prime = False

    while not prime:
        # choose odd prime p
        p = randrange(10**6, 10**9) | 1
        while not simplePrimeTrial(p):
            p += 2

        # choose 1 < k < 2(p + 1)
        k = 0
        while (k % p == 0):
            k = randrange(2, 2 * (p + 1))

        # get N
        N = 2 * k * p + 1

        # conditions
        for a in [2, 3]:
            if pow(a, (N - 1) // 2, N) == (N - 1):
                if gcd(a**k, N) == 1:
                    prime = True
                    break

        if prime:
            break

    print (p, k, N, prime)

Milbrae said
April 2, 2017 at 5:54 PM
Forgot to mention a possible result:
```
805103465 393971070 634374947133515101 True
```

Milbrae said

April 2, 2017 at 6:19 PM

The trial division routine is completely out of order. Sorry! Here’s the corrected procedure:

def simplePrimeTrial(n):
    if n < 6:
        return (n == 2) or (n == 3) or (n == 5)

    if (n % 2 == 0) or (n % 3 == 0) or (n % 5 == 0):
        return False

    d = 7
    while d * d <= n:
        if n % d == 0:
            return False
        d += 2

    return True

programmingpraxis said

April 3, 2017 at 12:01 AM

Here is my favorite simple primality checker:

def isPrime(n):
    wheel = [1,2,2,4,2,4,2,4,6,2,6]
    w, f, fs = 0, 2, []
    while f*f <= n:
        if n % f == 0:
            return False
        w, f = w + 1, f + wheel[w]
        if w == 11: w = 3
    return True

It uses a prime wheel, so it should be about twice as fast as yours.

Milbrae said
April 3, 2017 at 8:53 AM
@programmingpraxis: Awesome! Thanks you!

Also, I’vve found another bug in my source. The conditions check should read like this:
```
# conditions
        for a in [2, 3]:
            if pow(a, (N - 1) // 2, N) == (N - 1):
                if gcd(a**k + 1, N) == 1:
                    prime = True
                    break
```
The GCD bit still bothers me. If k and N are large and of different size, this will take forever and possibly raise a MemoryError.
programmingpraxis said
April 3, 2017 at 12:51 PM
Your gcd is fine, and will execute quickly; Knuth guarantees it. Most Python programmers would write it like this:
```
def gcd(m, n):
    while n <> 0:
        m, n = n, m % n
    return m
```

programmingpraxis said

April 3, 2017 at 12:58 PM

@Millbrae: By the way, I blindly copied the prime checker from my wheel factorization program, and it does too much work. Here’s a corrected version:

def isPrime(n):
    wheel = [1,2,2,4,2,4,2,4,6,2,6]
    w, f = 0, 2
    while f * f <= n:
        if n % f == 0:
            return False
        w, f = w+1, f + wheel[w]
        if w == 11: w = 3
    return True

And here is the corresponding function that factors integers, which is what I blindly copied from:

def factors(n):
    wheel = [1,2,2,4,2,4,2,4,6,2,6]
    w, f, fs = 0, 2, []
    while f * f <= n:
        while n % f == 0:
            fs.append(f)
            n /= f
        w, f = w + 1, f + wheel[w]
        if w == 11: w = 3
    fs.append(n)
    return fs

All these functions have appeared on Programming Praxis previously.

Milbrae said
April 3, 2017 at 1:13 PM
Thanks again, PP!

Yet there has to be some faster way to get GCD(a**k + 1, N) checked – if k and N are large. Python and PyPy take forever to finish for number 10**12 < n < 10**18.
[Quote from the description]Henry Pocklington’s Criterion, which dates to 1914, gives us a way to find such primes quickly.[/Quote]
… quickly … Quickly … QUICKLY
programmingpraxis said
April 3, 2017 at 1:18 PM
@Millbrae: If that gcd isn’t fast enough, you might want to look at this.
Paul said
April 3, 2017 at 2:07 PM
@Millbrae: If you do gcd(pow(a, k+1, N), N) there are no memory problems and it runs a lot faster.

Milbrae said

April 3, 2017 at 3:33 PM

Thanks a lot, Paul! I’ve been thinking abuot something like this:

def euclid_gcd_a_power_k_plus_1_n(a, k, n):
    op = (pow(a, k, n) + 1) % n
    while op:
        r = n % op
        n = op
        op = r
    return n

Haven’t tested it though…

Milbrae said

April 3, 2017 at 4:25 PM

Latest code in Python:

from __future__ import division
from random import randrange
from math import sqrt
from time import clock

def gcd(a, b):
    while b:
        r = a % b
        a = b
        b = r
    return a

def binary_gcd(a, b):
    shift = 0

    while a and b:
        if (a & 1):
            if (b & 1):
                a, b = (a - b, b) if a > b else (a, b - a)
            else:
                b >>= 1
        elif (b & 1):
            a >>= 1
        else:
            shift += 1
            a >>= 1
            b >>= 1

    return (a or b) << shift

def euclid_gcd_a_power_k_plus_1_n(a, k, n):
    op = (pow(a, k, n) + 1) % n
    while op:
        r = n % op
        n = op
        op = r
    return n

def isPrime1(n):
    wheel = [1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6]
    w, f = 0, 2
    while f * f <= n:
        if n % f == 0:
            return False
        w, f = w + 1, f + wheel[w]
        if w == 11: w = 3
    return True

def isPrime2(n):
    if n < 5:
        return (n == 2) or (n == 3)

    if (not (n & 1)) or (n % 3 == 0):
        return False

    limit = int(sqrt(n))
    d = 5
    i = 2
    while d <= limit:
        if n % d == 0: return False
        d += i
        i = 6 - i

    return True

#
if __name__ == "__main__":

    prime = False

    t = clock()

    while not prime:
        # choose odd prime p
        p = randrange(10**15, 10**18) | 1
        while not isPrime1(p):
            p += 2

        # choose 1 < k < 2(p + 1)
        k = 0
        while (k % p == 0):
            k = randrange(2, 2 * (p + 1))

        # get N
        N = 2 * k * p + 1

        # intermediate display
        print ("p = %d / k = %d" % (p, k))

        # conditions
        for a in [2, 3]:
            if pow(a, (N - 1) // 2, N) == (N - 1):
                #if euclid_gcd_a_power_k_plus_1_n(a, k, N) == 1:
                if gcd(pow(a, k + 1, N), N) == 1:
                    prime = True
                    break

        if prime:
            break

    print ()
    print ("--> p = %d, k = %d, N = %d, prime = %s" % (p, k, N, prime))
    print ("%f seconds" % (clock() - t))

Sample output:

[...]
--> p = 119053509950371751, k = 166652094649708243, N = 39681033617258670310871283078086987, prime = True
362.305565 seconds

I’ve tested N at https://www.alpertron.com.ar/ECM.HTM. Everything’s turned out fine.

Milbrae said

April 3, 2017 at 4:48 PM

Additionally, some numbers. If you manually set p = 4089548364052918515677660474680192167456016670586738083090737083455827, the above code will output…

--> p = 4089548364052918515677660474680192167456016670586738083090737083455827, k = 3954952269115446666063958491090918288547649049955194620953816304862234, N = 32347937164136905687977400799521711490238478133100320456741932622470671464996525396091254589705022471721754863572714447927950109639719075037, prime = True

N is then a 140-digit prime ;)

programmingpraxis said

April 3, 2017 at 4:59 PM

@Milbrae: It should not take six minutes to generate a random 35-digit prime. You seem to misunderstand the algorithm. Start with a small prime, and use it to generate a larger prime. Then recursively use the larger prime that you just generated to generate an even larger prime. And so on. Continue the “ratcheting” process until the end result is the desired size. My Scheme implementation generates a random 35-digit prime in 15ms, which is four orders of magnitude faster, and generates a 140-digit prime in a little bit less than a second:

> (time (rand-prime #e1e34 #e1e35 #t))
(607 1022 2)
(1240709 147645 2)
(366368960611 227126727978 2)
(166424366512554387349117 133785462153 2)
(time (rand-prime 10000000000000000000000000000000000 ...))
    no collections
    15 ms elapsed cpu time
    21 ms elapsed real time
    2632920 bytes allocated
44530321574804691176098788802937803

> (time (rand-prime #e1e139 #e1e140 #t))
(4957 2898 2)
(28730773 19790703 3)
(1137204390806839 1366720362160076 3)
(3108480793707103200437027119529 768273308507447641982667587584 3)
(4776325647626426015624444714829864966096790270290693688655873 16807863103533598
20583032452649851040365352536001193141338112 3)
(1605596552464028483036164591830457063463238223076044024630019667995991481919686
4066663901322728680984541663162851815063553 2679361547405450752 3)
(time (rand-prime 10000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000 .
..))
    20 collections
    968 ms elapsed cpu time, including 0 ms collecting
    981 ms elapsed real time, including 1 ms collecting
    83927816 bytes allocated, including 84370688 bytes reclaimed
86039473266377526955299780382198183524627420317908841567482428065070282459178371
593009228697984643733130344190244139690296472557146983283713

Milbrae said

April 3, 2017 at 6:03 PM

Well, then maybe this will suffice…

from __future__ import division
from random import randrange
from time import clock

def gcd(a, b):
    while b:
        r = a % b
        a = b
        b = r
    return a

def isPrime(n):
    wheel = [1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6]
    w, f = 0, 2
    while f * f <= n:
        if n % f == 0:
            return False
        w, f = w + 1, f + wheel[w]
        if w == 11: w = 3
    return True

def generate_large_prime(seed = 127, length = 30):
    p = seed | 1
    while not isPrime(p):
        p += 2

    cnt = 1
    while len(str(p)) < length:
        k = 0
        while k % p == 0:
            k = randrange(2, 2 * (p + 1))

        N = 2 * k * p + 1

        for a in [2, 3]:
            if pow(a, (N - 1) // 2, N) == (N - 1):
                if gcd(pow(a, k + 1, N), N) == 1:
                    print ("%3d. %d" % (cnt, N))
                    cnt += 1
                    p = N
                    break

    return None

#
if __name__ == "__main__":

    t = clock()
    generate_large_prime(1000, 200)
    print ("\n%f seconds" % (clock() - t))

Sample output:

  1. 2149171
  2. 9834713954551
  3. 104610237726493021241396621
  4. 34200932046364083964639294647733262301861498990032657
  5. 396323208934679755427233521579881191053730646550572161988585492317517817803
0298810631751649309406392290359
  6. 28909194903683721538230094312229715957688421171095500854838580090229909584920332496482513814271432996604585263726247834081029062687390514758431875294893376647359646225385355346599037139194196446020843723285340689

0.393686 seconds

programmingpraxis said
April 3, 2017 at 6:13 PM
@Milbrae: Much faster. But the result is 212 digits, not the requested 200. Look at my computation of k-lo and k-hi.

Milbrae said

April 3, 2017 at 6:19 PM

Or slightly changed, showing only the generated prime (and its length)

[...]
def generate_large_prime(seed = 127, length = 30):
    p = seed | 1
    while not isPrime(p):
        p += 2

    while len(str(p)) < length:
        k = 0
        while k % p == 0:
            k = randrange(2, 2 * (p + 1))

        N = 2 * k * p + 1

        for a in [2, 3]:
            if pow(a, (N - 1) // 2, N) == (N - 1):
                if gcd(pow(a, k + 1, N), N) == 1:
                    p = N
                    break

    return N

#
if __name__ == "__main__":

    t = clock()
    lp = generate_large_prime(100, 150)
    print ("%d (%d digits)" % (lp, len(str(lp))))
    print ("\n%f seconds" % (clock() - t))

Sample output:

49147198556389931870735198516291966016405112299420103881855365916653844771408104584372776813914932321850903847461556401822236286680702008668299053716265038489973603779363979922476355643397827803871258431258822366365509 (218 digits)

0.206664 seconds

Generating Large Random Primes | Programming Praxis said
June 9, 2017 at 2:01 PM
[…] studied Pocklington’s Criterion, which lets us quickly find large random primes, in a previous exercise. That algorithm generates a certified prime — a number that is proven to be prime — rather than […]

Programming Praxis