| /* 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; |