Re: KISS 64-bit pseudo-random number generator

On 9/12/24 21:10, Krishna Myneni wrote:
...

 I ran a simple Monte-Carlo integration of the area in a unit circle using the 64-bit LCG PRNG (from random.4th in kForth dist) and Marsaglia's 64-bit KISS PRNG. Hard to say which is better for this problem.
 -- Krishna
 === Comparison of 64-bit LCG and KISS PRNGs ===
 PRNG: RANDOM  (random.4th)
Ntrials   Area      log(error)
-------------------------------
10^2     3.28        -0.86
10^3     3.18        -1.42
10^4     3.1724      -1.51
10^5     3.14048     -2.95
10^6     3.14350     -2.72
10^7     3.14225     -3.18
10^8     3.14152     -4.14
 PRNG: RAN-KISS (kiss.4th)
Ntrials   Area      log(error)
-------------------------------
10^2     3.24       -1.01
10^3     3.10       -1.38
10^4     3.1252     -1.79
10^5     3.14084    -3.12
10^6     3.14073    -3.06
10^7     3.14184    -3.61
10^8     3.14167    -4.11
 

A bit more involved test of the same 64-bit PRNGs starts to show the possibility of defects in Marsaglia's 64-bit kiss prng (RAN-KISS) compared to a simple 64-bit LCG prng (RANDOM). Instead of drawing from a uniform distribution, as assumed for these generators, the tests below draw from a Maxwell-Boltzmann distribution, which describes the probability distribution for speeds of ideal gas atoms at a given temperature (non-interacting atoms except for hard sphere collisions in three dimensions). The different "moments" of the speed may be computed from the samples and compared with exact theoretical results: <v>, <v^2>, <v^3>, etc., where v is a random sample of the speed for the M-B distribution at given temperature and mass of atom. The parameter "a", given by
a = m/(2*kB*T)
where m is the mass, T is the temperature, and kB is the Boltzmann constant, determines the shape of the M-B distribution. Unsigned 64-bit random numbers, generated by the PRNG under test, are converted to floating point numbers in the interval [0, 1]. Ideally, the PRNG will distribute these samples uniformly over the interval. The samples are converted to random speed samples with an M-B distribution by an inverse mapping of the cumulative probability distribution for M-B to the speed. The moments for <v>, <v^2>, and <v^3> are computed from these samples.
In this simulation, the temperature is set at 300 K, and the mass of the He atom is used. Theory predictions for this m and T are
<v>_th = sqrt( 4/(pi*a) ) = 1258.727 m/s
<v^2>_th = 3/(2*a) = 1869538 (m/s)^2
<v^3>_th = sqrt( 16/(pi*a^3) ) = 3140144 x 10^3 (m/s)^3
Below are the results for the two PRNGs. For N trials of 10^5 through 10^7, the simple 64-bit LCG prng consistently gives less relative error than Marsaglia's 64-bit kiss prng with respect to the theory predictions.
Source code for these test results is given below.
--
Krishna
=== begin tests ===
Test of prng RANDOM 64-bit (random.4th):
' random test-prng
Moments of speed
  N       <v> (m/s)    <v^2> (m/s)^2    <v^3> (m/s)^3
10^2     1220.2444     1761665.1       2879481740.
10^3     1258.2315     1844889.5       3025976380.
10^4     1250.0696     1837249.9       3055346856.
10^5     1259.9899     1870367.4       3143506292.
10^6     1259.3923     1868383.4       3136761299.
10^7     1259.6343     1869366.9       3140010276.
Relative Errors
  N       |e1|       |e2|       |e3|
10^2  3.13e-02   5.77e-02   8.30e-02
10^3  1.19e-03   1.32e-02   3.64e-02
10^4  7.67e-03   1.73e-02   2.70e-02
10^5  2.08e-04   4.44e-04   1.07e-03
10^6  2.66e-04   6.18e-04   1.08e-03
10^7  7.39e-05   9.15e-05   4.27e-05
Test of prng RAN-KISS 64-bit (kiss.4th):
' ran-kiss test-prng
Moments of speed
  N       <v> (m/s)    <v^2> (m/s)^2    <v^3> (m/s)^3
