Modified MILC RNG to avoid extreme values in Gaussian according to "r…

…ng: avoid extreme values in gaussian jcosborn#8" commit in main branch

src/rng.nim

@@ -1,7 +1,7 @@
-import rng/milcrng_spv # Modified for Pauli-Villars
+import rng/milcrng
 import rng/mrg32k3a
 import rng/distribution
 
-export milcrng_spv
+export milcrng
 export mrg32k3a
 export distribution

src/rng/milcrng.nim

@@ -99,6 +99,8 @@ else:
     prn.multiplier = 100005'u32 + 8'u32 * index
     #prn.addend = 12345
     #prn.scale = 1.0 / float32(0x01000000)
+    #prn.iset = 1 # ~~~~~~~~~
+    #prn.gset = 0 # ~~~~~~~~~
 proc seedIndep*(prn: var RngMilc6; sed,index: auto) {.inline.} =
   seedX(prn, sed.uint32, index.uint32)
 proc seed*(prn: var RngMilc6; sed,index: auto) {.inline.} =
@@ -148,14 +150,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.

src/rng/mrg32k3a.nim

@@ -134,6 +134,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
@@ -160,10 +161,11 @@ 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