From e37110b08b8a22a05be78784126862c99b3b7baa Mon Sep 17 00:00:00 2001
From: Wirawan Purwanto <wirawan0@gmail.com>
Date: Thu, 21 May 2015 11:30:19 -0400
Subject: [PATCH] * Added more support for legacy 'rannyi' RNG from our Fortran
 codes,   including one from normal gaussian PDF.

---
 math/random/__init__.py  |  2 +-
 math/random/rng_lcg48.py | 42 +++++++++++++++++++++++++++++++++++++++-
 2 files changed, 42 insertions(+), 2 deletions(-)

diff --git a/math/random/__init__.py b/math/random/__init__.py
index d98b90a..f85bfcb 100644
--- a/math/random/__init__.py
+++ b/math/random/__init__.py
@@ -7,7 +7,7 @@
 import numpy
 
 class rng_base(object):
-  """Bas class for random number generator."""
+  """Base class for random number generator."""
 
 # Standard classes
 from wpylib.math.random.rng_lcg48 import lcg48
diff --git a/math/random/rng_lcg48.py b/math/random/rng_lcg48.py
index b2de4de..eff8834 100644
--- a/math/random/rng_lcg48.py
+++ b/math/random/rng_lcg48.py
@@ -10,7 +10,8 @@ Implementation of 48-bit Linear Congruential pseudorandom number generator
 
 This RNG is used only for testing and regression purposes.
 
-
+This module contains a python implementation of the Fortran 'rannyu' function
+(including the hardcoded seed library).
 """
 
 from wpylib.math.random import rng_base
@@ -43,12 +44,51 @@ class lcg48(rng_base):
     self.m = 34522712143931  # 11**13
     self.update_seed(rannyu_seed_lib[index])
 
+
+class rannyu(lcg48):
+  """Simulation of Fortran rannyu method with support for normal
+  (Gaussian) PDF."""
+  def __init__(self, seed=None):
+    if seed==None:
+      seed = numpy.random.randint(102)+1
+    self.use_seed_lib(seed)
+  def normal(self, size):
+    rslt = numpy.empty(size)
+    size = numpy.product(rslt.shape)
+    for i in xrange(size):
+      rslt.flat[i] = GaussianRandom(self)
+    return rslt
+
+
+def GaussianRandom(rnd):
+  """Same implementation as in my perl or Fortran code.
+  rnd is a random number generator object (i.e. the rannyu object)."""
+  from numpy import sqrt, log
+  if not hasattr(rnd, "_gr_next_rnd"):
+    while True:
+      x = 2.0 * rnd() - 1.0
+      y = 2.0 * rnd() - 1.0
+      r2 = x*x + y*y
+      if 0.0 < r2 < 1.0: break
+
+    z = sqrt(-2.0 * log(r2) / r2)
+    rnd._gr_next_rnd = z * x  # v
+    return z * y   # u
+  else:
+    v = rnd._gr_next_rnd
+    del rnd._gr_next_rnd
+    return v
+
+
 def mkint48(N):
   """Generates a 48-bit number from four integers.
   Used for encoding conversion for legacy rannyu library."""
   return N[3] + (N[2] << 12) + (N[1] << 24) + (N[0] << 36)
 
+
 def split48(value):
+  """Splits a 48-bit integer into a 4-tuple of 12-bit integers.
+  """
   i4 = value & 0xFFF
   value >>= 12
   i3 = value & 0xFFF