# Using random bits to generate random numbers

Problem: How do you randomly generate any possible floating point number, including NAs and Infs?

Details:

• I wanted everything from the super small 3.1e-15 to the super large 6.6e230
• I wanted Inf and NA to be possible

Thanks to Emil Hvitfeldt and Michael Chirico for their assistance, ideas and feedback!

# Attempt 1 - runif()

My first attempt using runif() to generate any possible number between 0 and the maximum double was never going to work:

rand_uniform <- function(N) {
runif(N, min = 0, max = .Machine\$double.xmax)
}

rand_uniform(10)
  4.773031e+307 6.689646e+307 1.029815e+308 1.632679e+308 3.625622e+307
 1.615029e+308 1.698236e+308 1.187912e+308 1.130954e+308 1.110728e+307

Given that the maximum double value on R is approx 2x10^308:

• 50% of all numbers in the range [0, max_double] are above 1x10^308
• 99% of all numbers in the range [0, max_double] are above 2x10^306

So 99.99% of all numbers sampled from a uniform distribution in the range [0, max_double] will be above 2x10^304.

This is never going to give me the diversity of floating point numbers I wanted.

Also, runif() will never generate an NA or Inf value.

# Attempt 2 - random magnitude and significand in base 10

Separately draw random numbers for the leading part and exponential part of the double

rand_power10 <- function(N) {
runif(N, -1, 1) * 10^runif(N, -308, 308)
}

rand_power10(10)
  -3.530508e+267 -3.094719e-178   1.009171e+93 -5.093762e-232  2.188117e-144
  -3.229847e-73  7.712072e-301   3.498596e-73 -1.287197e+227   2.507201e-99
• Gives a better spread of random doubles
• Not sure that it will ever generate NA or Inf
• Misses out some numbers - e.g. 1.5 x 10^308 will never be output by this method

# Attempt 3 - random magnitude and significand in base 2

Separately draw random numbers for the leading part and exponential part of the double, but do this as powers of 2, no powers of 10

rand_power2 <- function(N) {
runif(N, -1, 1) * 2^runif(N, -1024, 1024)
}

rand_power2(10)
  -2.635011e+196   9.186218e+89 -3.482499e+172  -3.126403e+32   1.375545e+18
  8.286859e+177  7.923594e-295  -7.186193e-15  7.487349e+142 -1.174630e+118
• Gives a great spread of random doubles
• I think this covers the entire range of values
• Not sure that it will ever generate NA or Inf

# Attempt 4 - Generate 64 random bits - 1 bit at atime

Generate any possible 64-bit pattern and cast it to a double precision floating points value.

rand_bits1 <- function(N) {
packBits(
sample(as.raw(c(0, 1)), 64 * N, replace = TRUE),
type = 'raw'
),
what = 'double',
n    = N
)
}

rand_bits1(10)
  -1.764955e-298 -2.887630e+176 -1.333649e+299  3.632060e-292   5.835171e-46
 -1.699497e-298  1.461211e-180 -1.370493e-195  5.631333e-256 -1.145778e-215
• Guaranteed to be able to generate every possible 64bit number
• Can generate NA and Inf

# Attempt 5 - Generate 64 random bits - 8 bits at a time

Generate any possible 64-bit pattern and cast it to a double precision floating points value.

rand_bits8 <- function(N) {
sample(as.raw(0:255), N * 8, replace = TRUE),
what = 'double',
n    = N
)
}

rand_bits8(10)
  -5.384031e+226  -4.136980e+04  4.671568e+278 -1.640025e-283   1.392010e-49
  7.343209e+243  -3.141708e-65 -2.390032e+198 -2.806815e-299  3.123026e+131
• Guaranteed to be able to generate every possible 64bit number
• Can generate NA and Inf

# Benchmark

Quick benchmark of the various methods 