The following code raises a ZeroDivisionError:

from random import Random

Random(14064741636871487939).paretovariate(0.01)

This raises:

random.py, line 692, in paretovariate
return 1.0 / u ** (1.0/alpha)
ZeroDivisionError: float division by zero

That specific stack trace is from 3.8.5 but I've tried this on 3.6 and 3.7 as well.

It's a little hard to tell what the intended correct range of parameter values is for paretovariate, and this may just be lack of validation for the alpha < 1 case - perhaps that was never intended to be supported?

Based on some very informal inspection, what seems to happen is that the probability of getting a ZeroDivisionError approaches one as you approach zero. They rarely occur at this level of alpha (0.01), but for alpha=0.001 they seem to occur just under half the time, and for alpha=0.0001 they occur more than 90% of the time.

(For the interested, this bug was found by Hypothesis as part of its own testing of our integration with the Random API)

I guess on actual inspection of the code (which I should have done before, sorry) it's obvious why this happens: u ** (1.0 / alpha) will round to 0.0 for small values of u, and the smaller alpha is the higher the probability of that happening, in that it happens with probability c ** alpha for some c.

I'm inclined to ignore this. No actual user has complained about this, and I doubt any ever will: it's got to be rare as hen's teeth to use a parameter outside of, say, [0.1, 10.0], in real life. The division error can't happen for those. For "small" alpha, the distribution simply does contain values too large to represent as binary doubles. So the only other defensible thing to do be done in this case is to return math.inf (or raise OverflowError) - but why slow it down for something no actual user cares about?

The point of the "u = 1.0 - self.random()" line was to prevent the case where *u* was exactly equal to zero; however, as you noted, raising a very small number to a small can round down to zero.

We could wrap the calculation in a try/except to catch the ZeroDivisionError but that is probably insufficient. The codomain can easily spew beyond the range of a C double:

   >>> alpha = 0.01
   >>> u = 0.005
   >>> 1.0 / u ** (1.0 / alpha)
   1.2676506002282268e+230             <== Huge!

We could clamp the values as is done in gammavariate(). Or we can raise a clean Overflow or Underflow exception but that would be unpleasant for users if it were to arise unexpectedly.

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)

I don't feel a need to change anything, my post hit a few seconds after yours :-)

That said, "return u ** (-1.0 / alpha)" would be a slight, but nice improvement.

It might also be reasonable to document a safe range of alpha values.

I should say, I don't actually care about this bug at all: I only ran into it because of trying to recreate the random API and testing that my recreation worked sensibly. I thought I'd do the good citizen thing and report it, but I'm totally fine with any resolution and given that that code has been unchanged for at least 18 years I doubt anyone else has ever run into this problem.

I do think it would be beneficial to document the range of reasonable alpha on paretovariate - the lack of such documentation is the only reason I hit this at all - but I think leaving the behaviour as is would be fine.

+1 for "u ** (-1.0 / alpha)"!

New changeset 5c32709 by Raymond Hettinger in branch 'master':
bpo-41421: Algebraic simplification for random.paretovariate() (GH-21695)
5c32709

How about this wording?

If alpha is smaller than 0.1, an occasional OverflowError can
arise when the variate exceeds the range of Python float.

I like using 0.1 because it's easy and gives us some wiggle room. The actual cutoff is lower:

>>> alpha = 0.05079
>>> sys.float_info.epsilon ** (-1.0 / alpha)

    >>> alpha = 0.05078
    >>> sys.float_info.epsilon ** (-1.0 / alpha)
    Traceback (most recent call last):
       File "<pyshell#48>", line 1, in <module>
         sys.float_info.epsilon ** (-1.0 / alpha)
    1.590779741838475e+308

That text is fine, if you feel something needs to be said at all. I really don't. A Pareto distribution of this kind with parameter <= 1.0 has infinite expected value - VERY long tail. Anyone who knows what they're doing already knows that. The reason the implementation can't "blow up" for parameters >= (roughly) 0.1 isn't that IEEE doubles have such a large dynamic range but rather that they can't represent a number < 1.0 larger than 1 - 2**-53 (so u = 1 - random.random() is always at least 2**-53). The actual distribution has infinite expected value nonetheless, but the implementation is incapable of generating any of its very large values (which, while very rare, are also very large).

Okay, marking this as closed.