fork(5) download 
  1. # lenstra's algorithm
  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 lenstra(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 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 gcd(q[2], n)
  62.             pp = p * pp
  63.     return False
  64.  
  65. # 523 * 769
  66. print lenstra(402187, 25)
  67. print lenstra(402187, 25)
  68. print lenstra(402187, 25)
  69. print lenstra(402187, 25)
  70. print lenstra(402187, 25)
  71. print lenstra(402187, 25)
  72. print lenstra(402187, 25)
  73. print lenstra(402187, 25)
  74. print lenstra(402187, 25)
  75. print lenstra(402187, 25)
  76.  
  77. print ""
  78.  
  79. # 4723 * 82193
  80. print lenstra(388197539, 100)
  81. print lenstra(388197539, 100)
  82. print lenstra(388197539, 100)
  83. print lenstra(388197539, 100)
  84. print lenstra(388197539, 100)
  85. print lenstra(388197539, 100)
  86. print lenstra(388197539, 100)
  87. print lenstra(388197539, 100)
  88. print lenstra(388197539, 100)
  89. print lenstra(388197539, 100)
  90.  
  91. print ""
  92.  
  93. # 678413 * 49153471
  94. print lenstra(33346353721523, 1000)
  95. print lenstra(33346353721523, 1000)
  96. print lenstra(33346353721523, 1000)
  97. print lenstra(33346353721523, 1000)
  98. print lenstra(33346353721523, 1000)
  99. print lenstra(33346353721523, 1000)
  100. print lenstra(33346353721523, 1000)
  101. print lenstra(33346353721523, 1000)
  102. print lenstra(33346353721523, 1000)
  103. print lenstra(33346353721523, 1000)
Success #stdin #stdout 0.76s 12784KB
comments (?)
stdin 
Standard input is empty
stdout 
False
523
False
523
False
769
523
False
769
523

4723
82193
False
4723
4723
False
False
4723
4723
False

False
678413
678413
678413
False
49153471
False
678413
678413
False
https://ideone.com/QEDmMY
language:
Python (cpython 2.7.16)
created:
12 years ago
visibility:
secret

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!