htswanalysis/MACS/lib/gsl/gsl-1.11/rng/uni32.c

   1 /* rng/uni32.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 /**
  21   This is a lagged Fibonacci generator which supposedly excellent
  22   statistical properties (I do not concur)
  23
  24   I got it from the net and translated into C.
  25
  26 * ======================================================================
  27 * NIST Guide to Available Math Software.
  28 * Fullsource for module UNI from package CMLIB.
  29 * Retrieved from CAMSUN on Tue Oct  8 14:04:10 1996.
  30 * ======================================================================
  31
  32 C***BEGIN PROLOGUE  UNI
  33 C***DATE WRITTEN   810915
  34 C***REVISION DATE  830805
  35 C***CATEGORY NO.  L6A21
  36 C***KEYWORDS  RANDOM NUMBERS, UNIFORM RANDOM NUMBERS
  37 C***AUTHOR    BLUE, JAMES, SCIENTIFIC COMPUTING DIVISION, NBS
  38 C             KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS
  39 C             MARSAGLIA, GEORGE, COMPUTER SCIENCE DEPT., WASH STATE UNIV
  40 C
  41 C***PURPOSE  THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON [0,1
  42 C             AND CAN BE USED ON ANY COMPUTER WITH WHICH ALLOWS INTEGERS
  43 C             AT LEAST AS LARGE AS 32767.
  44 C***DESCRIPTION
  45 C
  46 C       THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON THE INTER
  47 C       [0,1).  IT CAN BE USED WITH ANY COMPUTER WHICH ALLOWS
  48 C       INTEGERS AT LEAST AS LARGE AS 32767.
  49 C
  50 C
  51 C   USE
  52 C       FIRST TIME....
  53 C                   Z = UNI(JD)
  54 C                     HERE JD IS ANY  N O N - Z E R O  INTEGER.
  55 C                     THIS CAUSES INITIALIZATION OF THE PROGRAM
  56 C                     AND THE FIRST RANDOM NUMBER TO BE RETURNED AS Z.
  57 C       SUBSEQUENT TIMES...
  58 C                   Z = UNI(0)
  59 C                     CAUSES THE NEXT RANDOM NUMBER TO BE RETURNED AS Z.
  60 C
  61 C
  62 C..................................................................
  63 C   NOTE: USERS WHO WISH TO TRANSPORT THIS PROGRAM FROM ONE COMPUTER
  64 C         TO ANOTHER SHOULD READ THE FOLLOWING INFORMATION.....
  65 C
  66 C   MACHINE DEPENDENCIES...
  67 C      MDIG = A LOWER BOUND ON THE NUMBER OF BINARY DIGITS AVAILABLE
  68 C              FOR REPRESENTING INTEGERS, INCLUDING THE SIGN BIT.
  69 C              THIS VALUE MUST BE AT LEAST 16, BUT MAY BE INCREASED
  70 C              IN LINE WITH REMARK A BELOW.
  71 C
  72 C   REMARKS...
  73 C     A. THIS PROGRAM CAN BE USED IN TWO WAYS:
  74 C        (1) TO OBTAIN REPEATABLE RESULTS ON DIFFERENT COMPUTERS,
  75 C            SET 'MDIG' TO THE SMALLEST OF ITS VALUES ON EACH, OR,
  76 C        (2) TO ALLOW THE LONGEST SEQUENCE OF RANDOM NUMBERS TO BE
  77 C            GENERATED WITHOUT CYCLING (REPEATING) SET 'MDIG' TO THE
  78 C            LARGEST POSSIBLE VALUE.
  79 C     B. THE SEQUENCE OF NUMBERS GENERATED DEPENDS ON THE INITIAL
  80 C          INPUT 'JD' AS WELL AS THE VALUE OF 'MDIG'.
  81 C          IF MDIG=16 ONE SHOULD FIND THAT
  82    Editors Note: set the seed using 152 in order to get uni(305)
  83    -jt
  84 C            THE FIRST EVALUATION
  85 C              Z=UNI(305) GIVES Z=.027832881...
  86 C            THE SECOND EVALUATION
  87 C              Z=UNI(0) GIVES   Z=.56102176...
  88 C            THE THIRD EVALUATION
  89 C              Z=UNI(0) GIVES   Z=.41456343...
  90 C            THE THOUSANDTH EVALUATION
  91 C              Z=UNI(0) GIVES   Z=.19797357...
  92 C
  93 C***REFERENCES  MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM
  94 C                 NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U.
  95 C***ROUTINES CALLED  I1MACH,XERROR
  96 C***END PROLOGUE  UNI
  97
  98   **/
  99
 100 #include <config.h>
 101 #include <stdlib.h>
 102 #include <gsl/gsl_rng.h>
 103
 104 static inline unsigned long int uni32_get (void *vstate);
 105 static double uni32_get_double (void *vstate);
 106 static void uni32_set (void *state, unsigned long int s);
 107
 108 static const unsigned long int MDIG = 32;       /* Machine digits in int */
 109 static const unsigned long int m1 = 2147483647;         /* 2^(MDIG-1) - 1 */
 110 static const unsigned long int m2 = 65536;      /* 2^(MDIG/2) */
 111
 112 typedef struct
 113   {
 114     int i, j;
 115     unsigned long m[17];
 116   }
 117 uni32_state_t;
 118
 119 static inline unsigned long
 120 uni32_get (void *vstate)
 121 {
 122   uni32_state_t *state = (uni32_state_t *) vstate;
 123   const long int i = state->i;
 124   const long int j = state->j;
 125
 126   /* important k not be unsigned */
 127   long int k = state->m[i] - state->m[j];
 128
 129   if (k < 0)
 130     k += m1;
 131   state->m[j] = k;
 132
 133   if (i == 0)
 134     {
 135       state->i = 16;
 136     }
 137   else
 138     {
 139       (state->i)--;
 140     }
 141
 142   if (j == 0)
 143     {
 144       state->j = 16;
 145     }
 146   else
 147     {
 148       (state->j)--;
 149     }
 150
 151   return k;
 152 }
 153
 154 static double
 155 uni32_get_double (void *vstate)
 156 {
 157   return uni32_get (vstate) / 2147483647.0 ;
 158 }
 159
 160 static void
 161 uni32_set (void *vstate, unsigned long int s)
 162 {
 163   long int seed, k0, k1, j0, j1;
 164   int i;
 165
 166   uni32_state_t *state = (uni32_state_t *) vstate;
 167
 168   /* For this routine, the seeding is very elaborate! */
 169   /* A flaw in this approach is that seeds 1,2 give exactly the
 170      same random number sequence!  */
 171
 172   /*s = 2*s+1; *//* enforce seed be odd */
 173   seed = (s < m1 ? s : m1);     /* seed should be less than m1 */
 174   seed -= (seed % 2 == 0 ? 1 : 0);
 175
 176   k0 = 9069 % m2;
 177   k1 = 9069 / m2;
 178   j0 = seed % m2;
 179   j1 = seed / m2;
 180
 181   for (i = 0; i < 17; i++)
 182     {
 183       seed = j0 * k0;
 184       j1 = (seed / m2 + j0 * k1 + j1 * k0) % (m2 / 2);
 185       j0 = seed % m2;
 186       state->m[i] = j0 + m2 * j1;
 187     }
 188   state->i = 4;
 189   state->j = 16;
 190
 191   return;
 192 }
 193
 194 static const gsl_rng_type uni32_type =
 195 {"uni32",                       /* name */
 196  2147483646,                    /* RAND_MAX */
 197  0,                             /* RAND_MIN */
 198  sizeof (uni32_state_t),
 199  &uni32_set,
 200  &uni32_get,
 201  &uni32_get_double};
 202
 203 const gsl_rng_type *gsl_rng_uni32 = &uni32_type;