]> hackdaworld.org Git - physik/posic.git/blob - random/random.c
tersoff MAYBE fixed now!
[physik/posic.git] / 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=malloc(RAND_BUFSIZE*sizeof(unsigned int));
40         if(random->buffer==NULL) {
41                 perror("malloc random buffer");
42                 return -1;
43         }
44         random->b_ptr=random->buffer+RAND_BUFSIZE;
45
46         dprintf(random->logfd,"[random] init successfull\n");
47
48         return 0;
49 }
50
51
52 int rand_close(t_random *random) {
53
54         dprintf(random->logfd,"[random] shutdown\n");
55
56         if(random->buffer) free(random->buffer);
57         if(random->fd) close(random->fd);
58         if(random->logfd) close(random->logfd);
59
60         return 0;
61 }
62
63 unsigned int rand_get(t_random *random) {
64
65         int left;
66
67         if(random->b_ptr==random->buffer+RAND_BUFSIZE) {
68                 if(random->status&RAND_STAT_VERBOSE)
69                         dprintf(random->logfd,
70                                 "[random] getting new random numbers\n");
71                 if(!(random->status&RAND_STAT_UDEV)) {
72                         lseek(random->fd,0,SEEK_SET);
73                         dprintf(random->logfd,
74                                 "[random] warning, rereading random file\n");
75                 }
76                 left=RAND_BUFSIZE*sizeof(unsigned int);
77                 while(left) {
78                         left-=read(random->fd,
79                                    random->buffer+\
80                                    RAND_BUFSIZE*sizeof(unsigned int)-left,
81                                    left);
82                 }
83                 if(random->status&RAND_STAT_VERBOSE)
84                         dprintf(random->logfd,
85                                 "[random] got new random numbers\n");
86                 random->b_ptr=random->buffer;
87         }
88
89         return(*(random->b_ptr++));
90 }
91
92 double rand_get_double(t_random *random) {
93
94         return(1.0*rand_get(random)/((long long unsigned int)URAND_MAX+1));
95 }
96
97 double rand_get_double_linear(t_random *random,double a,double b) {
98
99         return((-b+sqrt(b*b+2*a*(b+a/2)*rand_get_double(random)))/a);
100 }
101
102 double rand_get_gauss(t_random *random) {
103
104         double a,b,w;
105
106         if(random->status&RAND_STAT_GAUSS) {
107                 random->status&=~RAND_STAT_GAUSS;
108                 return random->gauss;
109         }
110
111         w=0;
112         while((w>=1.0)||(w==0.0)) {
113                 a=-2.0*rand_get_double(random)+1.0;
114                 b=-2.0*rand_get_double(random)+1.0;
115                 w=a*a+b*b;
116         }
117
118         w=sqrt(-2.0*log(w)/w);
119
120         random->gauss=b*w;
121         random->status|=RAND_STAT_GAUSS;
122         
123         return(a*w);
124 }
125         
126 unsigned int rand_get_int_linear(t_random *random,unsigned int max,
127                              double a,double b) {
128
129         return((int)(rand_get_double_linear(random,a,b)*max));
130 }       
131
132 unsigned int rand_get_reject(t_random *random,unsigned int max_x,
133                              unsigned int max_y,unsigned int *graph) {
134
135         unsigned int x,y;
136         unsigned char ok;
137
138         ok=0;
139         while(!ok) {
140                 x=(max_x*rand_get_double(random));
141                 y=(max_y*rand_get_double(random));
142                 if(y<=graph[x]) ok=1;
143         }
144
145         return x;
146 }
147