3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
21 This is a lagged Fibonacci generator which supposedly excellent
22 statistical properties (I do not concur)
24 I got it from the net and translated into C.
26 * ======================================================================
27 * NIST Guide to Available Math Software.
28 * Fullsource for module UNI from package CMLIB.
29 * Retrieved from CAMSUN on Tue Oct 8 14:04:10 1996.
30 * ======================================================================
32 C***BEGIN PROLOGUE UNI
33 C***DATE WRITTEN 810915
34 C***REVISION DATE 830805
35 C***CATEGORY NO. L6A21
36 C***KEYWORDS RANDOM NUMBERS, UNIFORM RANDOM NUMBERS
37 C***AUTHOR BLUE, JAMES, SCIENTIFIC COMPUTING DIVISION, NBS
38 C KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS
39 C MARSAGLIA, GEORGE, COMPUTER SCIENCE DEPT., WASH STATE UNIV
41 C***PURPOSE THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON [0,1
42 C AND CAN BE USED ON ANY COMPUTER WITH WHICH ALLOWS INTEGERS
43 C AT LEAST AS LARGE AS 32767.
46 C THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON THE INTER
47 C [0,1). IT CAN BE USED WITH ANY COMPUTER WHICH ALLOWS
48 C INTEGERS AT LEAST AS LARGE AS 32767.
54 C HERE JD IS ANY N O N - Z E R O INTEGER.
55 C THIS CAUSES INITIALIZATION OF THE PROGRAM
56 C AND THE FIRST RANDOM NUMBER TO BE RETURNED AS Z.
59 C CAUSES THE NEXT RANDOM NUMBER TO BE RETURNED AS Z.
62 C..................................................................
63 C NOTE: USERS WHO WISH TO TRANSPORT THIS PROGRAM FROM ONE COMPUTER
64 C TO ANOTHER SHOULD READ THE FOLLOWING INFORMATION.....
66 C MACHINE DEPENDENCIES...
67 C MDIG = A LOWER BOUND ON THE NUMBER OF BINARY DIGITS AVAILABLE
68 C FOR REPRESENTING INTEGERS, INCLUDING THE SIGN BIT.
69 C THIS VALUE MUST BE AT LEAST 16, BUT MAY BE INCREASED
70 C IN LINE WITH REMARK A BELOW.
73 C A. THIS PROGRAM CAN BE USED IN TWO WAYS:
74 C (1) TO OBTAIN REPEATABLE RESULTS ON DIFFERENT COMPUTERS,
75 C SET 'MDIG' TO THE SMALLEST OF ITS VALUES ON EACH, OR,
76 C (2) TO ALLOW THE LONGEST SEQUENCE OF RANDOM NUMBERS TO BE
77 C GENERATED WITHOUT CYCLING (REPEATING) SET 'MDIG' TO THE
78 C LARGEST POSSIBLE VALUE.
79 C B. THE SEQUENCE OF NUMBERS GENERATED DEPENDS ON THE INITIAL
80 C INPUT 'JD' AS WELL AS THE VALUE OF 'MDIG'.
81 C IF MDIG=16 ONE SHOULD FIND THAT
82 Editors Note: set the seed using 152 in order to get uni(305)
84 C THE FIRST EVALUATION
85 C Z=UNI(305) GIVES Z=.027832881...
86 C THE SECOND EVALUATION
87 C Z=UNI(0) GIVES Z=.56102176...
88 C THE THIRD EVALUATION
89 C Z=UNI(0) GIVES Z=.41456343...
90 C THE THOUSANDTH EVALUATION
91 C Z=UNI(0) GIVES Z=.19797357...
93 C***REFERENCES MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM
94 C NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U.
95 C***ROUTINES CALLED I1MACH,XERROR
102 #include <gsl/gsl_rng.h>
104 static inline unsigned long int uni_get (void *vstate);
105 static double uni_get_double (void *vstate);
106 static void uni_set (void *state, unsigned long int s);
108 static const unsigned int MDIG = 16; /* Machine digits in int */
109 static const unsigned int m1 = 32767; /* 2^(MDIG-1) - 1 */
110 static const unsigned int m2 = 256; /* 2^(MDIG/2) */
119 static inline unsigned long
120 uni_get (void *vstate)
122 uni_state_t *state = (uni_state_t *) vstate;
123 const int i = state->i;
124 const int j = state->j;
126 /* important k not be unsigned */
127 long k = state->m[i] - state->m[j];
155 uni_get_double (void *vstate)
157 return uni_get (vstate) / 32767.0 ;
161 uni_set (void *vstate, unsigned long int s)
163 unsigned int i, seed, k0, k1, j0, j1;
165 uni_state_t *state = (uni_state_t *) vstate;
167 /* For this routine, the seeding is very elaborate! */
168 /* A flaw in this approach is that seeds 1,2 give exactly the
169 same random number sequence! */
171 s = 2 * s + 1; /* enforce seed be odd */
172 seed = (s < m1 ? s : m1); /* seed should be less than m1 */
179 for (i = 0; i < 17; ++i)
182 j1 = (seed / m2 + j0 * k1 + j1 * k0) % (m2 / 2);
184 state->m[i] = j0 + m2 * j1;
192 static const gsl_rng_type uni_type =
194 32766, /* RAND_MAX */
196 sizeof (uni_state_t),
201 const gsl_rng_type *gsl_rng_uni = &uni_type;