Summary: Random from Brent in Algol W

diff --git a/src/random.alw b/src/random.alw
new file mode 100644 (file)
index 0000000..052a308
--- /dev/null
+++ b/src/random.alw
@@ -0,0 +1,38 @@
+COMMENT:       RANDOM NUMBER GENERATOR
+                ***********************
+
+  PROCEDURE RANDOM RETURNS A LONG REAL RANDOM NUMBER UNIFORMLY
+      DISTRIBUTED IN (0,1) (INCLUDING 0 BUT NOT 1),
+      RANINIT(R) WITH R ANY INTEGER MUST BE CALLED FOR
+      INITIALIZATION BEFORE THE FIRST CALL TO RANDOM, AND THE
+      DECLARATIONS OF RAN1, RAN2 AND RAN3 MUST BE GLOBAL,
+      THE ALGORITHM RETURNS X(N)/2**56, WHERE
+                   X(N) = X(N-1) + X(N-127) (MOO 2**56),
+      SINCE 1 + X + X**127 IS PRIMITIVE (MOD 2), THE PERIOD IS AT
+      LEAST 2**127 - 1 > 10**38, SEE KNUTH (1969), PP. 26, 34, 464,
+      X(N) IS STORED IN A LONG REAL WORD AS
+      RAN3 = X(N)/2**56 - 1/2, AND ALL FLOATING POINT ARITHMETIC
+      IS EXACT;
+
+  LONG REAL PROCEDURE RANDOM(INTEGER VALUE NAUGHT);
+  BEGIN
+    LONG REAL RAN1; INTEGER RAN2; LONG REAL ARRAY RAN3 (0::126);
+    INTEGER R; LOGICAL INIT;
+    INIT := FALSE;
+    IF INIT THEN GO TO L3;
+    R := ABS(NAUGHT) REM 8190 + 1;
+    RAN2 := 127; WHILE RAN2 > 0 DO
+    BEGIN RAN2 := RAN2 - 1; RAN1 := -2L**55;
+      FOR I := 1 UNTIL 7 DO
+      BEGIN R := (1756*R) REM 8191;
+        RAN1 := (RAN1 + (R DIV 32) )*( 1/256) ;
+      END;
+      RAN3 (RAN2) := RAN1
+    END;
+    INIT := TRUE;
+  L3: RAN2 := IF RAN2 = 0 THEN 126 ELSE RAN2 - 1;
+    RAN1 := RAN1 + RAN3 (RAN2);
+    RAN3 (RAN2) := RAN1 := IF RAN1 < 0L THEN RAN1 + 0.5L
+                        ELSE RAN1 - 0.5L;
+    RAN1 + 0.5L
+  END RANDOM.