diff options
Diffstat (limited to 'gsl-1.9/rng/uni32.c')
| -rw-r--r-- | gsl-1.9/rng/uni32.c | 203 |
1 files changed, 203 insertions, 0 deletions
diff --git a/gsl-1.9/rng/uni32.c b/gsl-1.9/rng/uni32.c new file mode 100644 index 0000000..ba5536a --- /dev/null +++ b/gsl-1.9/rng/uni32.c @@ -0,0 +1,203 @@ +/* rng/uni32.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/** + This is a lagged Fibonacci generator which supposedly excellent + statistical properties (I do not concur) + + I got it from the net and translated into C. + +* ====================================================================== +* NIST Guide to Available Math Software. +* Fullsource for module UNI from package CMLIB. +* Retrieved from CAMSUN on Tue Oct 8 14:04:10 1996. +* ====================================================================== + +C***BEGIN PROLOGUE UNI +C***DATE WRITTEN 810915 +C***REVISION DATE 830805 +C***CATEGORY NO. L6A21 +C***KEYWORDS RANDOM NUMBERS, UNIFORM RANDOM NUMBERS +C***AUTHOR BLUE, JAMES, SCIENTIFIC COMPUTING DIVISION, NBS +C KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS +C MARSAGLIA, GEORGE, COMPUTER SCIENCE DEPT., WASH STATE UNIV +C +C***PURPOSE THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON [0,1 +C AND CAN BE USED ON ANY COMPUTER WITH WHICH ALLOWS INTEGERS +C AT LEAST AS LARGE AS 32767. +C***DESCRIPTION +C +C THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON THE INTER +C [0,1). IT CAN BE USED WITH ANY COMPUTER WHICH ALLOWS +C INTEGERS AT LEAST AS LARGE AS 32767. +C +C +C USE +C FIRST TIME.... +C Z = UNI(JD) +C HERE JD IS ANY N O N - Z E R O INTEGER. +C THIS CAUSES INITIALIZATION OF THE PROGRAM +C AND THE FIRST RANDOM NUMBER TO BE RETURNED AS Z. +C SUBSEQUENT TIMES... +C Z = UNI(0) +C CAUSES THE NEXT RANDOM NUMBER TO BE RETURNED AS Z. +C +C +C.................................................................. +C NOTE: USERS WHO WISH TO TRANSPORT THIS PROGRAM FROM ONE COMPUTER +C TO ANOTHER SHOULD READ THE FOLLOWING INFORMATION..... +C +C MACHINE DEPENDENCIES... +C MDIG = A LOWER BOUND ON THE NUMBER OF BINARY DIGITS AVAILABLE +C FOR REPRESENTING INTEGERS, INCLUDING THE SIGN BIT. +C THIS VALUE MUST BE AT LEAST 16, BUT MAY BE INCREASED +C IN LINE WITH REMARK A BELOW. +C +C REMARKS... +C A. THIS PROGRAM CAN BE USED IN TWO WAYS: +C (1) TO OBTAIN REPEATABLE RESULTS ON DIFFERENT COMPUTERS, +C SET 'MDIG' TO THE SMALLEST OF ITS VALUES ON EACH, OR, +C (2) TO ALLOW THE LONGEST SEQUENCE OF RANDOM NUMBERS TO BE +C GENERATED WITHOUT CYCLING (REPEATING) SET 'MDIG' TO THE +C LARGEST POSSIBLE VALUE. +C B. THE SEQUENCE OF NUMBERS GENERATED DEPENDS ON THE INITIAL +C INPUT 'JD' AS WELL AS THE VALUE OF 'MDIG'. +C IF MDIG=16 ONE SHOULD FIND THAT + Editors Note: set the seed using 152 in order to get uni(305) + -jt +C THE FIRST EVALUATION +C Z=UNI(305) GIVES Z=.027832881... +C THE SECOND EVALUATION +C Z=UNI(0) GIVES Z=.56102176... +C THE THIRD EVALUATION +C Z=UNI(0) GIVES Z=.41456343... +C THE THOUSANDTH EVALUATION +C Z=UNI(0) GIVES Z=.19797357... +C +C***REFERENCES MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM +C NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U. +C***ROUTINES CALLED I1MACH,XERROR +C***END PROLOGUE UNI + + **/ + +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_rng.h> + +static inline unsigned long int uni32_get (void *vstate); +static double uni32_get_double (void *vstate); +static void uni32_set (void *state, unsigned long int s); + +static const unsigned long int MDIG = 32; /* Machine digits in int */ +static const unsigned long int m1 = 2147483647; /* 2^(MDIG-1) - 1 */ +static const unsigned long int m2 = 65536; /* 2^(MDIG/2) */ + +typedef struct + { + int i, j; + unsigned long m[17]; + } +uni32_state_t; + +static inline unsigned long +uni32_get (void *vstate) +{ + uni32_state_t *state = (uni32_state_t *) vstate; + const long int i = state->i; + const long int j = state->j; + + /* important k not be unsigned */ + long int k = state->m[i] - state->m[j]; + + if (k < 0) + k += m1; + state->m[j] = k; + + if (i == 0) + { + state->i = 16; + } + else + { + (state->i)--; + } + + if (j == 0) + { + state->j = 16; + } + else + { + (state->j)--; + } + + return k; +} + +static double +uni32_get_double (void *vstate) +{ + return uni32_get (vstate) / 2147483647.0 ; +} + +static void +uni32_set (void *vstate, unsigned long int s) +{ + long int seed, k0, k1, j0, j1; + int i; + + uni32_state_t *state = (uni32_state_t *) vstate; + + /* For this routine, the seeding is very elaborate! */ + /* A flaw in this approach is that seeds 1,2 give exactly the + same random number sequence! */ + + /*s = 2*s+1; *//* enforce seed be odd */ + seed = (s < m1 ? s : m1); /* seed should be less than m1 */ + seed -= (seed % 2 == 0 ? 1 : 0); + + k0 = 9069 % m2; + k1 = 9069 / m2; + j0 = seed % m2; + j1 = seed / m2; + + for (i = 0; i < 17; i++) + { + seed = j0 * k0; + j1 = (seed / m2 + j0 * k1 + j1 * k0) % (m2 / 2); + j0 = seed % m2; + state->m[i] = j0 + m2 * j1; + } + state->i = 4; + state->j = 16; + + return; +} + +static const gsl_rng_type uni32_type = +{"uni32", /* name */ + 2147483646, /* RAND_MAX */ + 0, /* RAND_MIN */ + sizeof (uni32_state_t), + &uni32_set, + &uni32_get, + &uni32_get_double}; + +const gsl_rng_type *gsl_rng_uni32 = &uni32_type; |