Implemented Smith's method, which avoid overflow for all but
the most extreme exponents. This agrees with what is used in NumPy.