我已经尝试实施Baillie-PSW primality test几天了，但遇到了一些问题。尝试使用Lucas probable prime test时特别。我的问题不是关于Baile，而是关于如何以一些数为模生成正确的Lucas序列
对于前两个伪素，我的代码给出了正确的结果，例如对于323和377。但是，对于下一个psudoprime，标准实现和加倍版本都将失败。
尝试在V_1上进行模运算会完全破坏Luckas序列生成器的加倍版本。
关于如何正确地在Python中实现Lucas可能的主要测试的任何提示或建议？
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:
return jacobi(a % m, m)
elif a % 4 == 3 and m % 4 == 3:
return -jacobi(m, a)
return jacobi(m, a)
def smallest_D(num):
D = 5
k = 1
while k > 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)
最佳答案
这是我的卢卡斯伪素测验；您可以在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)
请注意，1159是已知的Lucas伪素数（A217120）。
关于python - 卢卡斯可能的主要测试，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/38343738/