examples/smpi/NAS/common/randdp.c

   1 //---------------------------------------------------------------------
   2 //   This function is C verson of random number generator randdp.f
   3 //---------------------------------------------------------------------
   4
   5 double  randlc(X, A)
   6 double *X;
   7 double *A;
   8 {
   9       static int        KS=0;
  10       static double  R23, R46, T23, T46;
  11       double    T1, T2, T3, T4;
  12       double    A1;
  13       double    A2;
  14       double    X1;
  15       double    X2;
  16       double    Z;
  17       int         i, j;
  18
  19       if (KS == 0)
  20       {
  21         R23 = 1.0;
  22         R46 = 1.0;
  23         T23 = 1.0;
  24         T46 = 1.0;
  25
  26         for (i=1; i<=23; i++)
  27         {
  28           R23 = 0.50 * R23;
  29           T23 = 2.0 * T23;
  30         }
  31         for (i=1; i<=46; i++)
  32         {
  33           R46 = 0.50 * R46;
  34           T46 = 2.0 * T46;
  35         }
  36         KS = 1;
  37       }
  38
  39 /*  Break A into two parts such that A = 2^23 * A1 + A2 and set X = N.  */
  40
  41       T1 = R23 * *A;
  42       j  = T1;
  43       A1 = j;
  44       A2 = *A - T23 * A1;
  45
  46 /*  Break X into two parts such that X = 2^23 * X1 + X2, compute
  47     Z = A1 * X2 + A2 * X1  (mod 2^23), and then
  48     X = 2^23 * Z + A2 * X2  (mod 2^46).                            */
  49
  50       T1 = R23 * *X;
  51       j  = T1;
  52       X1 = j;
  53       X2 = *X - T23 * X1;
  54       T1 = A1 * X2 + A2 * X1;
  55
  56       j  = R23 * T1;
  57       T2 = j;
  58       Z = T1 - T23 * T2;
  59       T3 = T23 * Z + A2 * X2;
  60       j  = R46 * T3;
  61       T4 = j;
  62       *X = T3 - T46 * T4;
  63       return(R46 * *X);
  64 }