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