def is_prime(n):
if type(n) != int and type(n) != long:
raise TypeError('must be integer')
if n < 2:
return False
ps = [2,3,5,7,11,13,17,19,23,29,31,37,41,
43,47,53,59,61,67,71,73,79,83,89,97]
def is_spsp(n, a):
d, s = n-1, 0
while d%2 == 0:
d /= 2; s += 1
t = pow(a,d,n)
if t == 1:
return True
while s > 0:
if t == n-1:
return True
t = (t*t) % n
s -= 1
return False
if n in ps: return True
for p in ps:
if not is_spsp(n,p):
return False
return True
def rho_factors(n, limit=100):
if type(n) != int and type(n) != long:
raise TypeError('must be integer')
def gcd(a,b):
while b: a, b = b, a%b
return abs(a)
def rho_factor(n, c, limit):
f = lambda(x): (x*x+c) % n
t, h, d = 2, 2, 1
while d == 1:
if limit == 0:
raise OverflowError('limit exceeded')
t = f(t); h = f(f(h)); d = gcd(t-h, n)
if d == n:
return rho_factor(n, c+1, limit)
if is_prime(d):
return d
return rho_factor(d, c+1, limit)
if -1 <= n <= 1: return [n]
if n < -1: return [-1] + rho_factors(-n, limit)
fs = []
while n % 2 == 0:
n = n // 2; fs = fs + [2]
if n == 1: return fs
while not is_prime(n):
f = rho_factor(n, 1, limit)
n = n / f
fs = fs + [f]
return sorted(fs + [n])
print rho_factors(99999640000243)
ZGVmIGlzX3ByaW1lKG4pOgogICAgaWYgdHlwZShuKSAhPSBpbnQgYW5kIHR5cGUobikgIT0gbG9uZzoKICAgICAgICByYWlzZSBUeXBlRXJyb3IoJ211c3QgYmUgaW50ZWdlcicpCiAgICBpZiBuIDwgMjoKICAgICAgICByZXR1cm4gRmFsc2UKICAgIHBzID0gWzIsMyw1LDcsMTEsMTMsMTcsMTksMjMsMjksMzEsMzcsNDEsCiAgICAgICAgIDQzLDQ3LDUzLDU5LDYxLDY3LDcxLDczLDc5LDgzLDg5LDk3XQogICAgZGVmIGlzX3Nwc3AobiwgYSk6CiAgICAgICAgZCwgcyA9IG4tMSwgMAogICAgICAgIHdoaWxlIGQlMiA9PSAwOgogICAgICAgICAgICBkIC89IDI7IHMgKz0gMQogICAgICAgIHQgPSBwb3coYSxkLG4pCiAgICAgICAgaWYgdCA9PSAxOgogICAgICAgICAgICByZXR1cm4gVHJ1ZQogICAgICAgIHdoaWxlIHMgPiAwOgogICAgICAgICAgICBpZiB0ID09IG4tMToKICAgICAgICAgICAgICAgIHJldHVybiBUcnVlCiAgICAgICAgICAgIHQgPSAodCp0KSAlIG4KICAgICAgICAgICAgcyAtPSAxCiAgICAgICAgcmV0dXJuIEZhbHNlCiAgICBpZiBuIGluIHBzOiByZXR1cm4gVHJ1ZQogICAgZm9yIHAgaW4gcHM6CiAgICAgICAgaWYgbm90IGlzX3Nwc3AobixwKToKICAgICAgICAgICAgcmV0dXJuIEZhbHNlCiAgICByZXR1cm4gVHJ1ZQoKZGVmIHJob19mYWN0b3JzKG4sIGxpbWl0PTEwMCk6CiAgICBpZiB0eXBlKG4pICE9IGludCBhbmQgdHlwZShuKSAhPSBsb25nOgogICAgICAgIHJhaXNlIFR5cGVFcnJvcignbXVzdCBiZSBpbnRlZ2VyJykKICAgIGRlZiBnY2QoYSxiKToKICAgICAgICB3aGlsZSBiOiBhLCBiID0gYiwgYSViCiAgICAgICAgcmV0dXJuIGFicyhhKQogICAgZGVmIHJob19mYWN0b3IobiwgYywgbGltaXQpOgogICAgICAgIGYgPSBsYW1iZGEoeCk6ICh4KngrYykgJSBuCiAgICAgICAgdCwgaCwgZCA9IDIsIDIsIDEKICAgICAgICB3aGlsZSBkID09IDE6CiAgICAgICAgICAgIGlmIGxpbWl0ID09IDA6CiAgICAgICAgICAgICAgICByYWlzZSBPdmVyZmxvd0Vycm9yKCdsaW1pdCBleGNlZWRlZCcpCiAgICAgICAgICAgIHQgPSBmKHQpOyBoID0gZihmKGgpKTsgZCA9IGdjZCh0LWgsIG4pCiAgICAgICAgaWYgZCA9PSBuOgogICAgICAgICAgICByZXR1cm4gcmhvX2ZhY3RvcihuLCBjKzEsIGxpbWl0KQogICAgICAgIGlmIGlzX3ByaW1lKGQpOgogICAgICAgICAgICByZXR1cm4gZAogICAgICAgIHJldHVybiByaG9fZmFjdG9yKGQsIGMrMSwgbGltaXQpCiAgICBpZiAtMSA8PSBuIDw9IDE6IHJldHVybiBbbl0KICAgIGlmIG4gPCAtMTogcmV0dXJuIFstMV0gKyByaG9fZmFjdG9ycygtbiwgbGltaXQpCiAgICBmcyA9IFtdCiAgICB3aGlsZSBuICUgMiA9PSAwOgogICAgICAgIG4gPSBuIC8vIDI7IGZzID0gZnMgKyBbMl0KICAgIGlmIG4gPT0gMTogcmV0dXJuIGZzCiAgICB3aGlsZSBub3QgaXNfcHJpbWUobik6CiAgICAgICAgZiA9IHJob19mYWN0b3IobiwgMSwgbGltaXQpCiAgICAgICAgbiA9IG4gLyBmCiAgICAgICAgZnMgPSBmcyArIFtmXQogICAgcmV0dXJuIHNvcnRlZChmcyArIFtuXSkKCnByaW50IHJob19mYWN0b3JzKDk5OTk5NjQwMDAwMjQzKQ==