random/random.c

   1 /*
   2  * random.c - basic random deviates
   3  *
   4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
   5  *
   6  */
   7
   8 #define _GNU_SOURCE
   9 #include <stdio.h>
  10 #include <sys/types.h>
  11 #include <sys/stat.h>
  12 #include <fcntl.h>
  13 #include <unistd.h>
  14 #include <stdlib.h>
  15 #include <math.h>
  16
  17 #include "random.h"
  18
  19 int rand_init(t_random *random,char *randomfile,int logfd) {
  20
  21         random->status=0;
  22
  23         if(logfd) random->logfd=logfd;
  24
  25         if(randomfile==NULL) {
  26                 random->fd=open("/dev/urandom",O_RDONLY);
  27                 if(random->fd<0) {
  28                         perror("open urandom");
  29                         return -1;
  30                 }
  31                 random->status|=RAND_STAT_UDEV;
  32         } else {
  33                 random->fd=open(randomfile,O_RDONLY);
  34                 if(random->fd<0) {
  35                         perror("open randomfile");
  36                         return -1;
  37                 }
  38         }
  39         random->buffer=NULL;
  40         random->buffer=malloc(RAND_BUFSIZE*sizeof(unsigned int));
  41         if(random->buffer==NULL) {
  42                 perror("malloc random buffer");
  43                 return -1;
  44         }
  45         random->b_ptr=random->buffer+RAND_BUFSIZE;
  46
  47         dprintf(random->logfd,"[random] init successfull\n");
  48
  49         return 0;
  50 }
  51
  52
  53 int rand_close(t_random *random) {
  54
  55         dprintf(random->logfd,"[random] shutdown\n");
  56
  57         if(random->buffer) free(random->buffer);
  58         if(random->fd) close(random->fd);
  59         if(random->logfd>2) close(random->logfd); // could be stdo/e
  60
  61         return 0;
  62 }
  63
  64 unsigned int rand_get(t_random *random) {
  65
  66         int left;
  67
  68         if(random->b_ptr==random->buffer+RAND_BUFSIZE) {
  69                 if(random->status&RAND_STAT_VERBOSE)
  70                         dprintf(random->logfd,
  71                                 "[random] getting new random numbers\n");
  72                 if(!(random->status&RAND_STAT_UDEV)) {
  73                         lseek(random->fd,0,SEEK_SET);
  74                         dprintf(random->logfd,
  75                                 "[random] warning, rereading random file\n");
  76                 }
  77                 left=RAND_BUFSIZE*sizeof(unsigned int);
  78                 while(left) {
  79                         left-=read(random->fd,
  80                                    random->buffer+\
  81                                    RAND_BUFSIZE*sizeof(unsigned int)-left,
  82                                    left);
  83                 }
  84                 if(random->status&RAND_STAT_VERBOSE)
  85                         dprintf(random->logfd,
  86                                 "[random] got new random numbers\n");
  87                 random->b_ptr=random->buffer;
  88         }
  89
  90         return(*(random->b_ptr++));
  91 }
  92
  93 double rand_get_double(t_random *random) {
  94
  95         return(1.0*rand_get(random)/((long long unsigned int)URAND_MAX+1));
  96 }
  97
  98 double rand_get_double_linear(t_random *random,double a,double b) {
  99
 100         return((-b+sqrt(b*b+2*a*(b+a/2)*rand_get_double(random)))/a);
 101 }
 102
 103 double rand_get_gauss(t_random *random) {
 104
 105         double a,b,w;
 106
 107         if(random->status&RAND_STAT_GAUSS) {
 108                 random->status&=~RAND_STAT_GAUSS;
 109                 return random->gauss;
 110         }
 111
 112         w=0;
 113         while((w>=1.0)||(w==0.0)) {
 114                 a=-2.0*rand_get_double(random)+1.0;
 115                 b=-2.0*rand_get_double(random)+1.0;
 116                 w=a*a+b*b;
 117         }
 118
 119         w=sqrt(-2.0*log(w)/w);
 120
 121         random->gauss=b*w;
 122         random->status|=RAND_STAT_GAUSS;
 123
 124         return(a*w);
 125 }
 126
 127 unsigned int rand_get_int_linear(t_random *random,unsigned int max,
 128                              double a,double b) {
 129
 130         return((int)(rand_get_double_linear(random,a,b)*max));
 131 }
 132
 133 unsigned int rand_get_reject(t_random *random,unsigned int max_x,
 134                              unsigned int max_y,unsigned int *graph) {
 135
 136         unsigned int x,y;
 137         unsigned char ok;
 138
 139         ok=0;
 140         while(!ok) {
 141                 x=(max_x*rand_get_double(random));
 142                 y=(max_y*rand_get_double(random));
 143                 if(y<=graph[x]) ok=1;
 144         }
 145
 146         return x;
 147 }
 148