examples/smpi/NAS/common/randdp.c

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