-//---------------------------------------------------------------------
-// This function is C verson of random number generator randdp.f
-//---------------------------------------------------------------------
-
-double randlc(X, A)
-double *X;
-double *A;
+/*
+ * FUNCTION RANDLC (X, A)
+ *
+ * This routine returns a uniform pseudorandom double precision number in the
+ * range (0, 1) by using the linear congruential generator
+ *
+ * x_{k+1} = a x_k (mod 2^46)
+ *
+ * where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers
+ * before repeating. The argument A is the same as 'a' in the above formula,
+ * and X is the same as x_0. A and X must be odd double precision integers
+ * in the range (1, 2^46). The returned value RANDLC is normalized to be
+ * between 0 and 1, i.e. RANDLC = 2^(-46) * x_1. X is updated to contain
+ * the new seed x_1, so that subsequent calls to RANDLC using the same
+ * arguments will generate a continuous sequence.
+ *
+ * This routine should produce the same results on any computer with at least
+ * 48 mantissa bits in double precision floating point data. On Cray systems,
+ * double precision should be disabled.
+ *
+ * David H. Bailey October 26, 1990
+ *
+ * IMPLICIT DOUBLE PRECISION (A-H, O-Z)
+ * SAVE KS, R23, R46, T23, T46
+ * DATA KS/0/
+ *
+ * If this is the first call to RANDLC, compute R23 = 2 ^ -23, R46 = 2 ^ -46,
+ * T23 = 2 ^ 23, and T46 = 2 ^ 46. These are computed in loops, rather than
+ * by merely using the ** operator, in order to insure that the results are
+ * exact on all systems. This code assumes that 0.5D0 is represented exactly.
+ */
+double randlc(double *X, double*A)
{
- static int KS=0;
- static double R23, R46, T23, T46;
- double T1, T2, T3, T4;
- double A1;
- double A2;
- double X1;
- double X2;
- double Z;
- int i, j;
+ static int KS=0;
+ static double R23, R46, T23, T46;
+ double T1, T2, T3, T4;
+ double A1, A2;
+ double X1, X2;
+ double Z;
+ int i, j;
- if (KS == 0)
- {
- R23 = 1.0;
- R46 = 1.0;
- T23 = 1.0;
- T46 = 1.0;
-
- for (i=1; i<=23; i++)
- {
- R23 = 0.50 * R23;
- T23 = 2.0 * T23;
- }
- for (i=1; i<=46; i++)
- {
- R46 = 0.50 * R46;
- T46 = 2.0 * T46;
- }
- KS = 1;
- }
+ if (KS == 0) {
+ R23 = 1.0;
+ R46 = 1.0;
+ T23 = 1.0;
+ T46 = 1.0;
-/* Break A into two parts such that A = 2^23 * A1 + A2 and set X = N. */
+ for (i=1; i<=23; i++) {
+ R23 = 0.50 * R23;
+ T23 = 2.0 * T23;
+ }
+ for (i=1; i<=46; i++) {
+ R46 = 0.50 * R46;
+ T46 = 2.0 * T46;
+ }
+ KS = 1;
+ }
- T1 = R23 * *A;
- j = T1;
- A1 = j;
- A2 = *A - T23 * A1;
+/* Break A into two parts such that A = 2^23 * A1 + A2 and set X = N. */
+ T1 = R23 * *A;
+ j = T1;
+ A1 = j;
+ A2 = *A - T23 * A1;
/* Break X into two parts such that X = 2^23 * X1 + X2, compute
- Z = A1 * X2 + A2 * X1 (mod 2^23), and then
- X = 2^23 * Z + A2 * X2 (mod 2^46). */
+ Z = A1 * X2 + A2 * X1 (mod 2^23), and then X = 2^23 * Z + A2 * X2 (mod 2^46). */
+ T1 = R23 * *X;
+ j = T1;
+ X1 = j;
+ X2 = *X - T23 * X1;
+ T1 = A1 * X2 + A2 * X1;
- T1 = R23 * *X;
- j = T1;
- X1 = j;
- X2 = *X - T23 * X1;
- T1 = A1 * X2 + A2 * X1;
-
- j = R23 * T1;
- T2 = j;
- Z = T1 - T23 * T2;
- T3 = T23 * Z + A2 * X2;
- j = R46 * T3;
- T4 = j;
- *X = T3 - T46 * T4;
- return(R46 * *X);
-}
+ j = R23 * T1;
+ T2 = j;
+ Z = T1 - T23 * T2;
+ T3 = T23 * Z + A2 * X2;
+ j = R46 * T3;
+ T4 = j;
+ *X = T3 - T46 * T4;
+ return(R46 * *X);
+}
