Uniform Generator hangs for certain limits.

[GUIpsp]

See here for a MWE:

https://github.com/GUIpsp/rand_uniform_hang/blob/main/src/main.rs

I also tested that #1287 does not fix this issue.

[WarrenWeckesser]

The "hang" is actually this loop

rand/src/distributions/uniform.rs

Lines 861 to 867 in 0f3eced

	loop {
	let mask = (scale * max_rand + low).ge_mask(high);
	if !mask.any() {
	break;
	}
	scale = scale.decrease_masked(mask);
	}

taking a very long time to complete. That loop will require roughly spacing(low)/(2*spacing(scale)) iterations to complete, where spacing(x) is the unit of least precision of x. For the values in the MWE, this is 17179869184 iterations! Here's the calculation in a Python session:

In [52]: import numpy as np

In [53]: low = 0.99999999997819644

In [54]: high = 1.0

In [55]: scale = high - low

In [56]: np.spacing(low)/(2*np.spacing(scale))
Out[56]: 17179869184.0

The problem is that in the loop, scale is repeatedly decreased by the smallest possible amount until scale*max_rand + low is less than high. For inputs where high - low is much smaller than low, this will require a large number of iterations.

[dhardy]

This is about float sampling so #1287 is unrelated; #1289 is more relevant (but I don't think it will solve this particular issue).

Decreasing the scale was added in #477 (@sicking).

[WarrenWeckesser]

FYI: NumPy has an issue about their uniform generator claiming to generate variates in [low, high) even though it can, in fact, generate high. The result of that discussion was to simply document that, because of floating point rounding, high might be generated in the output even though the underlying standard generator for [0, 1) guarantees that 1 is not generated.

[vks]

Floating point rounding is not a great excuse, it's easy to fix returning high by rejection sampling. You might still decide not to do that, i.e. for performance reasons, and users who cannot accept high being returned can do the rejection sampling themselves.

(FWIW, I ran into similar issues with the C++ standard libraries.)

[GUIpsp]

I agree with whatr @vks, but in addition, I'd rather it return high than hit this degenerate case, having to wait for two hours.

[dhardy]

but in addition, I'd rather it return high than hit this degenerate case, having to wait for two hours.

To be clear, this reported bug is not something we should just document and ignore.

The least intrusive fix would probably be to calculate a more precise decrement.

---

Uniform implements sampling for both [a, b) and [a, b], each for single-sample and distribution sampling. Perhaps we should simply by only supporting [a, b] (i.e. inclusive ranges)?

Complication: impl<X: SampleUniform> TryFrom<Range<X>> for Uniform<X> uses the half-open variant. We can simply subtract 1 for integer types — but there are no numeric traits in std and only rand_distr depends on num-traits.

[WarrenWeckesser]

The least intrusive fix would probably be to calculate a more precise decrement.

I think I have a fix. Once I work out one particular edge case, I'll submit a pull request.