Good idea, but yet another that really belongs in an `imath` module (which doesn't yet exist).

How ambitious should it be?  Sympy supplies a `factorint()` function for this, which uses 4 approaches under the covers:  perfect power, trial division, Pollard rho, and Pollard p-1.  All relatively simple to code with trivial memory burden, but not really practical (too slow) for "hard" composites well within the practical range of advanced methods.

But I'd be happy enough to settle for that.