# lenstra's algorithm per bosma/lenstra
from random import randint
from fractions import gcd
def primes(n):
b, p, ps = [True] * (n+1), 2, []
for p in xrange(2, n+1):
if b[p]:
ps.append(p)
for i in xrange(p, n+1, p):
b[i] = False
return ps
def bezout(a, b):
if b == 0: return 1, 0, a
q, r = divmod(a, b)
x, y, g = bezout(b, r)
return y, x-q*y, g
def add(p, q, a, b, m):
if p[2] == 0: return q
if q[2] == 0: return p
if p[0] == q[0]:
if (p[1] + q[1]) % m == 0:
return 0, 1, 0 # infinity
n = (3 * p[0] * p[0] + a) % m
d = (2 * p[1]) % m
else:
n = (q[1] - p[1]) % m
d = (q[0] - p[0]) % m
x, y, g = bezout(d, m)
if g > 1: return 0, 0, d # failure
z = (n*x*n*x - p[0] - q[0]) % m
return z, (n * x * (p[0] - z) - p[1]) % m, 1
def mul(k, p, a, b, m):
r = (0,1,0)
while k > 0:
if k % 2 == 1:
r = add(p, r, a, b, m)
if r[2] > 1: return r
k = k // 2
p = add(p, p, a, b, m)
if p[2] > 1: return p
return r
def lenstra1(n, limit):
g = n
while g == n:
q = randint(0, n-1), randint(0, n-1), 1
a = randint(0, n-1)
b = (q[1]*q[1] - q[0]*q[0]*q[0] - a*q[0]) % n
g = gcd(4*a*a*a + 27*b*b, n)
if g > 1: return 0, g # lucky factor
for p in primes(limit):
pp = p
while pp < limit:
q = mul(p, q, a, b, n)
if q[2] > 1:
return 1, gcd(q[2], n)
pp = p * pp
return False
def parms(b1):
b2 = 10 * b1
er = [(1,31), (2,63), (3,127), (6,255), (12,511), \
(18,511), (24,1023), (30,1023), (60,2047)]
prev = 1,31
for (e, r) in er:
if e*e > b1/1250: break
prev = e, r
e, r = prev
rBar = int(round(b2/r))
u = randint(0, pow(2,30)//(e+2))
v = randint(0, pow(2,30)//(e+2))
uBar = randint(0, pow(2,30)//(e+2))
vBar = randint(0, pow(2,30)//(e+2))
return b2, e, r, rBar, u, v, uBar, vBar
def lenstra2(n, b1):
g = n
while g == n:
q = randint(0, n-1), randint(0, n-1), 1
a = randint(0, n-1)
b = (q[1]*q[1] - q[0]*q[0]*q[0] - a*q[0]) % n
g = gcd(4*a*a*a + 27*b*b, n)
if g > 1: return 0, g # lucky factor
for p in primes(b1):
pp = p
while pp < b1:
q = mul(p, q, a, b, n)
if q[2] > 1: return 1, gcd(q[2], n)
pp = p * pp
b2, e, r, rBar, u, v, uBar, vBar = parms(b1)
f = [1] * (r+1)
for i in range(1, r):
p = mul(pow(u*i+v,e), q, a, b, n)
if p[2] > 1: return 2, gcd(p[2], n)
f[i] = (f[i-1] * (q[0] - p[0])) % n
d = 1
for j in range(1, rBar):
pBar = mul(pow(uBar*j+vBar,e), q, a, b, n)
if pBar[2] > 1: return 3, gcd(pBar[2], n)
t = 0
for i in range(0, r):
t = (t + p[0] * f[i]) % n
d = (d * t) % n
g = gcd(d, n)
if 1 < g < n: return 4, g
return False
print lenstra1(4636028740351, 500)
print lenstra2(4636028740351, 500)
print lenstra2(28021516895802718283942993, 1000)
IyBsZW5zdHJhJ3MgYWxnb3JpdGhtIHBlciBib3NtYS9sZW5zdHJhCgpmcm9tIHJhbmRvbSBpbXBvcnQgcmFuZGludApmcm9tIGZyYWN0aW9ucyBpbXBvcnQgZ2NkCgpkZWYgcHJpbWVzKG4pOgogICAgYiwgcCwgcHMgPSBbVHJ1ZV0gKiAobisxKSwgMiwgW10KICAgIGZvciBwIGluIHhyYW5nZSgyLCBuKzEpOgogICAgICAgIGlmIGJbcF06CiAgICAgICAgICAgIHBzLmFwcGVuZChwKQogICAgICAgICAgICBmb3IgaSBpbiB4cmFuZ2UocCwgbisxLCBwKToKICAgICAgICAgICAgICAgIGJbaV0gPSBGYWxzZQogICAgcmV0dXJuIHBzCgpkZWYgYmV6b3V0KGEsIGIpOgogICAgaWYgYiA9PSAwOiByZXR1cm4gMSwgMCwgYQogICAgcSwgciA9IGRpdm1vZChhLCBiKQogICAgeCwgeSwgZyA9IGJlem91dChiLCByKQogICAgcmV0dXJuIHksIHgtcSp5LCBnCgpkZWYgYWRkKHAsIHEsIGEsIGIsIG0pOgogICAgaWYgcFsyXSA9PSAwOiByZXR1cm4gcQogICAgaWYgcVsyXSA9PSAwOiByZXR1cm4gcAogICAgaWYgcFswXSA9PSBxWzBdOgogICAgICAgIGlmIChwWzFdICsgcVsxXSkgJSBtID09IDA6CiAgICAgICAgICAgIHJldHVybiAwLCAxLCAwICMgaW5maW5pdHkKICAgICAgICBuID0gKDMgKiBwWzBdICogcFswXSArIGEpICUgbQogICAgICAgIGQgPSAoMiAqIHBbMV0pICUgbQogICAgZWxzZToKICAgICAgICBuID0gKHFbMV0gLSBwWzFdKSAlIG0KICAgICAgICBkID0gKHFbMF0gLSBwWzBdKSAlIG0KICAgIHgsIHksIGcgPSBiZXpvdXQoZCwgbSkKICAgIGlmIGcgPiAxOiByZXR1cm4gMCwgMCwgZCAjIGZhaWx1cmUKICAgIHogPSAobip4Km4qeCAtIHBbMF0gLSBxWzBdKSAlIG0KICAgIHJldHVybiB6LCAobiAqIHggKiAocFswXSAtIHopIC0gcFsxXSkgJSBtLCAxCgpkZWYgbXVsKGssIHAsIGEsIGIsIG0pOgogICAgciA9ICgwLDEsMCkKICAgIHdoaWxlIGsgPiAwOgogICAgICAgIGlmIGsgJSAyID09IDE6CiAgICAgICAgICAgIHIgPSBhZGQocCwgciwgYSwgYiwgbSkKICAgICAgICAgICAgaWYgclsyXSA+IDE6IHJldHVybiByCiAgICAgICAgayA9IGsgLy8gMgogICAgICAgIHAgPSBhZGQocCwgcCwgYSwgYiwgbSkKICAgICAgICBpZiBwWzJdID4gMTogcmV0dXJuIHAKICAgIHJldHVybiByCgpkZWYgbGVuc3RyYTEobiwgbGltaXQpOgogICAgZyA9IG4KICAgIHdoaWxlIGcgPT0gbjoKICAgICAgICBxID0gcmFuZGludCgwLCBuLTEpLCByYW5kaW50KDAsIG4tMSksIDEKICAgICAgICBhID0gcmFuZGludCgwLCBuLTEpCiAgICAgICAgYiA9IChxWzFdKnFbMV0gLSBxWzBdKnFbMF0qcVswXSAtIGEqcVswXSkgJSBuCiAgICAgICAgZyA9IGdjZCg0KmEqYSphICsgMjcqYipiLCBuKQogICAgaWYgZyA+IDE6IHJldHVybiAwLCBnICMgbHVja3kgZmFjdG9yCiAgICBmb3IgcCBpbiBwcmltZXMobGltaXQpOgogICAgICAgIHBwID0gcAogICAgICAgIHdoaWxlIHBwIDwgbGltaXQ6CiAgICAgICAgICAgIHEgPSBtdWwocCwgcSwgYSwgYiwgbikKICAgICAgICAgICAgaWYgcVsyXSA+IDE6CiAgICAgICAgICAgICAgICByZXR1cm4gMSwgZ2NkKHFbMl0sIG4pCiAgICAgICAgICAgIHBwID0gcCAqIHBwCiAgICByZXR1cm4gRmFsc2UKCmRlZiBwYXJtcyhiMSk6CiAgICBiMiA9IDEwICogYjEKICAgIGVyID0gWygxLDMxKSwgKDIsNjMpLCAoMywxMjcpLCAoNiwyNTUpLCAoMTIsNTExKSwgXAogICAgICAgICAgKDE4LDUxMSksICgyNCwxMDIzKSwgKDMwLDEwMjMpLCAoNjAsMjA0NyldCiAgICBwcmV2ID0gMSwzMQogICAgZm9yIChlLCByKSBpbiBlcjoKICAgICAgICBpZiBlKmUgPiBiMS8xMjUwOiBicmVhawogICAgICAgIHByZXYgPSBlLCByCiAgICBlLCByID0gcHJldgogICAgckJhciA9IGludChyb3VuZChiMi9yKSkKICAgIHUgPSByYW5kaW50KDAsIHBvdygyLDMwKS8vKGUrMikpCiAgICB2ID0gcmFuZGludCgwLCBwb3coMiwzMCkvLyhlKzIpKQogICAgdUJhciA9IHJhbmRpbnQoMCwgcG93KDIsMzApLy8oZSsyKSkKICAgIHZCYXIgPSByYW5kaW50KDAsIHBvdygyLDMwKS8vKGUrMikpCiAgICByZXR1cm4gYjIsIGUsIHIsIHJCYXIsIHUsIHYsIHVCYXIsIHZCYXIKCmRlZiBsZW5zdHJhMihuLCBiMSk6CiAgICBnID0gbgogICAgd2hpbGUgZyA9PSBuOgogICAgICAgIHEgPSByYW5kaW50KDAsIG4tMSksIHJhbmRpbnQoMCwgbi0xKSwgMQogICAgICAgIGEgPSByYW5kaW50KDAsIG4tMSkKICAgICAgICBiID0gKHFbMV0qcVsxXSAtIHFbMF0qcVswXSpxWzBdIC0gYSpxWzBdKSAlIG4KICAgICAgICBnID0gZ2NkKDQqYSphKmEgKyAyNypiKmIsIG4pCiAgICBpZiBnID4gMTogcmV0dXJuIDAsIGcgIyBsdWNreSBmYWN0b3IKICAgIGZvciBwIGluIHByaW1lcyhiMSk6CiAgICAgICAgcHAgPSBwCiAgICAgICAgd2hpbGUgcHAgPCBiMToKICAgICAgICAgICAgcSA9IG11bChwLCBxLCBhLCBiLCBuKQogICAgICAgICAgICBpZiBxWzJdID4gMTogcmV0dXJuIDEsIGdjZChxWzJdLCBuKQogICAgICAgICAgICBwcCA9IHAgKiBwcAogICAgYjIsIGUsIHIsIHJCYXIsIHUsIHYsIHVCYXIsIHZCYXIgPSBwYXJtcyhiMSkKICAgIGYgPSBbMV0gKiAocisxKQogICAgZm9yIGkgaW4gcmFuZ2UoMSwgcik6CiAgICAgICAgcCA9IG11bChwb3codSppK3YsZSksIHEsIGEsIGIsIG4pCiAgICAgICAgaWYgcFsyXSA+IDE6IHJldHVybiAyLCBnY2QocFsyXSwgbikKICAgICAgICBmW2ldID0gKGZbaS0xXSAqIChxWzBdIC0gcFswXSkpICUgbgogICAgZCA9IDEKICAgIGZvciBqIGluIHJhbmdlKDEsIHJCYXIpOgogICAgICAgIHBCYXIgPSBtdWwocG93KHVCYXIqait2QmFyLGUpLCBxLCBhLCBiLCBuKQogICAgICAgIGlmIHBCYXJbMl0gPiAxOiByZXR1cm4gMywgZ2NkKHBCYXJbMl0sIG4pCiAgICAgICAgdCA9IDAKICAgICAgICBmb3IgaSBpbiByYW5nZSgwLCByKToKICAgICAgICAgICAgdCA9ICh0ICsgcFswXSAqIGZbaV0pICUgbgogICAgICAgIGQgPSAoZCAqIHQpICUgbgogICAgZyA9IGdjZChkLCBuKQogICAgaWYgMSA8IGcgPCBuOiByZXR1cm4gNCwgZwogICAgcmV0dXJuIEZhbHNlCgpwcmludCBsZW5zdHJhMSg0NjM2MDI4NzQwMzUxLCA1MDApCnByaW50IGxlbnN0cmEyKDQ2MzYwMjg3NDAzNTEsIDUwMCkKcHJpbnQgbGVuc3RyYTIoMjgwMjE1MTY4OTU4MDI3MTgyODM5NDI5OTMsIDEwMDAp