How to implement Frobenius pseudoprime algorithm?

Question

Someone told me that the Frobenius pseudoprime algorithm take three times longer to run than the Miller–Rabin primality test but has seven times the resolution. So then if one where to run the former ten times and the later thirty times, both would take the same time to run, but the former would provide about 233% more analyse power. In trying to find out how to perform the test, the following paper was discovered with the algorithm at the end:

A Simple Derivation for the Frobenius Pseudoprime Test

There is an attempt at implementing the algorithm below, but the program never prints out a number. Could someone who is more familiar with the math notation or algorithm verify what is going on please?

Edit 1: The code below has corrections added, but the implementation for compute_wm_wm1 is missing. Could someone explain the recursive definition from an algorithmic standpoint? It is not "clicking" for me.

Edit 2: The erroneous code has been removed, and an implementation of the compute_wm_wm1 function has been added below. It appears to work but may require further optimization to be practical.

from random import SystemRandom
from fractions import gcd
random = SystemRandom().randrange

def find_prime_number(bits, test):
    number = random((1 << bits - 1) + 1, 1 << bits, 2)
    while True:
        for _ in range(test):
            if not frobenius_pseudoprime(number):
                break
        else:
            return number
        number += 2

def frobenius_pseudoprime(integer):
    assert integer & 1 and integer >= 3
    a, b, d = choose_ab(integer)
    w1 = (a ** 2 * extended_gcd(b, integer)[0] - 2) % integer
    m = (integer - jacobi_symbol(d, integer)) >> 1
    wm, wm1 = compute_wm_wm1(w1, m, integer)
    if w1 * wm != 2 * wm1 % integer:
        return False
    b = pow(b, (integer - 1) >> 1, integer)
    return b * wm % integer == 2

def choose_ab(integer):
    a, b = random(1, integer), random(1, integer)
    d = a ** 2 - 4 * b
    while is_square(d) or gcd(2 * d * a * b, integer) != 1:
        a, b = random(1, integer), random(1, integer)
        d = a ** 2 - 4 * b
    return a, b, d

