BTW, if we have to "do something", how about changing

return 1.0 / u ** (1.0/alpha)

to the mathematically equivalent

return (1.0 / u) ** (1.0/alpha)

? Not sure about Linux-y boxes, but on Windows that would raise OverflowError instead of ZeroDivisionError. Which makes more sense, because the Pareto variable we're _trying_ to produce is too large to represent.  If it returns math.inf on some box instead, same thing.

Or, also faster, and suffering fewer rounding errors,

return  u ** (-1.0 / alpha)