diff --git a/src/rng/milcrng.nim b/src/rng/milcrng.nim
index 0dcb2c53..7bcf635d 100644
--- a/src/rng/milcrng.nim
+++ b/src/rng/milcrng.nim
@@ -148,14 +148,15 @@ proc gaussian*(prn: var RngMilc6): float32 =
       result = prn.gset
   else:
     const
-      TINY = 9.999999999999999e-308
+      MAX = 0x01000000.float
+      SCALE1 = 1.0 / (MAX + 1.0)
     var
       v: cdouble
       p: cdouble
       r: cdouble
-    v = prn.uniform
+    v = SCALE1 * (prn.uniform * MAX + 1.0)
     p = prn.uniform * 2.0 * PI
-    r = sqrt(-2.0 * ln(v + TINY))
+    r = sqrt(-2.0 * ln(v))
     result = r * cos(p)
 
 # Only needed for non-vectorized RNGs.
diff --git a/src/rng/mrg32k3a.nim b/src/rng/mrg32k3a.nim
index 2f572e67..dde41c8f 100644
--- a/src/rng/mrg32k3a.nim
+++ b/src/rng/mrg32k3a.nim
@@ -120,6 +120,7 @@ proc seed*(prn: var MRG32k3a; sed,index: auto) {.inline.} =
   seedIndep(prn, ss, index)
 
 proc uniform*(prn:var MRG32k3a):float =
+  ## Uniform random numbers from 0 to 1, excluding 0 and 1.
   var p1,p2:float
   p1 = a12 * prn.s1[1].float - a13n * prn.s1[0].float
   p1 = p1 mod m1
@@ -145,11 +146,11 @@ proc uniform*(prn:var MRG32k3a):float =
 proc gaussian*(prn: var MRG32k3a): float =
   ## Gaussian normal deviate
   ## Probability distribution exp( -x*x/2 ), so < x^2 > = 1
-  const TINY = 9.999999999999999e-308
+  # uniform excludes 0 and 1
   var v,p,r: float
   v = prn.uniform
   p = prn.uniform * 2.0 * PI
-  r = sqrt(-2.0 * ln(v + TINY))
+  r = sqrt(-2.0 * ln(v))
   result = r * cos(p)
 
 import maths/types