Random float in C using getrandom

If you need to generate doubles, the following algorithm could be of use:

CPython generates random numbers using the following algorithm (I changed the function name, typedefs and return values, but algorithm remains the same):

double get_random_double() {
    uint32_t a = get_random_uint32_t() >> 5;
    uint32_t b = get_random_uint32_t() >> 6;
    return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
}

The source of that algorithm is a Mersenne Twister 19937 random number generator by Takuji Nishimura and Makoto Matsumoto. Unfortunately the original link mentioned in the source is not available for download any longer.

The comment on this function in CPython notes the following:

[this function] is the function named genrand_res53 in the original code; generates a random number on [0,1) with 53-bit resolution; note that 9007199254740992 == 2**53; I assume they're spelling "/2**53" as multiply-by-reciprocal in the (likely vain) hope that the compiler will optimize the division away at compile-time. 67108864 is 2**26. In effect, a contains 27 random bits shifted left 26, and b fills in the lower 26 bits of the 53-bit numerator.

The orginal code credited Isaku Wada for this algorithm, 2002/01/09

---

Simplifying from that code, if you want to create a float fast, you should mask the bits of uint32_t with (1 << FLT_MANT_DIG) - 1 and divide by (1 << FLT_MANT_DIG) to get the proper [0, 1) interval:

#include <stdio.h>
#include <sys/syscall.h>
#include <unistd.h>
#include <stdint.h>
#include <float.h>

int main() {
    uint32_t r = 0;
    float result;
    for (int i = 0; i < 20; i++) {
        syscall(SYS_getrandom, &r, sizeof(uint32_t), 0);
        result = (float)(r & ((1 << FLT_MANT_DIG) - 1)) / (1 << FLT_MANT_DIG);
        printf("%f\n", result);
    }
    return 0;
}

Since it can be assumed that your Linux has a C99 compiler, we can use ldexpf instead of that division:

#include <math.h>

result = ldexpf(r & ((1 << FLT_MANT_DIG) - 1), -FLT_MANT_DIG);

To get the closed interval [0, 1], you can do the slightly less efficient

result = ldexpf(r % (1 << FLT_MANT_DIG), -FLT_MANT_DIG);

---

To generate lots of good quality random numbers fast, I'd just use the system call to fetch enough data to seed a PRNG or CPRNG, and proceed from there.

OP has 2 issues:

1. How to started the sequence very randomly.

2. How to generate a double on the [0...1) range.

The usual method is to take a very random source like /dev/urandom or the result from the syscall() or maybe even seed = time() ^ process_id; and seed via srand(). Then call rand() as needed.

Below includes a quickly turned method to generate a uniform [0.0 to 1.0) (linear distribution). But like all random generating functions, really good ones are base on extensive study. This one simply calls rand() a few times based on DBL_MANT_DIG and RAND_MAX,

[Edit] Original double rand_01(void) has a weakness in that it only generates a 2^52 different doubles rather than 2^53. It has been amended. Alternative: a double version of rand_01_ld(void) far below.

#include <assert.h>
#include <float.h>
#include <limits.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

double rand_01(void) {
  assert(FLT_RADIX == 2); // needed for DBL_MANT_DIG
  unsigned long long limit = (1ull << DBL_MANT_DIG) - 1;
  double r = 0.0;
  do {
    r += rand();
    // Assume RAND_MAX is a power-of-2 - 1
    r /= (RAND_MAX/2 + 1)*2.0;
    limit = limit / (RAND_MAX/2 + 1) / 2;
  } while (limit);

  // Use only DBL_MANT_DIG (53) bits of precision.
  if (r < 0.5) {
    volatile double sum = 0.5 + r;
    r = sum - 0.5;
  }
  return r;
}

int main(void) {
  FILE *istream = fopen("/dev/urandom", "rb");
  assert(istream);
  unsigned long seed = 0;
  for (unsigned i = 0; i < sizeof seed; i++) {
    seed *= (UCHAR_MAX + 1);
    int ch = fgetc(istream);
    assert(ch != EOF);
    seed += (unsigned) ch;
  }
  fclose(istream);
  srand(seed);

  for (int i=0; i<20; i++) {
    printf("%f\n", rand_01());
  }

  return 0;
}

---

If one wanted to extend to an even wider FP, unsigned wide integer types may be insufficient. Below is a portable method that does not have that limitation.

long double rand_01_ld(void) {
  // These should be calculated once rather than each function call
  // Leave that as a separate implementation problem
  // Assume RAND_MAX is power-of-2 - 1
  assert((RAND_MAX & (RAND_MAX + 1U)) == 0);
  double rand_max_p1 = (RAND_MAX/2 + 1)*2.0;
  unsigned BitsPerRand = (unsigned) round(log2(rand_max_p1));
  assert(FLT_RADIX != 10);
  unsigned BitsPerFP = (unsigned) round(log2(FLT_RADIX)*LDBL_MANT_DIG);

  long double r = 0.0;
  unsigned i;
  for (i = BitsPerFP; i >= BitsPerRand; i -= BitsPerRand) {
    r += rand();
    r /= rand_max_p1;
  }
  if (i) {
    r += rand() % (1 << i);
    r /= 1 << i;
  }
  return r;
}