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
89
90
91
92
93
94
|
/*
* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
* Copyright (C) 2008 - INRIA - Arnaud TORSET
*
* This file must be used under the terms of the CeCILL.
* This source file is licensed as described in the file COPYING, which
* you should have received as part of this distribution. The terms
* are also available at
* http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
*
*/
/*
URAND, A UNIVERSAL RANDOM NUMBER GENERATOR
BY, MICHAEL A. MALCOLM, CLEVE B. MOLER,
STAN-CS-73-334, JANUARY 1973,
COMPUTER SCIENCE DEPARTMENT,
School of Humanities and Sciences, STANFORD UNIVERSITY,
ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/73/334/CS-TR-73-334.pdf
*/
#include "rand.h"
#include <stdio.h>
double drands(void) {
int m=1;
const int itwo=2;
static int m2=0,halfm=0,ia=0, ic=0, mic=0, iy=0;
static double s=0.0;
if (m2==0){
/* if first entry, compute machine integer word length */
while (m>m2){
m2=m;
m=itwo*m2;
}
halfm = m2;
/* compute multiplier and increment for linear congruential method */
ia = 8*(int)(halfm*atan(1.0)/8.0) + 5;
ic = 2*(int)(halfm*(0.5-sqrt(3.0)/6.0)) + 1;
mic = (m2 - ic) + m2;
/* s is the scale factor for converting to floating point */
s = 0.5/halfm;
}
/* compute next random number */
iy = iy*ia;
/* the following statement is for computers which do not allow
integer overflow on addition */
if (iy > mic)
{
iy = (iy - m2) - m2;
}
iy = iy + ic;
/* the following statement is for computers where the
word length for addition is greater than for multiplication */
if (iy/2 > m2)
{
iy = (iy - m2) - m2;
}
/* the following statement is for computers where integer
overflow affects the sign bit */
if (iy < 0)
{
iy = (iy + m2) + m2;
}
return (double)iy*s;
}
|
