UPDATE: This is a documentation bug. See the last comment.
random(f)
, for f
a denormalized number, underrepresents 0.0 and can return f
, even though the spec says that the result is strictly less than the argument.
Simple demo:
multiset(list) :=
block([ar,idxlist,vallist,retval:[]],local(ar),
list: map('ratdisrep,list),
ar[x]:=0,
for i in list do ar[i]:ar[i]+1,
idxlist: rest(arrayinfo(ar),2),
vallist: listarray(ar),
kill(ar),
while idxlist # [] do
( push([first(first(idxlist)),first(vallist)],retval),
idxlist: rest(idxlist),
vallist: rest(vallist) ),
sort(retval))$
fpprintprec:6$
multiset(makelist(random(2e-323),i,1,1000)) =>
[[0.0, 113], [4.94066e-324, 262], [9.88131e-324, 247], [1.4822e-323, 259], [1.97626e-323, 119]]
2e-323 => 1.97626e-323
0.0 is about half as common as it should be, and 1.97626e-323 shouldn't be returned at all.
I have tested this with larger denormalized floats as well, e.g., random(4.94e-320)
.
This behavior is inherited from SBCL's random
function -- I will file a bug report with them.
Maxima 5.45.1 SBCL 2.0.0 Windows
Maxima's rand-mt19937.lisp says it's from cmucl (not that that matters much; cmucl doesn't use mt19937 anymore.)
Perhaps these are relevant: https://vigna.di.unimi.it/xorshift/, Look for the section titled "Generating uniform doubles in the unit interval". There's a link there to https://vigna.di.unimi.it/xorshift/random_real.c for generating these.
This seems like such a niche area.
Yes, it is a niche area, but I came across it honestly -- I wasn't looking for problems with random, but rather trying to use random for some testing. The current code for random(float) must be something like
random(2^53)/2.0^53*float
, which isn't quite right for denormalized numbers.I believe this is a correct fix (even if not elegant): if random(a)=a, return 0.0.
i think I'm at least partially wrong. More tomorrow.
The comment in the code for
%random-double-float
saysI believe this means the lowest bit of the fraction is always 0, so we're missing one bit.
This is then multiplied by the arg of
random
.The xoroshiro link says that to get a random float in the interval the same interval you should probably be
(x >> 1) * 2^-53
, werex
is a 64-bit unsigned integer.I'm guessing, but if this matters, we should use the algorithm given in
random_real
.I was wrong. The current distribution of values of
random(<denorm>)
is in fact correct. It is the documentation that is subtly wrong.random(x)
for float x in effect generates a random real (not float) number in the interval[0,x)
. It then rounds it to the nearest floating-point number.The interval
[0,x)
can be broken up into equivalence classes defined by "nearest floating point number". Calling the smallest float ε, those intervals are(0,ε/2)(ε/2,3*ε/2)(3*ε/2,5ε/2) ... (x-ε/2,x)
. That is, the first and last intervals are half the size of the other intervals, so should appear half as often as the others.The only error is that
random
is defined to return a number "less than <x>". In fact, it should be "less than or equal to <x>". For <x> in a normal range, this hardly matters, since the probability of returning <x> is tiny.</x></x></x></x>So I will change this to a documentation bug.
But based on the code,
random(x)
should return a value less thanx
. Unless I'm missing something. And we should be consistent. Ifrandom(x)
can returnx
for floating-point values ofx
, then it should also returnx
whenx
is an integer. For integers, I don't think that's possible because it does(rem bits x)
.Note that this method has a bias. Perhaps it should be replaced by a method that doesn't have the bias. But that's a different bug. Perhaps people aren't generating enough values to notice the bias.
The only time anyone will notice that
x
is a possible value forrandom(x)
(with float x) is for small denormalized numbers. For normalized numbers, the probability of generating exactlyx
is something like 2^-52 = 10^-16.This is very different from the integer case.
Ok, I think
random(x)
producingx
only happens whenx
is a denormal because there aren't enough bits left to represent the product of0.9999999999999998*x
. The number0.999...8
comes from how%random-double-float
generates a large integer and smashes that into a double. Whenx
is not a denormal, there are enough bits so that0.999...8*x
has enough bits so that the product is notx
. I think.