In some simulations I am encountering a problem due to machine precision. I have to conduct an exponentiation operation on a number that lies in (0, 1). That number comes from a random number generator, while the exponent is determined by an optimization algorithm.

My question is how I can determine the maximum exponent n such that rand()^n does not result in 0.0?

Or equivalently, how I can determine the maximum exponent n such that rand()^(-n) does not result in Inf?

These arenâ€™t actually equivalent questions due to the way subnormals work. The key things to realize are that rand always generates a number in (floatmin(),1), and therefore the answer is the power which will cause floatmin() to under/overflow. This number is log(nextfloat(0.0))/log(floatmax()) or log(floatmax())/log(floatmin()) respectively. (0.9534450651769087 or -1.0019569471624266)

Thank you. It makes sense that the numbers are so low.

I was thinking of using clamp to limit the range of n, but ~0.95 is very low for my application. I am expecting n be up to 100 to 150. I understand that the operation depends on a random number so I can never know for sure what number will be drawn, but do you have any suggestions regarding that?

Specifically, what I am trying to do is generate a correlated exponential distribution using the Clayton copula. The input is simply a matrix of random numbers in (0, 1):

randmat = rand(100, 2)
matkk = similar(randmat)
function claytonsample!(matkk, Ď„; randmat=randmat)
matkk .= randmat
Ď„ == 0 && return matkk
n = size(matkk, 1)
for i in 1:n
v = matkk[i, 2]
u = matkk[i, 1]
matkk[i, 2] = (1 - u^(-Ď„) + u^(-Ď„)*v^(-(Ď„/(1 + Ď„))))^(-1/Ď„)
end
return matkk
end
invsurvival(x) = -log(x)
function simulkk!(matkk, Ď„; randmat=randmat)
claytonsample!(matkk, Ď„; randmat=randmat)
matkk .= invsurvival.(matkk)
end
simulkk!(matkk, 300; randmat=randmat)

The problem is in the terms u^(-Ď„) above. I was thinking of clamping Ď„. Iâ€™d rather not mess with the â€śseedâ€ť random numbers in randmat.

Wait, is rand() in [0, 1) or [floatmin(), 1) (so basically, (0,1)) or (floatmin(), 1)? I canâ€™t find it in the docs, but numpyâ€™s np.random.rand() claims to be [0,1) and I imagine Julia would follow the same convention.

Numerically, no matter what youâ€™re doing, division by anything even close to floatmin will be unstable. It will however be â€śsafeâ€ť, in that it will not error out. Even division by 0.0 is not an error, it just returns Inf.

In my case, it arises in a unit test that uses randomly generated input data, so while we would expect higher numbers in â€śrealâ€ť input, all the test requires is that 1/rand() be finite.

Donâ€™t you want your unit tests to use realistic numbers? Or just add 1.0 to the random numbers to be safe (or maybe put in some handling of Inf in your program)

A secondary motivation for my question is so that I can write a sentence like "in this experiment, x was drawn randomly from [0,1)" in a scientific paper describing a computational experiment (one that doesnâ€™t involve computing 1/x). I would like my sentence to be accurate.