2 * random.c - basic random deviates
4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
10 #include <sys/types.h>
19 int rand_init(t_random *random,char *randomfile,int logfd) {
23 if(logfd) random->logfd=logfd;
25 if(randomfile==NULL) {
26 random->fd=open("/dev/urandom",O_RDONLY);
28 perror("open urandom");
31 random->status|=RAND_STAT_UDEV;
33 random->fd=open(randomfile,O_RDONLY);
35 perror("open randomfile");
40 random->buffer=malloc(RAND_BUFSIZE*sizeof(unsigned int));
41 if(random->buffer==NULL) {
42 perror("malloc random buffer");
45 random->b_ptr=random->buffer+RAND_BUFSIZE;
47 dprintf(random->logfd,"[random] init successfull\n");
53 int rand_close(t_random *random) {
55 dprintf(random->logfd,"[random] shutdown\n");
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
64 unsigned int rand_get(t_random *random) {
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");
77 left=RAND_BUFSIZE*sizeof(unsigned int);
79 left-=read(random->fd,
81 RAND_BUFSIZE*sizeof(unsigned int)-left,
84 if(random->status&RAND_STAT_VERBOSE)
85 dprintf(random->logfd,
86 "[random] got new random numbers\n");
87 random->b_ptr=random->buffer;
90 return(*(random->b_ptr++));
93 double rand_get_double(t_random *random) {
95 return(1.0*rand_get(random)/((long long unsigned int)URAND_MAX+1));
98 double rand_get_double_linear(t_random *random,double a,double b) {
100 return((-b+sqrt(b*b+2*a*(b+a/2)*rand_get_double(random)))/a);
103 double rand_get_gauss(t_random *random) {
107 if(random->status&RAND_STAT_GAUSS) {
108 random->status&=~RAND_STAT_GAUSS;
109 return random->gauss;
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;
119 w=sqrt(-2.0*log(w)/w);
122 random->status|=RAND_STAT_GAUSS;
127 unsigned int rand_get_int_linear(t_random *random,unsigned int max,
130 return((int)(rand_get_double_linear(random,a,b)*max));
133 unsigned int rand_get_reject(t_random *random,unsigned int max_x,
134 unsigned int max_y,unsigned int *graph) {
141 x=(max_x*rand_get_double(random));
142 y=(max_y*rand_get_double(random));
143 if(y<=graph[x]) ok=1;