fork(1) download
copy 
  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.         print ds
  24.         if gcd(ds, n) > 1:
  25.             print "ALPHA"
  26.             return ds, 0, 0
  27.         if jacobi(ds, n) == -1:
  28.             print "BETA"
  29.             return ds, 1, (1 - ds) / 4
  30.         d, s = d + 2, s * -1
  31.  
  32. def lucasPQ(p, q, m, n):
  33.     # nth element of lucas sequence with
  34.     # parameters p and q (mod m); ignore
  35.     # modulus operation when m is zero
  36.     def mod(x):
  37.         if m == 0: return x
  38.         return x % m
  39.     def half(x):
  40.         if x % 2 == 1: x = x + m
  41.         return mod(x / 2)
  42.     un, vn, qn = 1, p, q
  43.     u = 0 if n % 2 == 0 else 1
  44.     v = 2 if n % 2 == 0 else p
  45.     k = 1 if n % 2 == 0 else q
  46.     n, d = n // 2, p * p - 4 * q
  47.     while n > 0:
  48.         u2 = mod(un * vn)
  49.         v2 = mod(vn * vn - 2 * qn)
  50.         q2 = mod(qn * qn)
  51.         n2 = n // 2
  52.         if n % 2 == 1:
  53.             uu = half(u * v2 + u2 * v)
  54.             vv = half(v * v2 + d * u * u2)
  55.             u, v, k = uu, vv, k * q2
  56.         un, vn, qn, n = u2, v2, q2, n2
  57.     return u, v, k
  58.  
  59. def isLucasPseudoprime(n):
  60.     d, p, q = selfridge(n)
  61.     if p == 0: return n == d
  62.     u, v, k = lucasPQ(p, q, n, n+1)
  63.     return u == 0
  64.  
  65. print isLucasPseudoprime(211)
Success #stdin #stdout 0.01s 7080KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
5
-7
9
-11
BETA
True
https://ideone.com/HyiKuj
language:
Python (cpython 2.7.16)
created:
6 hours 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!