# nearly square divisors, revisited
def factors(n): # 2,3,5-wheel
ws = [1,2,2,4,2,4,2,4,6,2,6]
w, f, fs = 0, 2, []
while f*f <= n:
while n % f == 0:
fs.append(f)
n /= f
f, w = f + ws[w], w+1
if w == 11: w = 3
if n < f*f: fs.append(n)
return fs
def knap(target, xs):
if xs == []: return 0
if target < xs[0]:
return knap(target, xs[1:])
else:
return max(xs[0] + knap(target-xs[0],xs[1:]), \
knap(target, xs[1:]))
def nsd(n):
from math import log, exp
fs = map(log, factors(n)[::-1])
b = int(round(exp(knap(sum(fs)/2, fs))))
return n/b, b
print nsd(224403121196654400)
IyBuZWFybHkgc3F1YXJlIGRpdmlzb3JzLCByZXZpc2l0ZWQKCmRlZiBmYWN0b3JzKG4pOiAjIDIsMyw1LXdoZWVsCiAgICB3cyA9IFsxLDIsMiw0LDIsNCwyLDQsNiwyLDZdCiAgICB3LCBmLCBmcyA9IDAsIDIsIFtdCiAgICB3aGlsZSBmKmYgPD0gbjoKICAgICAgICB3aGlsZSBuICUgZiA9PSAwOgogICAgICAgICAgICBmcy5hcHBlbmQoZikKICAgICAgICAgICAgbiAvPSBmCiAgICAgICAgZiwgdyA9IGYgKyB3c1t3XSwgdysxCiAgICAgICAgaWYgdyA9PSAxMTogdyA9IDMKICAgIGlmIG4gPCBmKmY6IGZzLmFwcGVuZChuKQogICAgcmV0dXJuIGZzCgpkZWYga25hcCh0YXJnZXQsIHhzKToKICAgIGlmIHhzID09IFtdOiByZXR1cm4gMAogICAgaWYgdGFyZ2V0IDwgeHNbMF06CiAgICAgICAgcmV0dXJuIGtuYXAodGFyZ2V0LCB4c1sxOl0pCiAgICBlbHNlOgogICAgICAgIHJldHVybiBtYXgoeHNbMF0gKyBrbmFwKHRhcmdldC14c1swXSx4c1sxOl0pLCBcCiAgICAgICAgICAgICAgICAgICBrbmFwKHRhcmdldCwgeHNbMTpdKSkKCmRlZiBuc2Qobik6CiAgICBmcm9tIG1hdGggaW1wb3J0IGxvZywgZXhwCiAgICBmcyA9IG1hcChsb2csIGZhY3RvcnMobilbOjotMV0pCiAgICBiID0gaW50KHJvdW5kKGV4cChrbmFwKHN1bShmcykvMiwgZnMpKSkpCiAgICByZXR1cm4gbi9iLCBiCgpwcmludCBuc2QoMjI0NDAzMTIxMTk2NjU0NDAwKQ==