def is_square(integer):
    if integer < 0:
        return False
    if integer < 2:
        return True
    x = integer >> 1
    seen = set([x])
    while x * x != integer:
        x = (x + integer // x) >> 1
        if x in seen:
            return False
        seen.add(x)
    return True

def extended_gcd(n, d):
    x1, x2, y1, y2 = 0, 1, 1, 0
    while d:
        n, (q, d) = d, divmod(n, d)
        x1, x2, y1, y2 = x2 - q * x1, x1, y2 - q * y1, y1
    return x2, y2

def jacobi_symbol(n, d):
    j = 1
    while n:
        while not n & 1:
            n >>= 1
            if d & 7 in {3, 5}:
                j = -j
        n, d = d, n
        if n & 3 == 3 == d & 3:
            j = -j
        n %= d
    return j if d == 1 else 0

def compute_wm_wm1(w1, m, n):
    a, b = 2, w1
    for shift in range(m.bit_length() - 1, -1, -1):
        if m >> shift & 1:
            a, b = (a * b - w1) % n, (b * b - 2) % n
        else:
            a, b = (a * a - 2) % n, (a * b - w1) % n
    return a, b

print('Probably prime:\n', find_prime_number(300, 10))

Answer 1

You seem to have misunderstood the algorithm completely due to not being familiar with the notation.

def frobenius_pseudoprime(integer):
    assert integer & 1 and integer >= 3
    a, b, d = choose_ab(integer)
    w1 = (a ** 2 // b - 2) % integer

That comes from the line

W₀ ≡ 2 (mod n) and W₁ ≡ a²b⁻¹ − 2 (mod n)

But the b^-1 doesn't mean 1/b here, but the modular inverse of b modulo n, i.e. an integer c with b·c ≡ 1 (mod n). You can most easily find such a c by continued fraction expansion of b/n or, equivalently, but with slightly more computation, by the extended Euclidean algorithm. Since you're probably not familiar with continued fractions, I recommend the latter.

    m = (integer - d // integer) // 2

comes from

n − (∆/n) = 2m

and misunderstands the Jacobi symbol as a fraction/division (admittedly, I have displayed it here even more like a fraction, but since the site doesn't support LaTeX rendering, we'll have to make do). The Jacobi symbol is a generalisation of the Legendre symbol - denoted identically - which indicates whether a number is a quadratic residue modulo an odd prime (if n is a quadratic residue modulo p, i.e. there is a k with k^2 ≡ n (mod p) and n is not a multiple of p, then (n/p) = 1, if n is a multiple of p, then (n/p) = 0, otherwise (n/p) = -1). The Jacobi symbol lifts the restriction that the 'denominator' be an odd prime and allows arbitrary odd numbers as 'denominators'. Its value is the product of the Legendre symbols with the same 'numerator' for all primes dividing n (according to multiplicity). More on that, and how to compute Jacobi symbols efficiently in the linked article. The line should correctly read

m = (integer - jacobi_symbol(d,integer)) // 2

The following lines I completely fail to understand, logically, here should follow the calculation of W_m and W_m+1 using the recursion

W_2j ≡ W_j² − 2 (mod n)

W_2j+1 ≡ W_jW_j+1 − W₁ (mod n)

An efficient method of using that recursion to compute the required values is given around formula (11) of the PDF.

    w_m0 = w1 * 2 // m % integer
    w_m1 = w1 * 2 // (m + 1) % integer
    w_m2 = (w_m0 * w_m1 - w1) % integer

The remainder of the function is almost correct, except of course that it now gets the wrong data due to earlier misunderstandings.

    if w1 * w_m0 != 2 * w_m2:

The (in)equality here should be modulo integer, namely if (w1*w_m0 - 2*w_m2) % integer != 0.

        return False
    b = pow(b, (integer - 1) // 2, integer)
    return b * w_m0 % integer == 2

Note, however, that if n is a prime, then

b^((n-1)/2) ≡ (b/n) (mod n)

where (b/n) is the Legendre (or Jacobi) symbol (for prime 'denominators', the Jacobi symbol is the Legendre symbol), hence b^((n-1)/2) ≡ ±1 (mod n). So you could use that as an extra check, if W_m is not 2 or n-2, n can't be prime, nor can it be if b^((n-1)/2) (mod n) is not 1 or n-1.

Probably computing b^((n-1)/2) (mod n) first and checking whether that's 1 or n-1 is a good idea, since if that check fails (that is the Euler pseudoprime test, by the way) you don't need the other, no less expensive, computations anymore, and if it succeeds, it's very likely that you need to compute it anyway.

Regarding the corrections, they seem correct, except for one that made a glitch I previously overlooked possibly worse:

if w1 * wm != 2 * wm1 % integer:

That applies the modulus only to 2 * wm1.

Concerning the recursion for the W_j, I think it is best to explain with a working implementation, first in toto for easy copy and paste:

def compute_wm_wm1(w1,m,n):
    a, b = 2, w1
    bits = int(log(m,2)) - 2
    if bits < 0:
        bits = 0
    mask = 1 << bits
    while mask <= m:
        mask <<= 1
    mask >>= 1
    while mask > 0:
        if (mask & m) != 0:
            a, b = (a*b-w1)%n, (b*b-2)%n
        else:
            a, b = (a*a-2)%n, (a*b-w1)%n
        mask >>= 1
    return a, b

Then with explanations in between:

def compute_wm_wm1(w1,m,n):

We need the value of W₁, the index of the desired number, and the number by which to take the modulus as input. The value W₀ is always 2, so we don't need that as a parameter.

Call it as

wm, wm1 = compute_wm_wm1(w1,m,integer)

in frobenius_pseudoprime (aside: not a good name, most of the numbers returning True are real primes).

    a, b = 2, w1

We initialise a and b to W₀ and W₁ respectively. At each point, a holds the value of W_j and b the value of W_j+1, where j is the value of the bits of m so far consumed. For example, with m = 13, the values of j, a and b develop as follows:

consumed remaining  j    a    b
           1101     0   w_0  w_1
   1        101     1   w_1  w_2
   11        01     3   w_3  w_4
   110        1     6   w_6  w_7
   1101            13  w_13  w_14

The bits are consumed left-to-right, so we have to find the first set bit of m and place our 'pointer' right before it

    bits = int(log(m,2)) - 2
    if bits < 0:
        bits = 0
    mask = 1 << bits

I subtracted a bit from the computed logarithm just to be entirely sure that we don't get fooled by a floating point error (by the way, using log limits you to numbers of at most 1024 bits, about 308 decimal digits; if you want to treat larger numbers, you have to find the base-2 logarithm of m in a different way, using log was the simplest way, and it's just a proof of concept, so I used that here).

    while mask <= m:
        mask <<= 1

Shift the mask until it's greater than m,so the set bit points just before m's first set bit. Then shift one position back, so we point at the bit.

    mask >>= 1
    while mask > 0:
        if (mask & m) != 0:
            a, b = (a*b-w1)%n, (b*b-2)%n

If the next bit is set, the value of the initial portion of consumed bits of m goes from j to 2*j+1, so the next values of the W sequence we need are W_2j+1 for a and W_2j+2 for b. By the above recursion formula,

W_{2j+1} = W_j * W_{j+1} - W_1 (mod n)
W_{2j+2} = W_{j+1}^2 - 2 (mod n)

Since a was W_j and b was W_j+1, a becomes (a*b - W_1) % n and b becomes (b * b - 2) % n.

        else:
            a, b = (a*a-2)%n, (a*b-w1)%n

If the next bit is not set, the value of the initial portion of consumed bits of m goes from j to 2*j, so a becomes W_2j = (W_j² - 2) (mod n), and b becomes W_2j+1 = (W_j * W_j+1 - W₁) (mod n).

        mask >>= 1

Move the pointer to the next bit. When we have moved past the final bit, mask becomes 0 and the loop ends. The initial portion of consumed bits of m is now all of m's bits, so the value is of course m. Then we can

    return a, b

Some additional remarks:

def find_prime_number(bits, test):
    while True:
        number = random(3, 1 << bits, 2)
        for _ in range(test):
            if not frobenius_pseudoprime(number):
                break
        else:
            return number

Primes are not too frequent among the larger numbers, so just picking random numbers is likely to take a lot of attempts to hit one. You will probably find a prime (or probable prime) faster if you pick one random number and check candidates in order.

Another point is that such a test as the Frobenius test is disproportionally expensive to find that e.g. a multiple of 3 is composite. Before using such a test (or a Miller-Rabin test, or a Lucas test, or an Euler test, ...), you should definitely do a bit of trial division to weed out most of the composites and do the work only where it has a fighting chance of being worth it.

Oh, and the is_square function isn't prepared to deal with arguments less than 2, divide-by-zero errors lurk there,

def is_square(integer):
    if integer < 0:
        return False
    if integer < 2:
        return True
    x = integer // 2

should help.