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.

About these ads

Pages: 1 2

6 Responses to “Frobenius Primality Test”

  1. Paul said

    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

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

  3. Paul said

    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

    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

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

Follow

Get every new post delivered to your Inbox.

Join 616 other followers

%d bloggers like this: