Frobenius Primality Test

November 13, 2012

Our code is a straight forward translation of the algorithm:

(define (frobenius-pseudoprime? n)
  (let loop ((a (randint 1 n)) (b (randint 1 n)))
    (let ((d (- (* a a) (* 4 b))))
      (if (or (= a b) (square? d) (not (= (gcd (* 2 a b d) n) 1)))
          (loop (randint 1 n) (randint 1 n))
          (let* ((m (/ (- n (jacobi d n)) 2))
                 (w (modulo (- (* a a (inverse b n)) 2) n)))
            (let loop ((u 2) (v w) (ms (digits m 2)))
              (if (pair? ms)
                  (if (zero? (car ms))
                      (loop (modulo (- (* u u) 2) n)
                            (modulo (- (* u v) w) n) (cdr ms))
                      (loop (modulo (- (* u v) w) n)
                            (modulo (- (* v v) 2) n) (cdr ms)))
                  (if (positive? (modulo (- (* u w) (* 2 v)) n)) #f
                    (let ((x (expm b (/ (- n 1) 2) n)))
                      (zero? (modulo (- (* u x) 2) n)))))))))))

To turn this into a proper primality checker, we combine trial division by the primes less than a hundred, which weeds out many trivial composites, with Miller-style strong pseudoprime tests to bases 2 and 3, and then finally a Lucas/Frobenius pseudoprime test.

(define (prime? n)
  (let loop ((ps (list 2 3 5 7 11 13 17 19 23 29 31 37
              41 43 47 53 59 61 67 71 73 79 83 89 97)))
    (if (pair? n)
        (if (zero? (modulo n (car ps)))
            (= n (car ps))
            (loop (cdr ps)))
        (and (strong-pseudoprime? n 2)
             (strong-pseudoprime? n 3)
             (frobenius-pseudoprime? n)))))

Here are two examples:

> (prime? 3825123056546413051)
#f
> (prime? (- (expt 2 89) 1))
#t

You can run the program at http://programmingpraxis.codepad.org/mbEWx4S1. The version of the program shown there provides a host of supporting functions and gives the details of the calculation of the Lucas chain, which may be useful for debugging.

Pages: 1 2

Posted by programmingpraxis
Filed in Exercises
6 Comments »

6 Responses to “Frobenius Primality Test”

  1. Paul said
    November 13, 2012 at 3:34 PM

    A Python version. I had some problems with the formatting of the text in the exercise. Especially symbols as != seem to be missing. The debugging information given was most helpful. A full version is here.

    def frobenius(n):
    """ returns True if n is probably prime
    """
    if n % 2 == 0 or n < 2:
        return False
    #step 1
    while 1:
        a = randint(1, n-1)
        b = randint(1, n-1)
        d = a ** 2 - 4 * b
        if  not (is_square(d) or gcd(2 * a  * b * d, n) != 1):
            break
    #step 2
    m = (n - jacobi(d, n)) // 2
    v = w = (a ** 2  * inverse(b, n) - 2) % n
    u = 2
    # step 3
    for bit in bits_of(m):
        if bit:
            u = (u * v - w) % n 
            v = (v * v - 2) % n
        else:
            v = (u * v - w) % n
            u = (u * u - 2) % n
    # step 4
    if (w * u - 2 * v) % n != 0:
        return False
    return (pow(b, (n - 1) // 2, n) * u) % n ==  2
  2. programmingpraxis said
    November 13, 2012 at 3:40 PM

    Fixed the missing not-equals sign in Step 4. Thanks.

  3. Paul said
    November 13, 2012 at 3:50 PM

    There is also something missing at the end of Step 1 line 1.
    In step 2 v is set to 2. Should that not be v = w?

  4. programmingpraxis said
    November 13, 2012 at 4:22 PM

    Fixed both of those. Thanks. I didn’t do that very well, did I?

    I’m glad the debugging information was useful for you.

  5. Paul said
    November 13, 2012 at 4:38 PM

    Somtimes mistakes are made. Thanks a lot for the exercises. As you can see from my posts I am enjoying them a lot.

Leave a comment