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

Advertisements

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.

6. […] Pages: 1 2 […]