-
[physik/nlsop.git] / random.c
1 /*
2  * random functions
3  *
4  * author: hackbard@hackdaworld.dyndns.org
5  *
6  */
7
8 #include <sys/types.h>
9 #include <sys/stat.h>
10 #include <fcntl.h>
11 #include <unistd.h>
12 #include <stdlib.h>
13 #include <math.h>
14
15 #include "random.h"
16
17 static int rand_fd;
18 static u32 *c_ptr,*b_ptr;
19
20 int rand_init(char *rf)
21 {
22  if(rf==NULL)
23  {
24   if((rand_fd=open("/dev/urandom",O_RDONLY))<0)
25   {
26    puts("cannot open /dev/urandom");
27    return -1;
28   }
29  } else
30  {
31   if((rand_fd=open(rf,O_RDONLY))<0)
32   {
33    printf("cannot open %s\n",rf);
34    return -1;
35   }
36  }
37  if((b_ptr=(u32 *)(malloc(BUFSIZE*sizeof(u32))))==NULL)
38  {
39   puts("failed allocating random buffer");
40   return -1;
41  }
42  c_ptr=b_ptr+BUFSIZE;
43  
44  return 1;
45 }
46
47 int rand_close(void)
48 {
49  close(rand_fd);
50  
51  return 1;
52 }
53
54 u32 get_rand(u32 max)
55 {
56  if(c_ptr>=b_ptr+BUFSIZE)
57  {
58   if(read(rand_fd,b_ptr,BUFSIZE*sizeof(u32))<BUFSIZE*sizeof(u32))
59   {
60    /* -> assume random file, end reached */
61    puts("warning, will read random file from beginning ...");
62    lseek(rand_fd,0,SEEK_SET);
63    read(rand_fd,b_ptr,BUFSIZE*sizeof(u32));
64   }
65   c_ptr=b_ptr;
66  }
67
68  return((u32)(*(c_ptr++)*(max*1.0/((long long unsigned int)URAND_MAX+1))));
69 }
70
71 u32 get_rand_lgp(u32 max,double a,double b)
72 {
73  return((u32)(1.0*max*(-1.0*b+sqrt(b*b+2*a*((b+a/2)*get_rand(URAND_MAX)/((long long unsigned int)URAND_MAX+1))))/a));
74 }
75
76 u32 get_rand_reject(u32 max_x,u32 max_y,u32 *graph)
77 {
78  u32 x,y;
79  unsigned char ok;
80
81  ok=0;
82  while(!ok)
83  {
84   x=get_rand(max_x);
85   y=get_rand(max_y);
86   if(y<=graph[x]) ok=1;
87  }
88
89  return x;
90 }