10^2     1212.8900     1702106.0       2666594142.
10^3     1274.7097     1900519.4       3190481667.
10^4     1259.6892     1866276.3       3123781859.
10^5     1260.4578     1872500.5       3147907366.
10^6     1260.2589     1871249.6       3145073404.
10^7     1259.8296     1869882.9       3140954422.
Relative Errors
  N       |e1|       |e2|       |e3|
10^2  3.72e-02   8.96e-02   1.51e-01
10^3  1.19e-02   1.66e-02   1.60e-02
10^4  3.03e-05   1.74e-03   5.21e-03
10^5  5.80e-04   1.58e-03   2.47e-03
10^6  4.22e-04   9.16e-04   1.57e-03
10^7  8.11e-05   1.85e-04   2.58e-04
=== end tests ===
=== begin test-prng-mb.4th ===
\ test-prng-mb.4th
\
\ Test a pseudo-random number generator by computing the
\ moments of powers of speed using random samples from a
\ Maxwell-Boltzmann distribution of ideal gas speeds in 3D.
\
\ Copyright 2024 Krishna Myneni
\
\ Permission is granted to use this code for any purpose,
\ provided the original source is referenced.
\
\ Usage Example:
\
\   ' random TEST-PRNG
\
\ where RANDOM is a user-provided PRNG which returns
\ a cell-sized pseudo-random number.
\
\ Notes
\
\ 1. The xt for any PRNG may be passed to the word
\    TEST-PRNG provided it returns a cell-sized
\    bit pattern. All bits of the cell are assumed
\    to be part of the pseudo-random number.
\
\ 2. The default M-B distribution is computed for
\    helium atoms at 300K. If you change the default
\    temperature or atomic mass, the maximum velocity
\    used in the root finder call in CDF^-1 may need
\    to be changed.
\
include ans-words
include strings
include phyconsts
include fsl/fsl-util
include fsl/erf
include fsl/regfalsi
include random
include kiss
[UNDEFINED] FPICK [IF]
cr cr .( This program requires a separate FP stack! ) cr
quit
[THEN]
4.0e0 pi f* fconstant 4pi
-1 0 d>f fconstant F_MAX_U
300.0e0  fconstant DEF_TEMPERATURE  \ K
1.0e4    fconstant DEF_VMAX         \ m/s
1.0e-9   fconstant DEF_TOL   \ root-finder tolerance
4.002602e0  fconstant m_He   \ mass in amu for He atom
: fcube ( F: r -- r^3 )  fdup fsquare f* ;
: f3dup 2 fpick 2 fpick 2 fpick ;
[UNDEFINED] FTUCK [IF]
: ftuck ( F: a b -- b a b ) fswap fover ;
[THEN]
\ Relative error between a value and its reference value
: rel-error ( F: a aref -- e ) ftuck f- fswap f/ ;
\ The random number generator to be tested must return a
\ cell-sized bit pattern (unsigned single) -- all bits of
\ the cell are assumed to be randomized.
defer prngU
defer prngU-init
' random is prngU   \ default PRNG
:noname ( -- )  1000 0 do prngU drop loop ;  is prngU-init
\ Generate a floating point random number between 0 and 1
\ using the selected unsigned single cell prng.
: randf ( F: -- r ) prngU 0 d>f  F_MAX_U f/ ;
fvariable m_atom_kg
: set-atomic-mass ( F: m_amu -- )  kg/amu f* m_atom_kg f! ;
m_He set-atomic-mass  \ default value of mass of atom
fvariable a
: set-T ( F: Tk -- )
     kB f* 2.0e0 f* m_atom_kg f@ fswap f/ a f! ;
\ Return probability density, p(v), for the M-B distribution
\ with parameter "a"
: pdf ( F: v -- p )
     fsquare fdup a f@ f* fnegate fexp f*
     4pi f* a f@ pi f/ 1.5e0 f** f* ;
