fork download
copy 
  1. # n+1 primality prover
  2.  
  3. def gcd(a, b):
  4.     while b <> 0:
  5.         a, b = b, a % b
  6.     return a
  7.  
  8. def signum(n):
  9.     if n > 0: return 1
  10.     if n < 0: return -1
  11.     return 0
  12.  
  13. def jacobi(a, p):
  14.     a, t = a % p, 1
  15.     while a <> 0:
  16.         while a % 2 == 0:
  17.             a /= 2
  18.             if p % 8 in (3, 5): t *= -1
  19.         a, p = p, a
  20.         if a % 4 == 3 and p % 4 == 3: t *= -1
  21.         a %= p
  22.     return t if p == 1 else 0
  23.  
  24. def u(p, q, m, n):
  25.     # nth element of Lucas U(p,q) sequence mod m
  26.     # for the bits of n from right to left:
  27.     # u2, v2, and q2 double at each bit and are added
  28.     # to the accumulators u, v, and q at each 1-bit
  29.     d, u2, v2, q2 = p*p - 4*q, 1, p, q
  30.     if n%2==1: u,v,q = 1,p,q
  31.     else: u,v,q = 0,2,1
  32.     while n > 1:
  33.         u2, v2, q2, n = \
  34.             (u2*v2)%m, (v2*v2-2*q2)%m, (q2*q2)%m, n//2
  35.         if n%2 == 1:
  36.             u, v = u*v2+u2*v, v*v2+d*u*u2
  37.             u = (u/2)%m if u%2==0 else ((u+m)/2)%m
  38.             v = (v/2)%m if v%2==0 else ((v+m)/2)%m
  39.             q *= q2
  40.     return u
  41.  
  42. def provePrime(n):
  43.     f, fs, fp, r, d = 2, [], 1, n+1, 5
  44.     while fp*fp < n:
  45.     	if f*f > r:
  46.     		fs.append(r)
  47.     		break
  48.         while r % f == 0:
  49.             r /= f; fp *= f
  50.             if f not in fs:
  51.                 fs.append(f)
  52.         f += 1
  53.     while jacobi(d, n) <> -1:
  54.         d = (abs(d)+2) * signum(d) * -1
  55.     p, q = 1, (1-d)/4
  56.     if gcd(d, n) > 1: return False
  57.     if u(p, q, n, n+1) <> 0: return False
  58.     for x in fs:
  59.         while True:
  60.             if gcd(u(p,q,n,(n+1)/x), n) == 1:
  61.             	print d, p, q, x, u(p, q, n, (n+1)/x)
  62.             	break
  63.             p, q = p + 2, p + q + 1
  64.     return True
  65.  
  66. print provePrime(10**12+39)
  67. print provePrime(10**12+61)
Success #stdin #stdout 0.04s 7896KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
-7 19 92 2 309271071084
-7 19 92 5 748491620867
-7 19 92 17573 18065186228
-7 19 92 1422637 418811399431
True
-7 1 2 2 989025798365
-7 15 58 3 968440455086
-7 15 58 7 176705259344
-7 15 58 47 429745757234
-7 15 58 168861871 566495358887
True
https://ideone.com/HVXIGR
language:
Python (cpython 2.7.16)
created:
11 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!