-
Notifications
You must be signed in to change notification settings - Fork 0
/
gauss.c
88 lines (77 loc) · 1.63 KB
/
gauss.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include <stdlib.h>
#include <math.h>
#include <stddef.h>
#define REAL double
#define SEQUENCES 10000000
#ifndef RAND_MAX
#define RAND_MAX 2147483647
#endif
/************************************
SEQUENCES is the maximum number of random sequences
that will be produced. A sequence changes every second
RAND_MAX is the factor that keeps the numbers
between 0.0 and 1.0. If the numbers get outside
of this interval you'll have to change RAND_MAX.
Look in the "man random" about what is the maximum
number generated by random() and set RAND_MAX to it
*************************************/
//int seed=3765; //for randm_() function
void seed_()
/*
* Seeds the rnd_() function
*/
{
// REAL r;
// randm_(&seed,&r);
srandom(time(0) % SEQUENCES);
// srand(time(0) % SEQUENCES);
}
REAL rnd_()
/*
* Returns a random number between 0.0 and 1.0
*/
{
// REAL r;
// randm_(&seed,&r);
// return r;
return (REAL)random()/RAND_MAX;
// return (REAL)rand()/RAND_MAX;
}
void gaussn_(REAL *Y, REAL d, int n)
/* Generates an array of n random numbers with a
* Normal (Gaussian) distribution with a standard
* deviation of d, and stores them in array Y.
*/
{
extern REAL gauss_();
REAL *x;
for (x = Y; x - Y < n; x++)
*x = d*gauss_();
}
REAL gauss_()
/* From Numerical Recepes in C, Ch.7.2
* Zero mean and unit variance
*/
{
extern REAL rnd_();
static int iset = 0;
static REAL gset;
REAL fac,rsq,v1,v2;
if (iset == 0)
{ do
{
v1=2.*rnd_()-1.;
v2=2.*rnd_()-1.;
rsq=v1*v1+v2*v2;
} while (rsq >= 1.0 || rsq == 0);
fac=(REAL)sqrt(-2.*log(rsq)/rsq);
gset=v1*fac;
iset=1;
return v2*fac;
}
else
{
iset=0;
return gset;
}
}