1 /* rng/rng.c
2  *
3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
4  *
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.
9  *
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.
14  *
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.
18  */
19 
20 #include "gsl__config.h"
21 #include <stdlib.h>
22 #include <stdio.h>
23 #include <string.h>
24 #include "gsl_errno.h"
25 #include "gsl_rng.h"
26 
27 gsl_rng *
28 gsl_rng_alloc (const gsl_rng_type * T)
29 {
30 
31   gsl_rng *r = (gsl_rng *) malloc (sizeof (gsl_rng));
32 
33   if (r == 0)
34     {
35       GSL_ERROR_VAL ("failed to allocate space for rng struct",
36                         GSL_ENOMEM, 0);
37     };
38 
39   r->state = malloc (T->size);
40 
41   if (r->state == 0)
42     {
43       free (r);         /* exception in constructor, avoid memory leak */
44 
45       GSL_ERROR_VAL ("failed to allocate space for rng state",
46                         GSL_ENOMEM, 0);
47     };
48 
49   r->type = T;
50 
51   gsl_rng_set (r, gsl_rng_default_seed);        /* seed the generator */
52 
53   return r;
54 }
55 
56 int
57 gsl_rng_memcpy (gsl_rng * dest, const gsl_rng * src)
58 {
59   if (dest->type != src->type)
60     {
61       GSL_ERROR ("generators must be of the same type", GSL_EINVAL);
62     }
63 
64   memcpy (dest->state, src->state, src->type->size);
65 
66   return GSL_SUCCESS;
67 }
68 
69 gsl_rng *
70 gsl_rng_clone (const gsl_rng * q)
71 {
72   gsl_rng *r = (gsl_rng *) malloc (sizeof (gsl_rng));
73 
74   if (r == 0)
75     {
76       GSL_ERROR_VAL ("failed to allocate space for rng struct",
77                         GSL_ENOMEM, 0);
78     };
79 
80   r->state = malloc (q->type->size);
81 
82   if (r->state == 0)
83     {
84       free (r);         /* exception in constructor, avoid memory leak */
85 
86       GSL_ERROR_VAL ("failed to allocate space for rng state",
87                         GSL_ENOMEM, 0);
88     };
89 
90   r->type = q->type;
91 
92   memcpy (r->state, q->state, q->type->size);
93 
94   return r;
95 }
96 
97 void
98 gsl_rng_set (const gsl_rng * r, unsigned long int seed)
99 {
100   (r->type->set) (r->state, seed);
101 }
102 
103 #ifndef HIDE_INLINE_STATIC
104 unsigned long int
105 gsl_rng_get (const gsl_rng * r)
106 {
107   return (r->type->get) (r->state);
108 }
109 
110 double
111 gsl_rng_uniform (const gsl_rng * r)
112 {
113   return (r->type->get_double) (r->state);
114 }
115 
116 double
117 gsl_rng_uniform_pos (const gsl_rng * r)
118 {
119   double x ;
120   do
121     {
122       x = (r->type->get_double) (r->state) ;
123     }
124   while (x == 0) ;
125 
126   return x ;
127 }
128 
129 /* Note: to avoid integer overflow in (range+1) we work with scale =
130    range/n = (max-min)/n rather than scale=(max-min+1)/n, this reduces
131    efficiency slightly but avoids having to check for the out of range
132    value.  Note that range is typically O(2^32) so the addition of 1
133    is negligible in most usage. */
134 
135 unsigned long int
136 gsl_rng_uniform_int (const gsl_rng * r, unsigned long int n)
137 {
138   unsigned long int offset = r->type->min;
139   unsigned long int range = r->type->max - offset;
140   unsigned long int scale;
141   unsigned long int k;
142 
143   if (n > range || n == 0)
144     {
145       GSL_ERROR_VAL ("invalid n, either 0 or exceeds maximum value of generator",
146                      GSL_EINVAL, 0) ;
147     }
148 
149   scale = range / n;
150 
151   do
152     {
153       k = (((r->type->get) (r->state)) - offset) / scale;
154     }
155   while (k >= n);
156 
157   return k;
158 }
159 #endif
160 
161 unsigned long int
162 gsl_rng_max (const gsl_rng * r)
163 {
164   return r->type->max;
165 }
166 
167 unsigned long int
168 gsl_rng_min (const gsl_rng * r)
169 {
170   return r->type->min;
171 }
172 
173 const char *
174 gsl_rng_name (const gsl_rng * r)
175 {
176   return r->type->name;
177 }
178 
179 size_t
180 gsl_rng_size (const gsl_rng * r)
181 {
182   return r->type->size;
183 }
184 
185 void *
186 gsl_rng_state (const gsl_rng * r)
187 {
188   return r->state;
189 }
190 
191 void
192 gsl_rng_print_state (const gsl_rng * r)
193 {
194   size_t i;
195   unsigned char *p = (unsigned char *) (r->state);
196   const size_t n = r->type->size;
197 
198   for (i = 0; i < n; i++)
199     {
200       /* FIXME: we're assuming that a char is 8 bits */
201       printf ("%.2x", *(p + i));
202     }
203 
204 }
205 
206 void
207 gsl_rng_free (gsl_rng * r)
208 {
209   free (r->state);
210   free (r);
211 }
212