fork download
copy 
  1. # lenstra's algorithm per bosma/lenstra
  2.  
  3. from random import randint
  4. from fractions import gcd
  5.  
  6. def primes(n):
  7.     b, p, ps = [True] * (n+1), 2, []
  8.     for p in xrange(2, n+1):
  9.         if b[p]:
  10.             ps.append(p)
  11.             for i in xrange(p, n+1, p):
  12.                 b[i] = False
  13.     return ps
  14.  
  15. def bezout(a, b):
  16.     if b == 0: return 1, 0, a
  17.     q, r = divmod(a, b)
  18.     x, y, g = bezout(b, r)
  19.     return y, x-q*y, g
  20.  
  21. def add(p, q, a, b, m):
  22.     if p[2] == 0: return q
  23.     if q[2] == 0: return p
  24.     if p[0] == q[0]:
  25.         if (p[1] + q[1]) % m == 0:
  26.             return 0, 1, 0 # infinity
  27.         n = (3 * p[0] * p[0] + a) % m
  28.         d = (2 * p[1]) % m
  29.     else:
  30.         n = (q[1] - p[1]) % m
  31.         d = (q[0] - p[0]) % m
  32.     x, y, g = bezout(d, m)
  33.     if g > 1: return 0, 0, d # failure
  34.     z = (n*x*n*x - p[0] - q[0]) % m
  35.     return z, (n * x * (p[0] - z) - p[1]) % m, 1
  36.  
  37. def mul(k, p, a, b, m):
  38.     r = (0,1,0)
  39.     while k > 0:
  40.         if k % 2 == 1:
  41.             r = add(p, r, a, b, m)
  42.             if r[2] > 1: return r
  43.         k = k // 2
  44.         p = add(p, p, a, b, m)
  45.         if p[2] > 1: return p
  46.     return r
  47.  
  48. def lenstra1(n, limit):
  49.     g = n
  50.     while g == n:
  51.         q = randint(0, n-1), randint(0, n-1), 1
  52.         a = randint(0, n-1)
  53.         b = (q[1]*q[1] - q[0]*q[0]*q[0] - a*q[0]) % n
  54.         g = gcd(4*a*a*a + 27*b*b, n)
  55.     if g > 1: return 0, g # lucky factor
  56.     for p in primes(limit):
  57.         pp = p
  58.         while pp < limit:
  59.             q = mul(p, q, a, b, n)
  60.             if q[2] > 1:
  61.                 return 1, gcd(q[2], n)
  62.             pp = p * pp
  63.     return False
  64.  
  65. def parms(b1):
  66.     b2 = 10 * b1
  67.     er = [(1,31), (2,63), (3,127), (6,255), (12,511), \
  68.           (18,511), (24,1023), (30,1023), (60,2047)]
  69.     prev = 1,31
  70.     for (e, r) in er:
  71.         if e*e > b1/1250: break
  72.         prev = e, r
  73.     e, r = prev
  74.     rBar = int(round(b2/r))
  75.     u = randint(0, pow(2,30)//(e+2))
  76.     v = randint(0, pow(2,30)//(e+2))
  77.     uBar = randint(0, pow(2,30)//(e+2))
  78.     vBar = randint(0, pow(2,30)//(e+2))
  79.     return b2, e, r, rBar, u, v, uBar, vBar
  80.  
  81. def lenstra2(n, b1):
  82.     g = n
  83.     while g == n:
  84.         q = randint(0, n-1), randint(0, n-1), 1
  85.         a = randint(0, n-1)
  86.         b = (q[1]*q[1] - q[0]*q[0]*q[0] - a*q[0]) % n
  87.         g = gcd(4*a*a*a + 27*b*b, n)
  88.     if g > 1: return 0, g # lucky factor
  89.     for p in primes(b1):
  90.         pp = p
  91.         while pp < b1:
  92.             q = mul(p, q, a, b, n)
  93.             if q[2] > 1: return 1, gcd(q[2], n)
  94.             pp = p * pp
  95.     b2, e, r, rBar, u, v, uBar, vBar = parms(b1)
  96.     f = [1] * (r+1)
  97.     for i in range(1, r):
  98.         p = mul(pow(u*i+v,e), q, a, b, n)
  99.         if p[2] > 1: return 2, gcd(p[2], n)
  100.         f[i] = (f[i-1] * (q[0] - p[0])) % n
  101.     d = 1
  102.     for j in range(1, rBar):
  103.         pBar = mul(pow(uBar*j+vBar,e), q, a, b, n)
  104.         if pBar[2] > 1: return 3, gcd(pBar[2], n)
  105.         t = 0
  106.         for i in range(0, r):
  107.             t = (t + p[0] * f[i]) % n
  108.         d = (d * t) % n
  109.     g = gcd(d, n)
  110.     if 1 < g < n: return 4, g
  111.     return False
  112.  
  113. print lenstra1(4636028740351, 500)
  114. print lenstra2(4636028740351, 500)
  115. print lenstra2(28021516895802718283942993, 1000)
Success #stdin #stdout 1.77s 12792KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
False
(3, 9387467L)
False
https://ideone.com/vMXSab
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!