How to generate a uniformly distributed random double in C using libsodium in range [-a,a]?

Question

The libsodium library has a function

uint32_t randombytes_uniform(const uint32_t upper_bound);

but obviously this returns an unsigned integer. Can I somehow use this to generate a uniformly distributed random double in range [-a,a] where a is also a double given by the user ? I am especially focused on the result being uniformly distributed/unbiased, so that is why I would like to use the libsodium library.

Answer 1

const uint32_t mybound = 1000000000; // Example
const uint32_t x = randombytes_uniform(mybound);
const double a = 3.5; // Example
const double variate = a * ( (2.0 * x / mybound) - 1);

Answer 2

Let me try to to do it step-by-step.

First, you obviously need to combine two calls to get up to 64bit of randomness for one double value output.

Second, you convert it to [0...1] interval. There are several ways to do it, all of the are good in some sense or another, I prefer uniform random dyadic rationals in the form n*2^-53, see here for details. You could try other methods listed above as well. NB: methods in the link produce results in [0...1) range, I've tried to do acceptance/rejection to get closed [0...1] range.

Last, I scale result into desired range.

Sorry, C++ only but it is trivial to convert to C

#include <stdint.h>
#include <math.h>
#include <iostream>
#include <random>

// emulate libsodium RNG, valid for full 32bits result only!
static uint32_t randombytes_uniform(const uint32_t upper_bound) {
    static std::mt19937 mt{9876713};
    return mt();
}

// get 64bits from two 32bit numbers
static inline uint64_t rng() {
    return (uint64_t)randombytes_uniform(UINT32_MAX) << 32 | randombytes_uniform(UINT32_MAX);
}

const  int32_t bits_in_mantissa = 53;
const uint64_t max  = (1ULL << bits_in_mantissa);
const uint64_t mask = (1ULL << (bits_in_mantissa+1)) - 1;

static double rnd(double a, double b) {
    uint64_t r;
    do {
        r = rng() & mask; // get 54 random bits, need 53 or max
    } while (r > max);
    double v = ldexp( (double)r, -bits_in_mantissa ); // http://xoshiro.di.unimi.it/random_real.c
    return a + (b-a)*v;
}

int main() {
    double a = -3.5;
    double b = 3.5;

    for(int k = 0; k != 100; ++k)
        std::cout << rnd(a, b) << '\n';

    return 0;
}

Answer 3

First recognizing that finding a random number [0...a] is a sufficient step, followed by a coin flip for +/-.

Step 2. Find the expo such that a < 2**expo or ceil(log2(a)).

int sign;
do {
  int exp;
  frexp(a, &exp);

Step 3. Form an integral 63-bit random number [0...0x7FFF_FFFF_FFFF_FFFF] and random sign. The 63 should be at least as wide as the precision of a double - which is often 53 bits. At this point r is certainly uniform.

  unit64_t r = randombytes_uniform(0xFFFFFFFF);
  r <<= 32;
  r |= randombytes_uniform(0xFFFFFFFF);
  // peel off one bit for sign
  sign = r & 1;
  r >>= 1; 

Step 4. Scale and test if in range. Repeat as needed.

  double candidate = ldexp(r/pow(2 63), expo);
} while (candidate > a);

Step 5. Apply the sign.

if (sign) {
  candidate = -candidate;
}
return candidate;

---

Avoid (2.0 * x / a) - 1 as the calculation is not symmetric about 0.0.

---

Code would benefit with improvements to deal with a near DBL_MAX.

Some rounding issues apply that this answer glosses over, yet the distribution remains uniform - except potentially at the edges.