You seem to misunderstand the algorithm due to ignorance of the notation.
def frobenius_pseudoprime(integer): assert integer & 1 and integer >= 3 a, b, d = choose_ab(integer) w1 = (a ** 2 // b - 2) % integer
It comes from a line
W 0 ≡ 2 (mod n) and W 1 ≡ a 2 b -1 - 2 (mod n)
But b -1 does not mean 1/b , but the modular inverse of b modulo n , i.e. integer c with b·c ≡ 1 (mod n) . You can most easily find such a c by continuing the fractional expansion of b/n or, which is the same, but with a bit more calculations, using the extended Euclidean algorithm. Since you are probably not familiar with ongoing fractions, I recommend the latter.
m = (integer - d // integer) // 2
derived from
n - (Δ / n) = 2m
and misunderstands the Jacobi symbol as a fraction / division (admittedly, I showed it here even more as a fraction, but since the site does not support LaTeX rendering, we will need to do it). The Jacobi symbol is a generalization of the Legendre symbol, denoted identically, which indicates whether the number is a quadratic residue modulo an odd prime number (if n is a quadratic residue modulo p , i.e. there exists k with k^2 ≡ n (mod p) and n 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 removes the restriction that the “denominator” is an odd prime and allows arbitrary odd numbers to be “denominators”. Its value is a product of Legendre symbols with the same “numerator” for all primes dividing n (in multiplicity). Learn more about this and how to efficiently calculate Jacobi characters in a related article. The string must read correctly
m = (integer - jacobi_symbol(d,integer)) // 2
The following lines, which I do not fully understand, logically, should follow the calculations of W m and W m + 1 using recursion
W 2j ≡ W j 2 - 2 (mod n)
W 2j + 1 ≡ W j W j + 1 - W 1 (mod n)
An effective method of using this recursion to calculate 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 rest of the function is almost correct, except that now it receives incorrect data due to earlier misunderstandings.
if w1 * w_m0 != 2 * w_m2:
Here the equality (in) 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 simple "denominators", the Jacobi symbol is the Legendre symbol), therefore b^((n-1)/2) ≡ ±1 (mod n) . Thus, you can use this as an additional check if W m is not 2 or n-2 , n cannot be simple, and cannot be if b^((n-1)/2) (mod n) not 1 or n-1 .
Probably by first calculating b^((n-1)/2) (mod n) and checking if this 1 or n-1 good idea, because if this test fails (this is a test of the Euler pseudo-terminal , by the way) you we need other, no less expensive calculations, and if it succeeds, it is very likely that you still need to calculate it.
As for the corrections, they seem to be correct, except for what made the failure, which I previously overlooked, possibly worse:
if w1 * wm != 2 * wm1 % integer:
This module applies only to 2 * wm1 .
Regarding recursion for W j , I think that 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 between:
def compute_wm_wm1(w1,m,n):
We need the value of W 1 , the index of the required number, and the number with which the module should be entered as input. The value of W 0 is always equal to 2, so we do not need this as a parameter.
Name him
wm, wm1 = compute_wm_wm1(w1,m,integer)
in frobenius_pseudoprime (aside: not a good name, most of the numbers that return True are real touches).
a, b = 2, w1
Initialize a and b on W 0 and W 1, respectively. At each point, a contains the value of W j and b value of W j + 1 , where j is the value of the bits m are still consumed. For example, with m = 13 values of j , a and b develop as follows:
consumed remaining jab 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
Bits are consumed from left to right, so we need to find the first bit of the set m and put our "pointer" right in front of it
bits = int(log(m,2)) - 2 if bits < 0: bits = 0 mask = 1 << bits
I subtracted a bit from the calculated logarithm just to be sure that we are not cheating with a floating point error (by the way, using log limits you to no more than 1024 bits, about 308 decimal digits, if you want to process large numbers, you need finding the logarithm of base-2 m differently, using log was the easiest way and this is just a proof of concept, so I used it here)
while mask <= m: mask <<= 1
Move the mask until it becomes larger than m , therefore, the bit of the set value before the very value of m will be set first. Then slide one position back so that we point to a 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 part of the consumed bits m goes from j to 2*j+1 , so the next values of the sequence W that we need are W 2j + 1 for a and W 2j + 2 for b . According to the recursion formula above
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 part of the consumed bits m goes from j to 2*j , so a becomes W 2j = (W j 2 - 2) (mod n), and b becomes W 2j + 1 = (W j * W j + 1 - W 1 ) (mod n).
mask >>= 1
Move the pointer to the next bit. When we have passed the last bit, mask becomes 0, and the loop ends. The initial part of the consumed bits m now represents the entire bit m , so the value, of course, is m . Then we can
return a, b
Some additional notes:
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
Prime numbers are not too frequent among large numbers, so just picking random numbers is likely to take many attempts to hit them. You are likely to find prime (or likely prime) if you pick one random number and check the candidates in order.
Another point is that a test such as the Frobenius test is disproportionately expensive to find, for example, a multiple of 3 - compound. Before using such a test (or the Miller-Rabin test, or the Lucas test, or the Euler test, ...), you should definitely do a little test division to weed out most of the composites and only do work where he has a chance to fight him.
Oh, and the is_square function is_square not ready to handle arguments less than 2, there are div-by-zero lurk there errors,
def is_square(integer): if integer < 0: return False if integer < 2: return True x = integer // 2
should help.