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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
|
/* rng/uni.c
*
* Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/**
This is a lagged Fibonacci generator which supposedly excellent
statistical properties (I do not concur)
I got it from the net and translated into C.
* ======================================================================
* NIST Guide to Available Math Software.
* Fullsource for module UNI from package CMLIB.
* Retrieved from CAMSUN on Tue Oct 8 14:04:10 1996.
* ======================================================================
C***BEGIN PROLOGUE UNI
C***DATE WRITTEN 810915
C***REVISION DATE 830805
C***CATEGORY NO. L6A21
C***KEYWORDS RANDOM NUMBERS, UNIFORM RANDOM NUMBERS
C***AUTHOR BLUE, JAMES, SCIENTIFIC COMPUTING DIVISION, NBS
C KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS
C MARSAGLIA, GEORGE, COMPUTER SCIENCE DEPT., WASH STATE UNIV
C
C***PURPOSE THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON [0,1
C AND CAN BE USED ON ANY COMPUTER WITH WHICH ALLOWS INTEGERS
C AT LEAST AS LARGE AS 32767.
C***DESCRIPTION
C
C THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON THE INTER
C [0,1). IT CAN BE USED WITH ANY COMPUTER WHICH ALLOWS
C INTEGERS AT LEAST AS LARGE AS 32767.
C
C
C USE
C FIRST TIME....
C Z = UNI(JD)
C HERE JD IS ANY N O N - Z E R O INTEGER.
C THIS CAUSES INITIALIZATION OF THE PROGRAM
C AND THE FIRST RANDOM NUMBER TO BE RETURNED AS Z.
C SUBSEQUENT TIMES...
C Z = UNI(0)
C CAUSES THE NEXT RANDOM NUMBER TO BE RETURNED AS Z.
C
C
C..................................................................
C NOTE: USERS WHO WISH TO TRANSPORT THIS PROGRAM FROM ONE COMPUTER
C TO ANOTHER SHOULD READ THE FOLLOWING INFORMATION.....
C
C MACHINE DEPENDENCIES...
C MDIG = A LOWER BOUND ON THE NUMBER OF BINARY DIGITS AVAILABLE
C FOR REPRESENTING INTEGERS, INCLUDING THE SIGN BIT.
C THIS VALUE MUST BE AT LEAST 16, BUT MAY BE INCREASED
C IN LINE WITH REMARK A BELOW.
C
C REMARKS...
C A. THIS PROGRAM CAN BE USED IN TWO WAYS:
C (1) TO OBTAIN REPEATABLE RESULTS ON DIFFERENT COMPUTERS,
C SET 'MDIG' TO THE SMALLEST OF ITS VALUES ON EACH, OR,
C (2) TO ALLOW THE LONGEST SEQUENCE OF RANDOM NUMBERS TO BE
C GENERATED WITHOUT CYCLING (REPEATING) SET 'MDIG' TO THE
C LARGEST POSSIBLE VALUE.
C B. THE SEQUENCE OF NUMBERS GENERATED DEPENDS ON THE INITIAL
C INPUT 'JD' AS WELL AS THE VALUE OF 'MDIG'.
C IF MDIG=16 ONE SHOULD FIND THAT
Editors Note: set the seed using 152 in order to get uni(305)
-jt
C THE FIRST EVALUATION
C Z=UNI(305) GIVES Z=.027832881...
C THE SECOND EVALUATION
C Z=UNI(0) GIVES Z=.56102176...
C THE THIRD EVALUATION
C Z=UNI(0) GIVES Z=.41456343...
C THE THOUSANDTH EVALUATION
C Z=UNI(0) GIVES Z=.19797357...
C
C***REFERENCES MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM
C NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U.
C***ROUTINES CALLED I1MACH,XERROR
C***END PROLOGUE UNI
**/
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
static inline unsigned long int uni_get (void *vstate);
static double uni_get_double (void *vstate);
static void uni_set (void *state, unsigned long int s);
static const unsigned int MDIG = 16; /* Machine digits in int */
static const unsigned int m1 = 32767; /* 2^(MDIG-1) - 1 */
static const unsigned int m2 = 256; /* 2^(MDIG/2) */
typedef struct
{
int i, j;
unsigned long m[17];
}
uni_state_t;
static inline unsigned long
uni_get (void *vstate)
{
uni_state_t *state = (uni_state_t *) vstate;
const int i = state->i;
const int j = state->j;
/* important k not be unsigned */
long k = state->m[i] - state->m[j];
if (k < 0)
k += m1;
state->m[j] = k;
if (i == 0)
{
state->i = 16;
}
else
{
(state->i)--;
}
if (j == 0)
{
state->j = 16;
}
else
{
(state->j)--;
}
return k;
}
static double
uni_get_double (void *vstate)
{
return uni_get (vstate) / 32767.0 ;
}
static void
uni_set (void *vstate, unsigned long int s)
{
unsigned int i, seed, k0, k1, j0, j1;
uni_state_t *state = (uni_state_t *) vstate;
/* For this routine, the seeding is very elaborate! */
/* A flaw in this approach is that seeds 1,2 give exactly the
same random number sequence! */
s = 2 * s + 1; /* enforce seed be odd */
seed = (s < m1 ? s : m1); /* seed should be less than m1 */
k0 = 9069 % m2;
k1 = 9069 / m2;
j0 = seed % m2;
j1 = seed / m2;
for (i = 0; i < 17; ++i)
{
seed = j0 * k0;
j1 = (seed / m2 + j0 * k1 + j1 * k0) % (m2 / 2);
j0 = seed % m2;
state->m[i] = j0 + m2 * j1;
}
state->i = 4;
state->j = 16;
return;
}
static const gsl_rng_type uni_type =
{"uni", /* name */
32766, /* RAND_MAX */
0, /* RAND_MIN */
sizeof (uni_state_t),
&uni_set,
&uni_get,
&uni_get_double};
const gsl_rng_type *gsl_rng_uni = &uni_type;
|
