examples/smpi/NAS/EP-sampling/randlc.c

   1
   2 /*
   3  *    FUNCTION RANDLC (X, A)
   4  *
   5  *  This routine returns a uniform pseudorandom double precision number in the
   6  *  range (0, 1) by using the linear congruential generator
   7  *
   8  *  x_{k+1} = a x_k  (mod 2^46)
   9  *
  10  *  where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
  11  *  before repeating.  The argument A is the same as 'a' in the above formula,
  12  *  and X is the same as x_0.  A and X must be odd double precision integers
  13  *  in the range (1, 2^46).  The returned value RANDLC is normalized to be
  14  *  between 0 and 1, i.e. RANDLC = 2^(-46) * x_1.  X is updated to contain
  15  *  the new seed x_1, so that subsequent calls to RANDLC using the same
  16  *  arguments will generate a continuous sequence.
  17  *
  18  *  This routine should produce the same results on any computer with at least
  19  *  48 mantissa bits in double precision floating point data.  On Cray systems,
  20  *  double precision should be disabled.
  21  *
  22  *  David H. Bailey     October 26, 1990
  23  *
  24  *     IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  25  *     SAVE KS, R23, R46, T23, T46
  26  *     DATA KS/0/
  27  *
  28  *  If this is the first call to RANDLC, compute R23 = 2 ^ -23, R46 = 2 ^ -46,
  29  *  T23 = 2 ^ 23, and T46 = 2 ^ 46.  These are computed in loops, rather than
  30  *  by merely using the ** operator, in order to insure that the results are
  31  *  exact on all systems.  This code assumes that 0.5D0 is represented exactly.
  32  */
  33
  34
  35 /*****************************************************************/
  36 /*************           R  A  N  D  L  C             ************/
  37 /*************                                        ************/
  38 /*************    portable random number generator    ************/
  39 /*****************************************************************/
  40
  41 double  randlc( double *X, double *A )
  42 {
  43       static int        KS=0;
  44       static double  R23, R46, T23, T46;
  45       double    T1, T2, T3, T4;
  46       double    A1;
  47       double    A2;
  48       double    X1;
  49       double    X2;
  50       double    Z;
  51       int         i, j;
  52
  53       if (KS == 0)
  54       {
  55         R23 = 1.0;
  56         R46 = 1.0;
  57         T23 = 1.0;
  58         T46 = 1.0;
  59
  60         for (i=1; i<=23; i++)
  61         {
  62           R23 = 0.50 * R23;
  63           T23 = 2.0 * T23;
  64         }
  65         for (i=1; i<=46; i++)
  66         {
  67           R46 = 0.50 * R46;
  68           T46 = 2.0 * T46;
  69         }
  70         KS = 1;
  71       }
  72
  73 /*  Break A into two parts such that A = 2^23 * A1 + A2 and set X = N.  */
  74
  75       T1 = R23 * *A;
  76       j  = T1;
  77       A1 = j;
  78       A2 = *A - T23 * A1;
  79
  80 /*  Break X into two parts such that X = 2^23 * X1 + X2, compute
  81     Z = A1 * X2 + A2 * X1  (mod 2^23), and then
  82     X = 2^23 * Z + A2 * X2  (mod 2^46).                            */
  83
  84       T1 = R23 * *X;
  85       j  = T1;
  86       X1 = j;
  87       X2 = *X - T23 * X1;
  88       T1 = A1 * X2 + A2 * X1;
  89
  90       j  = R23 * T1;
  91       T2 = j;
  92       Z = T1 - T23 * T2;
  93       T3 = T23 * Z + A2 * X2;
  94       j  = R46 * T3;
  95       T4 = j;
  96       *X = T3 - T46 * T4;
  97       return(R46 * *X);
  98 }
  99
 100
 101
 102 /*****************************************************************/
 103 /************   F  I  N  D  _  M  Y  _  S  E  E  D    ************/
 104 /************                                         ************/
 105 /************ returns parallel random number seq seed ************/
 106 /*****************************************************************/
 107