fork(11) download 
  1. # lucas pseudoprimality test
  2.  
  3. def gcd(a,b): # euclid's algorithm
  4.     if b == 0: return a
  5.     return gcd(b, a%b)
  6.  
  7. def jacobi(a, m):
  8.     # assumes a an integer and
  9.     # m an odd positive integer
  10.     a, t = a % m, 1
  11.     while a <> 0:
  12.         z = -1 if m % 8 in [3,5] else 1
  13.         while a % 2 == 0:
  14.             a, t = a / 2, t * z
  15.         if a%4 == 3 and m%4 == 3: t = -t
  16.         a, m = m % a, a
  17.     return t if m == 1 else 0
  18.  
  19. def selfridge(n):
  20.     d, s = 5, 1
  21.     while True:
  22.         ds = d * s
  23.         if gcd(ds, n) > 1:
  24.             return ds, 0, 0
  25.         if jacobi(ds, n) == -1:
  26.             return ds, 1, (1 - ds) / 4
  27.         d, s = d + 2, s * -1
  28.  
  29. def lucasPQ(p, q, m, n):
  30.     # nth element of lucas sequence with
  31.     # parameters p and q (mod m); ignore
  32.     # modulus operation when m is zero
  33.     def mod(x):
  34.         if m == 0: return x
  35.         return x % m
  36.     def half(x):
  37.         if x % 2 == 1: x = x + m
  38.         return mod(x / 2)
  39.     un, vn, qn = 1, p, q
  40.     u = 0 if n % 2 == 0 else 1
  41.     v = 2 if n % 2 == 0 else p
  42.     k = 1 if n % 2 == 0 else q
  43.     n, d = n // 2, p * p - 4 * q
  44.     while n > 0:
  45.         u2 = mod(un * vn)
  46.         v2 = mod(vn * vn - 2 * qn)
  47.         q2 = mod(qn * qn)
  48.         n2 = n // 2
  49.         if n % 2 == 1:
  50.             uu = half(u * v2 + u2 * v)
  51.             vv = half(v * v2 + d * u * u2)
  52.             u, v, k = uu, vv, k * q2
  53.         un, vn, qn, n = u2, v2, q2, n2
  54.     return u, v, k
  55.  
  56. def isLucasPseudoprime(n):
  57.     d, p, q = selfridge(n)
  58.     if p == 0: return n == d
  59.     u, v, k = lucasPQ(p, q, n, n+1)
  60.     return u == 0
  61.  
  62. print isLucasPseudoprime(1159)
Success #stdin #stdout 0.01s 9016KB
comments (?)
stdin 
Standard input is empty
stdout 
True
https://ideone.com/57Iayq
language:
Python (cpython 2.7.16)
created:
9 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!