盧卡斯很可能的素數測試

Question

我一直試圖實施Baillie-PSW primality test幾天，並遇到了一些問題。 Sepcifically當試圖使用Lucas probable prime test。 我的問題是不是伯樂，但對如何生成正確的盧卡斯序列模一定數量

對於前兩種psudoprimes我的代碼給出正確的結果，例如用於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) 

Answer 1

這是我的盧卡斯僞原理測試;你可以運行它在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是已知的盧卡斯pseudoprime（A217120）。