Xorshift128+ : quick yet high quality 64-bit pseudo random number generation

xorshift128plus.fsx

(* 
   The Xorshift128+ generator generates a high-quality 64-bit value very quickly (see reference). 
   It is presently used in the JavaScript engines of Chrome, Firefox and Safari, and it is the default generator in Erlang. 
   If you find its period too short for large-scale parallel simulations, use xorshift1024*.

   Note that the lowest bit of this generator is an LSFR, 
   and thus it is slightly less random than the other bits. 
   We suggest to use a sign test to extract a random Boolean value.

   Reference : http://xoroshiro.di.unimi.it/
               http://xoroshiro.di.unimi.it/xorshift128plus.c
*)

namespace Random

open System

/// The Xorshift128+ generator generates a high-quality 64-bit value very quickly
/// The state must be seeded so that it is not everywhere zero. 
type Xorshift128plus(seed0,seed1) =

   let mutable seed0 = seed0
   let mutable seed1 = seed1

   /// The Xorshift128+ generator generates a high-quality 64-bit value very quickly
   new() =
      // If you have a 64-bit seed, we suggest to seed a splitmix64 generator.
      let stdRng = Random()
      let seed0 = stdRng.Next() |> uint64
      let seed1 = stdRng.Next() |> uint64
      // The state must be seeded so that it is not everywhere zero. 
      if seed0 = 0UL && seed1 = 0UL then Xorshift128plus() else 
         Xorshift128plus(seed0,seed1)

   //---------

   /// Produce a pseudo-random uint64 between 0 and UInt64.MaxValue, update the generator
   member this.NextUInt64() =
      let result = seed1 + seed0
      let s1 = seed0 ^^^ (seed0 <<< 23) // a
      seed0 <- seed1
      seed1 <- s1 ^^^ seed1 ^^^ (s1 >>> 18) ^^^ (seed1 >>> 5) // b, c
      result

   /// Produce a pseudo-random int32 between 0 and Int.MaxValue, update the generator
   member this.NextInt() =
      let k = this.NextUInt64()
      int (k / 8589934596UL) // UInt64.MaxValue / Int32.MaxValue

   /// Produce a pseudo-random int32 between 0 and borneSup, update the generator
   member this.NextInt(borneSup) =
      let k = this.NextUInt64()
      k * ( (uint64 borneSup) / 8589934596UL) |> int // UInt64.MaxValue / Int32.MaxValue

   /// Produce a pseudo-random float between 0 and 1, update the generator
   member this.NextFloat() =
      let k = this.NextUInt64()
      (float k) / 18446744073709551615. // UInt64.MaxValue

   /// Produce a pseudo-random float between 0 and borneSup, update the generator
   member this.NextFloat(borneSup) =
      let k = this.NextUInt64()
      (float k) * (borneSup / 18446744073709551615.) // UInt64.MaxValue

   /// Produce a pseudo-random bool, update the generator
   member this.NextBool() =
      let k = this.NextUInt64()
      k > 9223372036854775807UL // UInt64.MaxValue / 2UL

   /// Produce a pseudo-random bool with a given probability, update the generator
   member this.NextBool(trueProba) =
      this.NextFloat() < trueProba

   //----------

   /// Equivalent to 2^64 call to next
   /// Can be used to generate 2^64 non-overlapping subsequences for parallel computations.
   member this.Jump() =
      let jump = [|9970080660168977919UL; 1305993406145048470UL|] // 0x8a5cd789635d2dff, 0x121fd2155c472f96
      let mutable s0 = 0UL
      let mutable s1 = 0UL
      for i = 0 to 1 do
         for b = 0 to 63 do
            if (jump.[i] &&& 1UL <<< b) = 1UL then
               s0 <- s0 ^^^ seed0
               s1 <- s1 ^^^ seed1
            this.NextUInt64() |> ignore
      seed0 <- s0
      seed1 <- s1

   /// Create a new generator identical to the current generator
   member this.Copy() = Xorshift128plus(seed0,seed1)

   /// Create n non overlapping generators
   /// Can be used to generate up to 2^64 non-overlapping subsequences for parallel computations.
   static member CreateNonOverlapping n =
      let gen = Xorshift128plus()
      Array.init n (fun _ -> gen.Jump() ; gen.Copy() )