`random.uniform(a, b)` can overflow to values outside of `[a, b]`

[Liam-DeVoe]

Bug report

Bug description:

random.uniform can overflow to inf outside the range of [a, b] in some cases:

import random

random.uniform(-1e308, 1e308)  # math.inf

The implementation of random.uniform is:

def uniform(self, a, b):
    return a + (b - a) * self.random()

where (b - a) can overflow to math.inf: assert 1e308 * 2 == math.inf.

Since the docstring of random.uniform guarantees that the return value will be in the range [a, b], I would classify this as a bug.

I see a few solutions:

• Error on overflow.
  • numpy does this: numpy.random.uniform(-1e308, 1e308) raises OverflowError.
• Document that random.uniform may overflow to math.inf, even if b < math.inf.
• Clamp out-of-bounds values (only if b - a == math.inf, for performance).

CPython versions tested on:

3.12

Operating systems tested on:

macOS

[ZeroIntensity]

cc @skirpichev

[rhettinger]

Sorry, this isn't a corner case that I care about. In the long history of the random module, this has not been a problem in actual practice. I don't think it even warrants a note in the docs.

In a way, it is no more interesting than hypothesis or another fuzzer revealing that operator.sub(1e308, -1e308) returns inf. There isn't any real value in documenting that every float subtraction in Python can be overflowed. That wouldn't help users and would just be a distraction. This applies broadly. Hypothesis will find overflows and numerical issues in a lot of floating point code (sin/asin roundtrip, basic trig identities, or even simple algebraic relations like a + b - b == 0 or a + b + c == c + b + a.

[Liam-DeVoe]

I understand this perspective. In my opinion though, from a developer's perspective, it's nice to be able to work with rock-solid guarantees from the standard library. So if Python says that random.uniform(a, b) is in the range [a, b], I can trust that really is true for all 64-bit floating point numbers, even in unlikely cases.

I agree "floating point operations can overflow" is an uninteresting disclaimer. But "this function in the standard library can return values outside of its promised range" is more interesting.

Because python doesn't document that sin and asin roundtrip for all 64-bit floating point numbers, I wouldn't consider that a bug if I found a counterexample. But python does document that all values of random.uniform(a, b) are in the range [a, b]. That's not true (albeit only in extreme cases!), so I would consider that a bug, either in the implementation, or in the docs.

Here is the documentation for random.uniform, for reference:

Return a random floating-point number N such that a <= N <= b for a <= b and b <= N <= a for b < a.

The end-point value b may or may not be included in the range depending on floating-point rounding in the expression `a + (b-a) * random()`.

Noting that it already contains a special-case for the inclusion of b based on floating-point details.

[skirpichev]

There isn't any real value in documenting that every float subtraction in Python can be overflowed.

Though, from docs - users can only guess that subtraction happens. OverflowError might be better here, like in most math's module functions in similar cases.

Note that other variants of linear interpolation doesn't overflow here:

def uniform2(x, a, b):
    return a + b*x - a*x
def uniform3(x, a, b):
    return fma(x, b, fma(-x, a, a))

But it's easy to find examples, when uniform2 is not monotonic. I suspect same is true for the uniform3. Edit: indeed, for a, b = 1.2, 1.7 and x1, x2 = 0.49581224138185065, 0.4958122413818507.

this has not been a problem in actual practice

@tybug, did your example arise in practice? Assuming float's are IEEE doubles, uniform(-1e308, 1e308) for me looks too far from approximation of uniformly distributed real numbers in [-1e308, 1e308], i.e. near the ends - all available floats are integers.

On this ground, I would rather agree that this corner case not worth efforts.

[rhettinger]

If this ever arises in a real-world context (not out of a fuzzer), I will add a note to the docs.

In the meantime, it is my opinion that a doc edit would just be distracting clutter.

We definitely don't want to slow down and clutter the code itself to test for inf and error out. That would make everyone pay for a problem that no one has. I know that this might feel "impure" or less-than-perfect, but it is practical. I really don't think there is a real world issue here that we care about. It hasn't arisen in the more than two and a half decade history of the module.

Thank you for the suggestion, but I will pass on it for now.

[Liam-DeVoe]

I ran into this while testing my code with Hypothesis, which I suspect you'll consider "not real-world use", though I would say that testing code is very real world :)

To clarify, I was not intentionally searching for bugs in the standard library. I simply had some code that called random.uniform somewhere, and was testing a higher-level invariant that my code did not produce out-of-bounds values. It did, not because my code was incorrect, but because I was relying on the promised bounds of random.uniform. My first thought was certainly not "the standard library guarantees are wrong".

Thanks all!

[Zac-HD]

I encouraged @tybug to open an issue because we ran into this as a practical problem, due to random.uniform() explicitly violating the documented bounds. The documentation explictly mentions that b may or may not be included due to floating-point rounding, which made it particularly surprising that we could generate a < b < N.

Of course it's reasonable to decide that clamping would not be a net improvement, but in that case I think the documentation should be corrected.

[skirpichev]

This "issue" happens for triangular(), again with meaningless parameters.

Edit:
Something like this should work without overflows:

def uniform(self, a, b):
    full_range = b - a
    if full_range != _inf:
        return a + full_range * self.random()
    a2, b2 = a/2, b/2
    mid = b2 + a2
    half_range = b2 - a2
    return mid + half_range * (2*self.random() - 1)

Though, it costs x1.2-1.3 speed on small ranges.

[mdickinson]

@skirpichev Simpler to use the usual expression applied to the halved inputs in the exceptional case, and then double the result: return (a2 + (b2 - a2) * self.random()) * 2.0. That way you keep the guarantee that (assuming IEEE 754, round-ties-to-even) the result is again in [a, b] - the halving operations and the doubling operation are exact, so the same analysis as for the main case carries through. (Corner case: halving isn't necessarily exact if either a or b is subnormal, but in that case it's impossible for b - a to overflow in the first place, so we don't end up in the exceptional branch of the code.)

The analysis of the midpoint-based code (mid + half_range * ...) is more complicated, and I don't think it's true any more that the output of that code is guaranteed to be in the closed interval [a, b]: in the case where self.random() returns 0.0, you end up with (b2 + a2) - (b2 - a2), which could be smaller than a.

[skirpichev]

(b2 + a2) - (b2 - a2), which could be smaller than a

Yep, e.g. a, b = 0.2550690257394217, 0.763774618976614. Thanks for the fix.

Though, I worry that overflow testing is too costly in non-exceptional case:

Benchmark	ref	patch
uniform(2.7, 3.14)	596 ns	697 ns: 1.17x slower

Using uniform() (or triangular()) in cases, when range close to overflow (or might overflow) doesn't make much sense for me. Did I miss something here?