# 100011322432435545
def isqrt(n): # newton
x = n; y = (x + 1) // 2
while y < x:
x = y; y = (x + n // x) // 2
return x
def isSquare(n):
# def q(n):
# from sets import Set
# s, sum = Set(), 0
# for x in xrange(0,n):
# t = pow(x,2,n)
# if t not in s:
# s.add(t)
# sum += pow(2,t)
# return sum
# q(32) => 33751571
# q(27) => 38348435
# q(25) => 19483219
# q(19) => 199411
# q(17) => 107287
# q(13) => 5659
# q(11) => 571
# q(7) => 23
# 99.82% of non-squares caught by bloom
# filters before square root calculation
if 33751571>>(n%32)&1==0: return False
if 38348435>>(n%27)&1==0: return False
if 19483219>>(n%25)&1==0: return False
if 199411>>(n%19)&1==0: return False
if 107287>>(n%17)&1==0: return False
if 5659>>(n%13)&1==0: return False
if 571>>(n%11)&1==0: return False
if 23>>(n% 7)&1==0: return False
s = isqrt(n)
if s * s == n: return s
return False
def isPrime(n, k=5): # miller-rabin
from random import randint
ps = [2,3,5,7,11,13,17,19,23,29]
if n < 2: return False
for p in ps:
if n%p == 0: return n==p
s, d = 0, n-1
while d%2 == 0:
s, d = s+1, d/2
for i in xrange(k):
x = pow(randint(2, n-1), d, n)
if x == 1 or x == n-1: continue
for r in xrange(1, s):
x = (x * x) % n
if x == 1: return False
if x == n-1: break
else: return False
return True
def factors(n):
def gcd(a, b): # euclid
if b == 0: return a
return gcd(b, a%b)
# 2,3,5-wheel to cube root of n
wheel = [1,2,2,4,2,4,2,4,6,2,6]
w, f, fs = 0, 2, []
while f*f*f <= n:
# print f, n, fs
while n % f == 0:
fs.append(f); n /= f
if n == 1: return fs
if isPrime(n):
fs.append(n)
return fs
f, w = f+wheel[w], w+1
if w == 11: w = 3
# n must be semi-prime, so use
# william hart's fermat variant
if isSquare(n):
# print "square", n
f = isqrt(n);
fs.append(f); fs.append(f)
return fs
for i in xrange(1,n):
s = isqrt(n*i)
if s*s <= n*i: s += 1
m = pow(s,2,n)
# print i, s, m
if isSquare(m):
t = isqrt(m); f = gcd(n, s-t)
fs.append(f); fs.append(n/f)
return sorted(fs)
print factors(100011322432435545)
IyAxMDAwMTEzMjI0MzI0MzU1NDUKCmRlZiBpc3FydChuKTogIyBuZXd0b24KICAgIHggPSBuOyB5ID0gKHggKyAxKSAvLyAyCiAgICB3aGlsZSB5IDwgeDoKICAgICAgICB4ID0geTsgeSA9ICh4ICsgbiAvLyB4KSAvLyAyCiAgICByZXR1cm4geAoKZGVmIGlzU3F1YXJlKG4pOgogICAgIyBkZWYgcShuKToKICAgICMgICAgIGZyb20gc2V0cyBpbXBvcnQgU2V0CiAgICAjICAgICBzLCBzdW0gPSBTZXQoKSwgMAogICAgIyAgICAgZm9yIHggaW4geHJhbmdlKDAsbik6CiAgICAjICAgICAgICAgdCA9IHBvdyh4LDIsbikKICAgICMgICAgICAgICBpZiB0IG5vdCBpbiBzOgogICAgIyAgICAgICAgICAgICBzLmFkZCh0KQogICAgIyAgICAgICAgICAgICBzdW0gKz0gcG93KDIsdCkKICAgICMgICAgIHJldHVybiBzdW0KICAgICMgcSgzMikgPT4gMzM3NTE1NzEKICAgICMgcSgyNykgPT4gMzgzNDg0MzUKICAgICMgcSgyNSkgPT4gMTk0ODMyMTkKICAgICMgcSgxOSkgPT4gICAxOTk0MTEKICAgICMgcSgxNykgPT4gICAxMDcyODcKICAgICMgcSgxMykgPT4gICAgIDU2NTkKICAgICMgcSgxMSkgPT4gICAgICA1NzEKICAgICMgcSg3KSAgPT4gICAgICAgMjMKICAgICMgOTkuODIlIG9mIG5vbi1zcXVhcmVzIGNhdWdodCBieSBibG9vbQogICAgIyBmaWx0ZXJzIGJlZm9yZSBzcXVhcmUgcm9vdCBjYWxjdWxhdGlvbgogICAgaWYgMzM3NTE1NzE+PihuJTMyKSYxPT0wOiByZXR1cm4gRmFsc2UKICAgIGlmIDM4MzQ4NDM1Pj4obiUyNykmMT09MDogcmV0dXJuIEZhbHNlCiAgICBpZiAxOTQ4MzIxOT4+KG4lMjUpJjE9PTA6IHJldHVybiBGYWxzZQogICAgaWYgICAxOTk0MTE+PihuJTE5KSYxPT0wOiByZXR1cm4gRmFsc2UKICAgIGlmICAgMTA3Mjg3Pj4obiUxNykmMT09MDogcmV0dXJuIEZhbHNlCiAgICBpZiAgICAgNTY1OT4+KG4lMTMpJjE9PTA6IHJldHVybiBGYWxzZQogICAgaWYgICAgICA1NzE+PihuJTExKSYxPT0wOiByZXR1cm4gRmFsc2UgCiAgICBpZiAgICAgICAyMz4+KG4lIDcpJjE9PTA6IHJldHVybiBGYWxzZQogICAgcyA9IGlzcXJ0KG4pCiAgICBpZiBzICogcyA9PSBuOiByZXR1cm4gcwogICAgcmV0dXJuIEZhbHNlCgpkZWYgaXNQcmltZShuLCBrPTUpOiAjIG1pbGxlci1yYWJpbgogICAgZnJvbSByYW5kb20gaW1wb3J0IHJhbmRpbnQKICAgIHBzID0gWzIsMyw1LDcsMTEsMTMsMTcsMTksMjMsMjldCiAgICBpZiBuIDwgMjogcmV0dXJuIEZhbHNlCiAgICBmb3IgcCBpbiBwczoKICAgICAgICBpZiBuJXAgPT0gMDogcmV0dXJuIG49PXAKICAgIHMsIGQgPSAwLCBuLTEKICAgIHdoaWxlIGQlMiA9PSAwOgogICAgICAgIHMsIGQgPSBzKzEsIGQvMgogICAgZm9yIGkgaW4geHJhbmdlKGspOgogICAgICAgIHggPSBwb3cocmFuZGludCgyLCBuLTEpLCBkLCBuKQogICAgICAgIGlmIHggPT0gMSBvciB4ID09IG4tMTogY29udGludWUKICAgICAgICBmb3IgciBpbiB4cmFuZ2UoMSwgcyk6CiAgICAgICAgICAgIHggPSAoeCAqIHgpICUgbgogICAgICAgICAgICBpZiB4ID09IDE6IHJldHVybiBGYWxzZQogICAgICAgICAgICBpZiB4ID09IG4tMTogYnJlYWsKICAgICAgICBlbHNlOiByZXR1cm4gRmFsc2UKICAgIHJldHVybiBUcnVlCgpkZWYgZmFjdG9ycyhuKToKICAgIGRlZiBnY2QoYSwgYik6ICMgZXVjbGlkCiAgICAgICAgaWYgYiA9PSAwOiByZXR1cm4gYQogICAgICAgIHJldHVybiBnY2QoYiwgYSViKQogICAgIyAyLDMsNS13aGVlbCB0byBjdWJlIHJvb3Qgb2YgbgogICAgd2hlZWwgPSBbMSwyLDIsNCwyLDQsMiw0LDYsMiw2XQogICAgdywgZiwgZnMgPSAwLCAyLCBbXQogICAgd2hpbGUgZipmKmYgPD0gbjoKICAgICAgICAjIHByaW50IGYsIG4sIGZzCiAgICAgICAgd2hpbGUgbiAlIGYgPT0gMDoKICAgICAgICAgICAgZnMuYXBwZW5kKGYpOyBuIC89IGYKICAgICAgICBpZiBuID09IDE6IHJldHVybiBmcwogICAgICAgIGlmIGlzUHJpbWUobik6CiAgICAgICAgICAgIGZzLmFwcGVuZChuKQogICAgICAgICAgICByZXR1cm4gZnMKICAgICAgICBmLCB3ID0gZit3aGVlbFt3XSwgdysxCiAgICAgICAgaWYgdyA9PSAxMTogdyA9IDMKICAgICMgbiBtdXN0IGJlIHNlbWktcHJpbWUsIHNvIHVzZQogICAgIyB3aWxsaWFtIGhhcnQncyBmZXJtYXQgdmFyaWFudAogICAgaWYgaXNTcXVhcmUobik6CiAgICAgICAgIyBwcmludCAic3F1YXJlIiwgbgogICAgICAgIGYgPSBpc3FydChuKTsKICAgICAgICBmcy5hcHBlbmQoZik7IGZzLmFwcGVuZChmKQogICAgICAgIHJldHVybiBmcwogICAgZm9yIGkgaW4geHJhbmdlKDEsbik6CiAgICAgICAgcyA9IGlzcXJ0KG4qaSkKICAgICAgICBpZiBzKnMgPD0gbippOiBzICs9IDEKICAgICAgICBtID0gcG93KHMsMixuKQogICAgICAgICMgcHJpbnQgaSwgcywgbQogICAgICAgIGlmIGlzU3F1YXJlKG0pOgogICAgICAgICAgICB0ID0gaXNxcnQobSk7IGYgPSBnY2Qobiwgcy10KQogICAgICAgICAgICBmcy5hcHBlbmQoZik7IGZzLmFwcGVuZChuL2YpCiAgICAgICAgICAgIHJldHVybiBzb3J0ZWQoZnMpCgpwcmludCBmYWN0b3JzKDEwMDAxMTMyMjQzMjQzNTU0NSk=