website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/rand.f

   1 *DECK RAND
   2       FUNCTION RAND (R)
   3 C***BEGIN PROLOGUE  RAND
   4 C***PURPOSE  Generate a uniformly distributed random number.
   5 C***LIBRARY   SLATEC (FNLIB)
   6 C***CATEGORY  L6A21
   7 C***TYPE      SINGLE PRECISION (RAND-S)
   8 C***KEYWORDS  FNLIB, RANDOM NUMBER, SPECIAL FUNCTIONS, UNIFORM
   9 C***AUTHOR  Fullerton, W., (LANL)
  10 C***DESCRIPTION
  11 C
  12 C      This pseudo-random number generator is portable among a wide
  13 C variety of computers.  RAND(R) undoubtedly is not as good as many
  14 C readily available installation dependent versions, and so this
  15 C routine is not recommended for widespread usage.  Its redeeming
  16 C feature is that the exact same random numbers (to within final round-
  17 C off error) can be generated from machine to machine.  Thus, programs
  18 C that make use of random numbers can be easily transported to and
  19 C checked in a new environment.
  20 C
  21 C      The random numbers are generated by the linear congruential
  22 C method described, e.g., by Knuth in Seminumerical Methods (p.9),
  23 C Addison-Wesley, 1969.  Given the I-th number of a pseudo-random
  24 C sequence, the I+1 -st number is generated from
  25 C             X(I+1) = (A*X(I) + C) MOD M,
  26 C where here M = 2**22 = 4194304, C = 1731 and several suitable values
  27 C of the multiplier A are discussed below.  Both the multiplier A and
  28 C random number X are represented in double precision as two 11-bit
  29 C words.  The constants are chosen so that the period is the maximum
  30 C possible, 4194304.
  31 C
  32 C      In order that the same numbers be generated from machine to
  33 C machine, it is necessary that 23-bit integers be reducible modulo
  34 C 2**11 exactly, that 23-bit integers be added exactly, and that 11-bit
  35 C integers be multiplied exactly.  Furthermore, if the restart option
  36 C is used (where R is between 0 and 1), then the product R*2**22 =
  37 C R*4194304 must be correct to the nearest integer.
  38 C
  39 C      The first four random numbers should be .0004127026,
  40 C .6750836372, .1614754200, and .9086198807.  The tenth random number
  41 C is .5527787209, and the hundredth is .3600893021 .  The thousandth
  42 C number should be .2176990509 .
  43 C
  44 C      In order to generate several effectively independent sequences
  45 C with the same generator, it is necessary to know the random number
  46 C for several widely spaced calls.  The I-th random number times 2**22,
  47 C where I=K*P/8 and P is the period of the sequence (P = 2**22), is
  48 C still of the form L*P/8.  In particular we find the I-th random
  49 C number multiplied by 2**22 is given by
  50 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
  51 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
  52 C Thus the 4*P/8 = 2097152 random number is 2097152/2**22.
  53 C
  54 C      Several multipliers have been subjected to the spectral test
  55 C (see Knuth, p. 82).  Four suitable multipliers roughly in order of
  56 C goodness according to the spectral test are
  57 C    3146757 = 1536*2048 + 1029 = 2**21 + 2**20 + 2**10 + 5
  58 C    2098181 = 1024*2048 + 1029 = 2**21 + 2**10 + 5
  59 C    3146245 = 1536*2048 +  517 = 2**21 + 2**20 + 2**9 + 5
  60 C    2776669 = 1355*2048 + 1629 = 5**9 + 7**7 + 1
  61 C
  62 C      In the table below LOG10(NU(I)) gives roughly the number of
  63 C random decimal digits in the random numbers considered I at a time.
  64 C C is the primary measure of goodness.  In both cases bigger is better.
  65 C
  66 C                   LOG10 NU(I)              C(I)
  67 C       A       I=2  I=3  I=4  I=5    I=2  I=3  I=4  I=5
  68 C
  69 C    3146757    3.3  2.0  1.6  1.3    3.1  1.3  4.6  2.6
  70 C    2098181    3.3  2.0  1.6  1.2    3.2  1.3  4.6  1.7
  71 C    3146245    3.3  2.2  1.5  1.1    3.2  4.2  1.1  0.4
  72 C    2776669    3.3  2.1  1.6  1.3    2.5  2.0  1.9  2.6
  73 C   Best
  74 C    Possible   3.3  2.3  1.7  1.4    3.6  5.9  9.7  14.9
  75 C
  76 C             Input Argument --
  77 C R      If R=0., the next random number of the sequence is generated.
  78 C        If R .LT. 0., the last generated number will be returned for
  79 C          possible use in a restart procedure.
  80 C        If R .GT. 0., the sequence of random numbers will start with
  81 C          the seed R mod 1.  This seed is also returned as the value of
  82 C          RAND provided the arithmetic is done exactly.
  83 C
  84 C             Output Value --
  85 C RAND   a pseudo-random number between 0. and 1.
  86 C
  87 C***REFERENCES  (NONE)
  88 C***ROUTINES CALLED  (NONE)
  89 C***REVISION HISTORY  (YYMMDD)
  90 C   770401  DATE WRITTEN
  91 C   890531  Changed all specific intrinsics to generic.  (WRB)
  92 C   890531  REVISION DATE from Version 3.2
  93 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  94 C***END PROLOGUE  RAND
  95       SAVE IA1, IA0, IA1MA0, IC, IX1, IX0
  96       DATA IA1, IA0, IA1MA0 /1536, 1029, 507/
  97       DATA IC /1731/
  98       DATA IX1, IX0 /0, 0/
  99 C***FIRST EXECUTABLE STATEMENT  RAND
 100       IF (R.LT.0.) GO TO 10
 101       IF (R.GT.0.) GO TO 20
 102 C
 103 C           A*X = 2**22*IA1*IX1 + 2**11*(IA1*IX1 + (IA1-IA0)*(IX0-IX1)
 104 C                   + IA0*IX0) + IA0*IX0
 105 C
 106       IY0 = IA0*IX0
 107       IY1 = IA1*IX1 + IA1MA0*(IX0-IX1) + IY0
 108       IY0 = IY0 + IC
 109       IX0 = MOD (IY0, 2048)
 110       IY1 = IY1 + (IY0-IX0)/2048
 111       IX1 = MOD (IY1, 2048)
 112 C
 113  10   RAND = IX1*2048 + IX0
 114       RAND = RAND / 4194304.
 115       RETURN
 116 C
 117  20   IX1 = MOD(R,1.)*4194304. + 0.5
 118       IX0 = MOD (IX1, 2048)
 119       IX1 = (IX1-IX0)/2048
 120       GO TO 10
 121 C
 122       END