\ Return cumulative probability, C(v1) i.e. P(0 <= v < v1),
\ the integral of p(v) from 0 to v1 for the M-B distribution.
: cdf ( F: v1 -- P )
     fdup fsquare a f@ f* fdup fsqrt erf
     fswap fnegate fexp frot f* 2e f* a f@ pi f/ fsqrt f* f- ;
\ Invert cdf to obtain velocity from value of cdf
fvariable pv1
: cdf^-1 ( F: P -- v )
    pv1 f!
    [: cdf pv1 f@ f- ;] 0.0e0  DEF_VMAX  DEF_TOL )root ;
\ N random trials of velocities following the M-B distribution
\ are sampled and the moments v, v^2, and v^3 are computed
\ from the samples and compared to theory.
fvariable vsum
fvariable v2sum
fvariable v3sum
: mb-prng ( Ntrials -- ) ( F: -- <v> <v^2> <v^3> )
     DEF_TEMPERATURE set-T
     0.0e0 fdup fdup vsum f! v2sum f! v3sum f!
     dup
     0 ?DO
       randf   \ random draw of cdf value
       cdf^-1  \ invert cdf to get random v in M-B distribution
       fdup         vsum f+!
       fdup fsquare v2sum f+!
            fcube   v3sum f+!
     LOOP
     s>f
     vsum  f@ fover f/ fswap
     v2sum f@ fover f/ fswap
     v3sum f@ fswap f/ ;
\ Theoretical moments of v for the M-B distribution at given
\ temperature and mass, specified by parameter "a".
: <v>_th   ( F: -- <v> )   4.0e0 pi    a f@ f* f/ fsqrt ;
: <v^2>_th ( F: -- <v^2> ) 3.0e0 2.0e0 a f@ f* f/ ;
: <v^3>_th ( F: -- <v^3> ) 16.0e0 pi   a f@ fcube f* f/ fsqrt ;
\ Compute relative errors of randomly sampled averages with
\ respect to the theoretical moments.
: moment-errors ( F: <v> <v^2> <v^3> -- e1 e2 e3 )
     frot <v>_th   rel-error
     frot <v^2>_th rel-error
     frot <v^3>_th rel-error ;
9 constant MAX-POW
MAX-POW 1+  3 float matrix v_mom{{
MAX-POW 1+  3 float matrix m_rel{{
: print-tables ( -- )
     cr ." Moments of speed"
     cr ."  N       <v> (m/s)    <v^2> (m/s)^2    <v^3> (m/s)^3" cr
     MAX-POW 1- 2 DO
       v_mom{{ I 2 }} f@ v_mom{{ I 1 }} f@ v_mom{{ I 0 }} f@
       ." 10^" I 1 .r 2 spaces
       12 4 f.rd 2 spaces
       12 1 f.rd 6 spaces
       12 0 f.rd cr
     LOOP
     cr ." Relative Errors"
     cr ."  N       |e1|       |e2|       |e3|"  cr
     precision  \ save precision
     3 set-precision
     MAX-POW 1- 2 DO
       m_rel{{ I 2 }} f@ m_rel{{ I 1 }} f@ m_rel{{ I 0 }} f@
       ." 10^" I 1 .r 2 spaces
       fabs fs. 2 spaces
       fabs fs. 2 spaces
       fabs fs. cr
     LOOP
     set-precision \ restore precision
     cr
;
: test-prng ( xt -- )
     is prngU
     prngU-init
     MAX-POW 1- 2 DO
       I s>f falog fround>s
       mb-prng
       f3dup
       v_mom{{ I 2 }} f!  v_mom{{ I 1 }} f! v_mom{{ I 0 }} f!
       moment-errors
       m_rel{{ I 2 }} f!  m_rel{{ I 1 }} f!  m_rel{{ I 0 }} f!
     LOOP
     print-tables
;
cr cr .( Usage: ' <prng> test-prng )
cr cr .(   e.g. ' random test-prng ) cr cr
=== end test-prng-mb.4th ===