/*
   Generates random uniform, normal deviates by reading random bytes
   from /dev/random.  When block is exhausted, refill. rand_data
   is of type uint32_t to ensure that 4 bytes are available for
   filling, which uses most of the precision available with doubles.
   Use built-in C standard library constant BUFSIZ for buffer size.
   This can be customized to the system for added speed, most likely
   by making it larger.

   Written by David Smith, 5 Jul 2006. 
*/

#include <stdio.h>
#include <stdlib.h>
#include <err.h>
#include <fcntl.h>
#include <stdint.h>
#include <sysexits.h>
#include <unistd.h>

double unidev()
{
   /*   Read random data from /dev/random into a block of memory
      and dole out random numbers as needed.  When block is
      exhausted, refill. rdata is of type uint32_t to ensure
      that 4 bytes are available for filling, which uses 
      most of the precision available with doubles.  */
   static uint32_t rdata[BUFSIZ]; 
   static int fd = -1, loc = BUFSIZ;
   static double factor;
   const char *rdev = "/dev/urandom";

   if (fd < 0) 
   {
      if ((fd = open(rdev, O_RDONLY)) < 0) 
         err(EX_NOINPUT, "%s open", rdev);
      factor = 1. / (UINT32_MAX + 1.);  /* range (0,1] */
      /* for a range of [0,1], factor should be just 1/UINT32_MAX */
   }

   if (BUFSIZ <= loc)   
   {
      loc = 0;   
      if (read(fd, &rdata, 4*BUFSIZ) != 4*BUFSIZ) 
         err(EX_NOINPUT, "%s read", rdev);
   }

   return 1. - rdata[loc++]*factor; 
}