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

   1 /* rng/slatec.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
  22 * ======================================================================
  23 * NIST Guide to Available Math Software.
  24 * Source for module RAND from package CMLIB.
  25 * Retrieved from TIBER on Fri Oct 11 11:43:42 1996.
  26 * ======================================================================
  27       FUNCTION RAND(R)
  28 C***BEGIN PROLOGUE  RAND
  29 C***DATE WRITTEN   770401   (YYMMDD)
  30 C***REVISION DATE  820801   (YYMMDD)
  31 C***CATEGORY NO.  L6A21
  32 C***KEYWORDS  RANDOM NUMBER,SPECIAL FUNCTION,UNIFORM
  33 C***AUTHOR  FULLERTON, W., (LANL)
  34 C***PURPOSE  Generates a uniformly distributed random number.
  35 C***DESCRIPTION
  36 C
  37 C      This pseudo-random number generator is portable among a wide
  38 C variety of computers.  RAND(R) undoubtedly is not as good as many
  39 C readily available installation dependent versions, and so this
  40 C routine is not recommended for widespread usage.  Its redeeming
  41 C feature is that the exact same random numbers (to within final round-
  42 C off error) can be generated from machine to machine.  Thus, programs
  43 C that make use of random numbers can be easily transported to and
  44 C checked in a new environment.
  45 C      The random numbers are generated by the linear congruential
  46 C method described, e.g., by Knuth in Seminumerical Methods (p.9),
  47 C Addison-Wesley, 1969.  Given the I-th number of a pseudo-random
  48 C sequence, the I+1 -st number is generated from
  49 C             X(I+1) = (A*X(I) + C) MOD M,
  50 C where here M = 2**22 = 4194304, C = 1731 and several suitable values
  51 C of the multiplier A are discussed below.  Both the multiplier A and
  52 C random number X are represented in double precision as two 11-bit
  53 C words.  The constants are chosen so that the period is the maximum
  54 C possible, 4194304.
  55 C      In order that the same numbers be generated from machine to
  56 C machine, it is necessary that 23-bit integers be reducible modulo
  57 C 2**11 exactly, that 23-bit integers be added exactly, and that 11-bit
  58 C integers be multiplied exactly.  Furthermore, if the restart option
  59 C is used (where R is between 0 and 1), then the product R*2**22 =
  60 C R*4194304 must be correct to the nearest integer.
  61 C      The first four random numbers should be .0004127026,
  62 C .6750836372, .1614754200, and .9086198807.  The tenth random number
  63 C is .5527787209, and the hundredth is .3600893021 .  The thousandth
  64 C number should be .2176990509 .
  65 C      In order to generate several effectively independent sequences
  66 C with the same generator, it is necessary to know the random number
  67 C for several widely spaced calls.  The I-th random number times 2**22,
  68 C where I=K*P/8 and P is the period of the sequence (P = 2**22), is
  69 C still of the form L*P/8.  In particular we find the I-th random
  70 C number multiplied by 2**22 is given by
  71 C I   =  0  1*P/8  2*P/8  3*P/8  4*P/8  5*P/8  6*P/8  7*P/8  8*P/8
  72 C RAND=  0  5*P/8  2*P/8  7*P/8  4*P/8  1*P/8  6*P/8  3*P/8  0
  73 C Thus the 4*P/8 = 2097152 random number is 2097152/2**22.
  74 C      Several multipliers have been subjected to the spectral test
  75 C (see Knuth, p. 82).  Four suitable multipliers roughly in order of
  76 C goodness according to the spectral test are
  77 C    3146757 = 1536*2048 + 1029 = 2**21 + 2**20 + 2**10 + 5
  78 C    2098181 = 1024*2048 + 1029 = 2**21 + 2**10 + 5
  79 C    3146245 = 1536*2048 +  517 = 2**21 + 2**20 + 2**9 + 5
  80 C    2776669 = 1355*2048 + 1629 = 5**9 + 7**7 + 1
  81 C
  82 C      In the table below LOG10(NU(I)) gives roughly the number of
  83 C random decimal digits in the random numbers considered I at a time.
  84 C C is the primary measure of goodness.  In both cases bigger is better.
  85 C
  86 C                   LOG10 NU(I)              C(I)
  87 C       A       I=2  I=3  I=4  I=5    I=2  I=3  I=4  I=5
  88 C
  89 C    3146757    3.3  2.0  1.6  1.3    3.1  1.3  4.6  2.6
  90 C    2098181    3.3  2.0  1.6  1.2    3.2  1.3  4.6  1.7
  91 C    3146245    3.3  2.2  1.5  1.1    3.2  4.2  1.1  0.4
  92 C    2776669    3.3  2.1  1.6  1.3    2.5  2.0  1.9  2.6
  93 C   Best
  94 C    Possible   3.3  2.3  1.7  1.4    3.6  5.9  9.7  14.9
  95 C
  96 C             Input Argument --
  97 C R      If R=0., the next random number of the sequence is generated.
  98 C        If R .LT. 0., the last generated number will be returned for
  99 C          possible use in a restart procedure.
 100 C        If R .GT. 0., the sequence of random numbers will start with
 101 C          the seed R mod 1.  This seed is also returned as the value of
 102 C          RAND provided the arithmetic is done exactly.
 103 C
 104 C             Output Value --
 105 C RAND   a pseudo-random number between 0. and 1.
 106 C***REFERENCES  (NONE)
 107 C***ROUTINES CALLED  (NONE)
 108 C***END PROLOGUE  RAND
 109       DATA IA1, IA0, IA1MA0 /1536, 1029, 507/
 110       DATA IC /1731/
 111       DATA IX1, IX0 /0, 0/
 112 C***FIRST EXECUTABLE STATEMENT  RAND
 113       IF (R.LT.0.) GO TO 10
 114       IF (R.GT.0.) GO TO 20
 115 C
 116 C           A*X = 2**22*IA1*IX1 + 2**11*(IA1*IX1 + (IA1-IA0)*(IX0-IX1)
 117 C                   + IA0*IX0) + IA0*IX0
 118 C
 119       IY0 = IA0*IX0
 120       IY1 = IA1*IX1 + IA1MA0*(IX0-IX1) + IY0
 121       IY0 = IY0 + IC
 122       IX0 = MOD (IY0, 2048)
 123       IY1 = IY1 + (IY0-IX0)/2048
 124       IX1 = MOD (IY1, 2048)
 125 C
 126  10   RAND = IX1*2048 + IX0
 127       RAND = RAND / 4194304.
 128       RETURN
 129 C
 130  20   IX1 = AMOD(R,1.)*4194304. + 0.5
 131       IX0 = MOD (IX1, 2048)
 132       IX1 = (IX1-IX0)/2048
 133       GO TO 10
 134 C
 135       END
 136
 137   **/
 138
 139 #include <config.h>
 140 #include <stdlib.h>
 141 #include <gsl/gsl_rng.h>
 142
 143 static inline unsigned long int slatec_get (void *vstate);
 144 static double slatec_get_double (void *vstate);
 145 static void slatec_set (void *state, unsigned long int s);
 146
 147 typedef struct
 148   {
 149     long int x0, x1;
 150   }
 151 slatec_state_t;
 152
 153 static const long P = 4194304;
 154 static const long a1 = 1536;
 155 static const long a0 = 1029;
 156 static const long a1ma0 = 507;
 157 static const long c = 1731;
 158
 159 static inline unsigned long int
 160 slatec_get (void *vstate)
 161 {
 162   long y0, y1;
 163   slatec_state_t *state = (slatec_state_t *) vstate;
 164
 165   y0 = a0 * state->x0;
 166   y1 = a1 * state->x1 + a1ma0 * (state->x0 - state->x1) + y0;
 167   y0 = y0 + c;
 168   state->x0 = y0 % 2048;
 169   y1 = y1 + (y0 - state->x0) / 2048;
 170   state->x1 = y1 % 2048;
 171
 172   return state->x1 * 2048 + state->x0;
 173 }
 174
 175 static double
 176 slatec_get_double (void *vstate)
 177 {
 178   return slatec_get (vstate) / 4194304.0 ;
 179 }
 180
 181 static void
 182 slatec_set (void *vstate, unsigned long int s)
 183 {
 184   slatec_state_t *state = (slatec_state_t *) vstate;
 185
 186   /* Only eight seeds are permitted.  This is pretty limiting, but
 187      at least we are guaranteed that the eight sequences are different */
 188
 189   s = s % 8;
 190   s *= P / 8;
 191
 192   state->x0 = s % 2048;
 193   state->x1 = (s - state->x0) / 2048;
 194 }
 195
 196 static const gsl_rng_type slatec_type =
 197 {"slatec",                      /* name */
 198  4194303,                       /* RAND_MAX */
 199  0,                             /* RAND_MIN */
 200  sizeof (slatec_state_t),
 201  &slatec_set,
 202  &slatec_get,
 203  &slatec_get_double};
 204
 205 const gsl_rng_type *gsl_rng_slatec = &slatec_type;