# Exponentiation and machine precision
**URL:** https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024
**Category:** Numerics
**Tags:** precision
**Created:** [April 5, 2022, 12:48am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024 "2022-04-05T00:48:00Z")
**Posts on this page:** 19
**Page:** 1
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 5, 2022, 12:48am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/1 "2022-04-05T00:48:00Z")
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`?
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [April 5, 2022, 1:02am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/2 "2022-04-05T01:02:00Z")
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`)
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 5, 2022, 1:12am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/3 "2022-04-05T01:12:53Z")
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?
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [April 5, 2022, 1:14am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/4 "2022-04-05T01:14:15Z")
What you actually want to do is generate a random number with the desired distribution.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 5, 2022, 1:20am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/5 "2022-04-05T01:20:36Z")
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):
```julia
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`.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [April 5, 2022, 1:32am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/6 "2022-04-05T01:32:08Z")
Herbie gives some useful suggestions here. [https://herbie.uwplse.org/demo/e7768dd3187b2711e9ab6825ebe15f0057554539.ea9754c765d6f2e3443236047de835c7bf56f0de/graph.html](https://herbie.uwplse.org/demo/e7768dd3187b2711e9ab6825ebe15f0057554539.ea9754c765d6f2e3443236047de835c7bf56f0de/graph.html) (although the last branch is incorrect, it should just be return `t_0`.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 5, 2022, 1:40am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/7 "2022-04-05T01:40:26Z")
This is unrelated, but I get a very different answer from yours:
[https://herbie.uwplse.org/demo/1b508faa0f15ca8e4e244eadb39ba19c403cdcf4.ea9754c765d6f2e3443236047de835c7bf56f0de/graph.html](https://herbie.uwplse.org/demo/1b508faa0f15ca8e4e244eadb39ba19c403cdcf4.ea9754c765d6f2e3443236047de835c7bf56f0de/graph.html)
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [April 5, 2022, 1:41am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/8 "2022-04-05T01:41:18Z")
I might have typed something in wrong, never mind, herbie just got drunk when evaluating yours.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 5, 2022, 1:43am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/9 "2022-04-05T01:43:13Z")
Oh man, that is some serious drinking problem. Hopefully both are equivalent.
---
### Author:  [@maxkapur](https://discourse.julialang.org/u/maxkapur)
#### Post date: [April 5, 2022, 3:23am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/10 "2022-04-05T03:23:24Z")
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.
[https://numpy.org/devdocs/reference/random/generated/numpy.random.rand.html](https://numpy.org/devdocs/reference/random/generated/numpy.random.rand.html)
Posting in the same topic because I have a fundamentally similar question–I want to know whether I can safely do `1 / rand(Float16)`.
---
### Author:  [@gustaphe](https://discourse.julialang.org/u/gustaphe)
#### Post date: [April 5, 2022, 4:07am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/11 "2022-04-05T04:07:22Z")
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`.
---
### Author:  [@maxkapur](https://discourse.julialang.org/u/maxkapur)
#### Post date: [April 5, 2022, 4:18am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/12 "2022-04-05T04:18:39Z")
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.
---
### Author:  [@gustaphe](https://discourse.julialang.org/u/gustaphe)
#### Post date: [April 5, 2022, 4:34am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/13 "2022-04-05T04:34:34Z")
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)
---
### Author:  [@maxkapur](https://discourse.julialang.org/u/maxkapur)
#### Post date: [April 5, 2022, 4:43am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/14 "2022-04-05T04:43:22Z")
Sure, depends on the context.
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.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [April 5, 2022, 4:51am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/15 "2022-04-05T04:51:27Z")
`rand()` is uniform over `[0,1)` but the posible return values are not all floating point numbers in `[0,1)`.
---
### Author:  [@maxkapur](https://discourse.julialang.org/u/maxkapur)
#### Post date: [April 5, 2022, 5:05am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/16 "2022-04-05T05:05:50Z")
Which numbers are left out?
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [April 5, 2022, 5:09am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/17 "2022-04-05T05:09:46Z")
`rand()` generates numbers that are `i*2^-53` where `i` is a uniform integer in `[0,2^53)`
---
### Author:  [@gustaphe](https://discourse.julialang.org/u/gustaphe)
#### Post date: [April 5, 2022, 5:42am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/18 "2022-04-05T05:42:15Z")
Eh, it’s `[0, 1)` with `p<0.05`. You’re more likely to have a cosmic bit flip place you outside the interval than for it to literally pick `0.0`
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [June 8, 2022, 9:16am UTC](https://discourse.julialang.org/t/exponentiation-and-machine-precision/79024/19 "2022-06-08T09:16:10Z")
Why not use the [`Copulas.jl`](https://github.com/lrnv/Copulas.jl) package ?
```julia
using Copulas, Distributions, Random
C = ClaytonCopula(2,300)
X₁ = Exponential()
X₂ = Exponential()
X = SklarDist(C,(X₁,X₂))
Random.rand!(X,matkk')
```
The same overflows appends of course, but here returns zeros and not NaN. If you dont want the overflows, however, moving to `BigFloats` sampling is really easy:
```julia
bigC = ClaytonCopula(2,big(300))
bigX = SklarDist(bigC,(X₁,X₂))
Random.rand!(bigX,matkk')
```
Note that \theta=300 is insanely close to a pure comonotony, especially in the lower tail, which is why your inverse Rosenblat transformation is failling when `u` is too small. You could also just pass bigfloats random numbers to your rosenblat transformation, it’ll work too.
See on [the github page](https://github.com/lrnv/Copulas.jl) for more details: many other copulas are available.