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_sys.h>
24 #include <gsl/gsl_rng.h>
26 /* This is the Unix rand48() generator. The generator returns the
27 upper 32 bits from each term of the sequence,
29 x_{n+1} = (a x_n + c) mod m
31 using 48-bit unsigned arithmetic, with a = 0x5DEECE66D , c = 0xB
32 and m = 2^48. The seed specifies the upper 32 bits of the initial
33 value, x_1, with the lower 16 bits set to 0x330E.
35 The theoretical value of x_{10001} is 244131582646046.
37 The period of this generator is ? FIXME (probably around 2^48). */
39 static inline void rand48_advance (void *vstate);
40 static unsigned long int rand48_get (void *vstate);
41 static double rand48_get_double (void *vstate);
42 static void rand48_set (void *state, unsigned long int s);
44 static const unsigned short int a0 = 0xE66D ;
45 static const unsigned short int a1 = 0xDEEC ;
46 static const unsigned short int a2 = 0x0005 ;
48 static const unsigned short int c0 = 0x000B ;
52 unsigned short int x0, x1, x2;
57 rand48_advance (void *vstate)
59 rand48_state_t *state = (rand48_state_t *) vstate;
61 /* work with unsigned long ints throughout to get correct integer
62 promotions of any unsigned short ints */
64 const unsigned long int x0 = (unsigned long int) state->x0 ;
65 const unsigned long int x1 = (unsigned long int) state->x1 ;
66 const unsigned long int x2 = (unsigned long int) state->x2 ;
71 state->x0 = (a & 0xFFFF) ;
75 /* although the next line may overflow we only need the top 16 bits
76 in the following stage, so it does not matter */
78 a += a0 * x1 + a1 * x0 ;
79 state->x1 = (a & 0xFFFF) ;
82 a += a0 * x2 + a1 * x1 + a2 * x0 ;
83 state->x2 = (a & 0xFFFF) ;
86 static unsigned long int
87 rand48_get (void *vstate)
89 unsigned long int x1, x2;
91 rand48_state_t *state = (rand48_state_t *) vstate;
92 rand48_advance (state) ;
94 x2 = (unsigned long int) state->x2;
95 x1 = (unsigned long int) state->x1;
97 return (x2 << 16) + x1;
101 rand48_get_double (void * vstate)
103 rand48_state_t *state = (rand48_state_t *) vstate;
105 rand48_advance (state) ;
107 return (ldexp((double) state->x2, -16)
108 + ldexp((double) state->x1, -32)
109 + ldexp((double) state->x0, -48)) ;
113 rand48_set (void *vstate, unsigned long int s)
115 rand48_state_t *state = (rand48_state_t *) vstate;
117 if (s == 0) /* default seed */
126 state->x1 = s & 0xFFFF ;
127 state->x2 = (s >> 16) & 0xFFFF ;
133 static const gsl_rng_type rand48_type =
134 {"rand48", /* name */
135 0xffffffffUL, /* RAND_MAX */
137 sizeof (rand48_state_t),
143 const gsl_rng_type *gsl_rng_rand48 = &rand48_type;