SplitMix implementation

splitmix.u

-- porting https://hackage.haskell.org/package/splitmix-0.1.0.1/docs/src/System.Random.SplitMix.html
-- could use a native popcount / hamming weight instruction

roll : Nat -> Nat ->{Random} Nat
roll n sides =
  Nat.sum (List.replicate n '(Random.natIn sides + 1))

-- example, rolling 3 x d8
> Random.splitmix 10349 '(roll 3 8)

> Random.splitmix 9383938 '(List.replicate 10 '(Random.natIn 50))

ability Random where
  nextNat : Nat
  split : Boolean -> () 

Nat.MaxValue = 1 `shiftLeft` 63
Float.Ulp = 1.0 / Nat.toFloat (1 `Nat.shiftLeft` 53)

Random.nat : '{Random} Nat
Random.nat _ = Random.nextNat

Random.natIn : Nat ->{Random} Nat
Random.natIn range = 
  if range == 0 then bug (Random.natIn, range)
  else Random.natTo (drop range 1)

Random.natTo : Nat ->{Random} Nat
Random.natTo range =
  mask = complement 0 `shiftRight` leadingZeros range
  go _ = 
    x = !Random.nat
    x' = x `Nat.and` mask
    if x' > range then !go
    else x'
  !go

Random.float : '{Random} Float
Random.float _ = Nat.toFloat (!Random.nat `shiftRight` 11) * Float.Ulp   

Random.splitmix : Nat -> '{g, Random} a ->{g} a
Random.splitmix seed c = handle !c with Random.splitmix.handler1 seed

Random.splitmix.handler1 seed = 
  Random.splitmix.handler (mix64 seed) (mixGamma (seed + goldenGamma))   

Random.splitmix.handler : Nat ->{} Nat ->{} Request {Random} a -> a
Random.splitmix.handler seed gamma = cases
  { a } -> a
  { Random.split isLeft -> resume } -> 
    seed' = seed + gamma
    seed'' = seed' + gamma
    if isLeft then 
      handle !resume with Random.splitmix.handler seed'' gamma
    else 
      handle !resume with Random.splitmix.handler (mix64 seed') (mixGamma seed'') 
  { Random.nextNat -> resume } -> 
    seed' = seed + gamma
    handle resume (mix64 seed')
    with Random.splitmix.handler seed' gamma

Random.splitmix.internals.goldenGamma = 11400714819323198485
Random.splitmix.internals.shiftXor n w = w `Nat.xor` (w `Nat.shiftRight` n)
Random.splitmix.internals.shiftXorMultiply n k w = shiftXor n w * k
Random.splitmix.internals.mix64 z0 =
   -- MurmurHash3Mixer
   z1 = shiftXorMultiply 33 18397679294719823053 z0
   z2 = shiftXorMultiply 33 14181476777654086739 z1
   shiftXor 33 z2

-- used only in mixGamma
Random.splitmix.internals.mix64variant13 z0 =
   -- Better Bit Mixing - Improving on MurmurHash3's 64-bit Finalizer
   -- http://zimbry.blogspot.fi/2011/09/better-bit-mixing-improving-on.html
   --
   -- Stafford's Mix13
   z1 = shiftXorMultiply 30 13787848793156543929 z0 -- MurmurHash3 mix constants
   z2 = shiftXorMultiply 27 10723151780598845931 z1
   shiftXor 31 z2

Random.splitmix.internals.mixGamma z0 =
  z1 = mix64variant13 z0 `Nat.or` 1             -- force to be odd
  n  = hammingWeight (z1 `xor` (z1 `shiftRight` 1))
  -- see: http://www.pcg-random.org/posts/bugs-in-splitmix.html
  -- let's trust the text of the paper, not the code.
  if n >= 24 then z1
  else z1 `xor` 12297829382473034410