Prueba probable de Lucas

He estado tratando de implementar la prueba de primalidad Baillie-PSW durante unos días y he tenido algunos problemas. Sepcíficamente al tratar de usar la prueba probable de Lucas . Mi pregunta no es sobre Baile, sino sobre cómo generar la secuencia de Lucas correcta.

Para los dos primeros psudoprimes mi código da el resultado correcto, por ejemplo, para 323 y 377 . Sin embargo, para el próximo psudoprime, tanto la implementación estándar como la versión de duplicación fallan.

Intentar realizar operaciones de módulo en V_1 rompe por completo la versión duplicada del generador de secuencias de Luckas.

¿Algún consejo o sugerencia sobre cómo implementar correctamente la prueba Prime probable de Lucas en Python?

 from fractions import gcd from math import log def luckas_sequence_standard(num, D=0): if D == 0: D = smallest_D(num) P = 1 Q = (1-D)/4 V0 = 2 V1 = P U0 = 0 U1 = 1 for _ in range(num): U2 = (P*U1 - Q*U0) % num U1, U0 = U2, U1 V2 = (P*V1 - Q*V0) % num V1, V0 = V2, V1 return U2%num, V2%num def luckas_sequence_doubling(num, D=0): if D == 0: D = smallest_D(num) P = 1 Q = (1 - D)/4 V0 = P U0 = 1 temp_num = num + 1 double = [] while temp_num > 1: if temp_num % 2 == 0: double.append(True) temp_num //= 2 else: double.append(False) temp_num += -1 k = 1 double.reverse() for is_double in double: if is_double: U1 = (U0*V0) % num V1 = V0**2 - 2*Q**k U0 = U1 V0 = V1 k *= 2 elif not is_double: U1 = ((P*U0 + V0)/2) % num V1 = (D*U0 + P*V0)/2 U0 = U1 V0 = V1 k += 1 return U1%num, V1%num def jacobi(a, m): if a in [0, 1]: return a elif gcd(a, m) != 1: return 0 elif a == 2: if m % 8 in [3, 5]: return -1 elif m % 8 in [1, 7]: return 1 if a % 2 == 0: return jacobi(2,m)*jacobi(a/2, m) elif a >= m or a  0 and jacobi(k*D, num) != -1: D += 2 k *= -1 return k*D if __name__ == '__main__': print luckas_sequence_standard(323) print luckas_sequence_doubling(323) print print luckas_sequence_standard(377) print luckas_sequence_doubling(377) print print luckas_sequence_standard(1159) print luckas_sequence_doubling(1159) 

Aquí está mi prueba de pseudoprimality de Lucas; Puede ejecutarlo en ideone.com/57Iayq .

 # lucas pseudoprimality test def gcd(a,b): # euclid's algorithm if b == 0: return a return gcd(b, a%b) def jacobi(a, m): # assumes a an integer and # m an odd positive integer a, t = a % m, 1 while a <> 0: z = -1 if m % 8 in [3,5] else 1 while a % 2 == 0: a, t = a / 2, t * z if a%4 == 3 and m%4 == 3: t = -t a, m = m % a, a return t if m == 1 else 0 def selfridge(n): d, s = 5, 1 while True: ds = d * s if gcd(ds, n) > 1: return ds, 0, 0 if jacobi(ds, n) == -1: return ds, 1, (1 - ds) / 4 d, s = d + 2, s * -1 def lucasPQ(p, q, m, n): # nth element of lucas sequence with # parameters p and q (mod m); ignore # modulus operation when m is zero def mod(x): if m == 0: return x return x % m def half(x): if x % 2 == 1: x = x + m return mod(x / 2) un, vn, qn = 1, p, q u = 0 if n % 2 == 0 else 1 v = 2 if n % 2 == 0 else p k = 1 if n % 2 == 0 else q n, d = n // 2, p * p - 4 * q while n > 0: u2 = mod(un * vn) v2 = mod(vn * vn - 2 * qn) q2 = mod(qn * qn) n2 = n // 2 if n % 2 == 1: uu = half(u * v2 + u2 * v) vv = half(v * v2 + d * u * u2) u, v, k = uu, vv, k * q2 un, vn, qn, n = u2, v2, q2, n2 return u, v, k def isLucasPseudoprime(n): d, p, q = selfridge(n) if p == 0: return n == d u, v, k = lucasPQ(p, q, n, n+1) return u == 0 print isLucasPseudoprime(1159) 

Tenga en cuenta que 1159 es una pseudoprima de Lucas conocida ( A217120 ).