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.
23 #include <gsl/gsl_rng.h>
25 /* This is the CRAY RANF generator. The generator returns the
26 upper 32 bits from each term of the sequence,
28 x_{n+1} = (a x_n) mod m
30 using 48-bit unsigned arithmetic, with a = 0x2875A2E7B175 and m =
31 2^48. The seed specifies the lower 32 bits of the initial value,
32 x_1, with the lowest bit set (to prevent the seed taking an even
33 value), and the upper 16 bits set to 0.
35 There is a subtlety in the implementation of the seed. The initial
36 state is put one step back by multiplying by the modular inverse of
37 a mod m. This is done for compatibility with the original CRAY
40 Note, you can only seed the generator with integers up to 2^32,
41 while the CRAY uses wide integers which can cover all 2^48 states
44 The theoretical value of x_{10001} is 141091827447341.
46 The period of this generator is 2^{46}. */
48 static inline void ranf_advance (void *vstate);
49 static unsigned long int ranf_get (void *vstate);
50 static double ranf_get_double (void *vstate);
51 static void ranf_set (void *state, unsigned long int s);
53 static const unsigned short int a0 = 0xB175 ;
54 static const unsigned short int a1 = 0xA2E7 ;
55 static const unsigned short int a2 = 0x2875 ;
59 unsigned short int x0, x1, x2;
64 ranf_advance (void *vstate)
66 ranf_state_t *state = (ranf_state_t *) vstate;
68 const unsigned long int x0 = (unsigned long int) state->x0 ;
69 const unsigned long int x1 = (unsigned long int) state->x1 ;
70 const unsigned long int x2 = (unsigned long int) state->x2 ;
75 state->x0 = (r & 0xFFFF) ;
78 r += a0 * x1 + a1 * x0 ;
79 state->x1 = (r & 0xFFFF) ;
82 r += a0 * x2 + a1 * x1 + a2 * x0 ;
83 state->x2 = (r & 0xFFFF) ;
86 static unsigned long int
87 ranf_get (void *vstate)
89 unsigned long int x1, x2;
91 ranf_state_t *state = (ranf_state_t *) vstate;
92 ranf_advance (state) ;
94 x1 = (unsigned long int) state->x1;
95 x2 = (unsigned long int) state->x2;
97 return (x2 << 16) + x1;
101 ranf_get_double (void * vstate)
103 ranf_state_t *state = (ranf_state_t *) vstate;
105 ranf_advance (state) ;
107 return (ldexp((double) state->x2, -16)
108 + ldexp((double) state->x1, -32)
109 + ldexp((double) state->x0, -48)) ;
113 ranf_set (void *vstate, unsigned long int s)
115 ranf_state_t *state = (ranf_state_t *) vstate;
117 unsigned short int x0, x1, x2 ;
118 unsigned long int r ;
120 unsigned long int b0 = 0xD6DD ;
121 unsigned long int b1 = 0xB894 ;
122 unsigned long int b2 = 0x5CEE ;
124 if (s == 0) /* default seed */
132 x0 = (s | 1) & 0xFFFF ;
133 x1 = s >> 16 & 0xFFFF ;
138 state->x0 = (r & 0xFFFF) ;
141 r += b0 * x1 + b1 * x0 ;
142 state->x1 = (r & 0xFFFF) ;
145 r += b0 * x2 + b1 * x1 + b2 * x0 ;
146 state->x2 = (r & 0xFFFF) ;
151 static const gsl_rng_type ranf_type =
153 0xffffffffUL, /* RAND_MAX */
155 sizeof (ranf_state_t),
161 const gsl_rng_type *gsl_rng_ranf = &ranf_type;