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