Lucas wahrscheinlich prime test

Question

Ich habe versucht, die Baillie-PSW primality test für ein paar Tage zu implementieren, und habe einige Probleme festgestellt. Sepzifisch beim Versuch, das Lucas probable prime test zu verwenden. Meine Frage ist nicht über Baile, sondern darauf, wie die richtige Lucas-Sequenz zu erzeugen Modulo gewisse Anzahl

Für die ersten beiden psudoprimes meinen Code gibt das korrekte Ergebnis, zB für 323 und 377. Für die nächste psudoprime schlägt jedoch sowohl die Standardimplementierung als auch die Doubling-Version fehl.

Der Versuch, Modulo-Operationen an V_1 durchzuführen, bricht die verdoppelnde Version des Luckas-Sequenzgenerators vollständig.

Irgendwelche Tipps oder Vorschläge, wie Sie den Lucas-wahrscheinlichen Haupttest in Python korrekt implementieren?

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

Hier ist mein Lucas Pseudoprimalitätstest; Sie können es unter ideone.com/57Iayq ausführen.

# 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) 

Beachten Sie, dass 1159 ist eine bekannte Lucas pseudoprim (A217120).