diff --git a/lib/bitgen.ml b/lib/bitgen.ml index 60fdc8e..37ead76 100644 --- a/lib/bitgen.ml +++ b/lib/bitgen.ml @@ -30,8 +30,9 @@ ]} *) +module SeedSequence = Seed.SeedSequence + module SFC64 = Sfc.SFC64 module PCG64 = Pcg.PCG64 -module SeedSequence = Seed.SeedSequence module Xoshiro256 = Xoshiro.Xoshiro256StarStar module Philox64 = Philox.Philox diff --git a/lib/common.ml b/lib/common.ml new file mode 100644 index 0000000..3eda26d --- /dev/null +++ b/lib/common.ml @@ -0,0 +1,51 @@ +open Stdint + + +type 'a store = + | Empty + | Value of 'a + + +let max_int = Int.max_int |> Uint32.of_int + +(* [to_uint32 nxt s store] advances state [s] and optains a random uint64 integers + using function [nxt]. The resulting int is split into 2 uint32 integers where + one of them is returned to the caller while the next is put into [store]. If + [store] is not empty then the [nxt] is not called and the integer stored in + [store] is returned and an empty store is returned to the caller. A 3-tuple is + returned of the form: (new uint32, new state s, new store). *) +let next_uint32 ~next s = function + | Value x -> x, s, Empty + | Empty -> + let uint, s' = next s in + Uint64.(logand uint max_int |> to_uint32), s', + Value Uint64.(shift_right uint 32 |> to_uint32) + + +(* [next_double nextu64 t] returns a random float from bitgenerator [t] using + function [nextu64] and a new bitgenerator with its internal state advanced forward by a step. *) +let next_double ~nextu64 t = match nextu64 t with + | u, t' -> Uint64.(shift_right u 11 |> to_int) |> Float.of_int + |> ( *. ) (1.0 /. 9007199254740992.0), t' + + +module type BITGEN = sig + type t + (** [t] is the state of the bitgenerator. *) + + val next_uint64 : t -> uint64 * t + (** [next_uint64 t] Generates a random unsigned 64-bit integer and a state + of the generator advanced forward by one step. *) + + val next_uint32 : t -> uint32 * t + (** [next_uint32 t] Generates a random unsigned 32-bit integer and a state + of the generator advanced forward by one step. *) + + val next_double : t -> float * t + (** [next_double t] Generates a random 64 bit float and a state of the + generator advanced forward by one step. *) + + val initialize : Seed.SeedSequence.t -> t + (** [initialize s] Returns the initial state of the generator. The random stream + is determined by the initialization of the seed sequence [s] of {!SeedSequence.t} type. *) +end diff --git a/lib/pcg.ml b/lib/pcg.ml index 7a1d12a..6154fb3 100644 --- a/lib/pcg.ml +++ b/lib/pcg.ml @@ -18,23 +18,7 @@ module PCG64 : sig The input seed is processed by {!SeedSequence} to generate both values. *) - type t - (** [t] is the state of the PCG64 bitgenerator *) - - val next_uint64 : t -> uint64 * t - (** Generate a random unsigned 64-bit integer and return a state of the - generator advanced by one step forward *) - - val next_uint32 : t -> uint32 * t - (** Generate a random unsigned 32-bit integer and return a state of the - generator advanced by one step forward *) - - val next_double : t -> float * t - (** Generate a random 64 bit float and return a state of the - generator advanced by one step forward *) - - val initialize : Seed.SeedSequence.t -> t - (** Get the initial state of the generator using a {!SeedSequence} type as input *) + include Common.BITGEN val advance : int128 -> t -> t (** [advance delta] Advances the underlying RNG as if [delta] draws have been made. @@ -44,7 +28,7 @@ module PCG64 : sig (** [next_bounded_uint64 bound t] returns an unsigned 64bit integers in the range (0, bound) as well as the state of the generator advanced one step forward. *) end = struct - type t = {s : setseq; has_uint32 : bool; uinteger : uint32} + type t = {s : setseq; ustore : uint32 Common.store} and setseq = {state : uint128; increment : uint128} let multiplier = Uint128.of_string "0x2360ed051fc65da44385df649fccf645" @@ -55,7 +39,7 @@ end = struct and r = Uint128.(shift_right state 122 |> to_int) in let nr = Uint32.(of_int r |> neg |> logand (of_int 63) |> to_int) in Uint64.(logor (shift_left v nr) (shift_right v r)) - + let next {state; increment} = let state' = Uint128.(state * multiplier + increment) in @@ -66,22 +50,15 @@ end = struct | u, s' -> u, {t with s = s'} - let next_uint32 t = match t.has_uint32 with - | true -> t.uinteger, {t with has_uint32 = false} - | false -> - let uint, s' = next t.s in - Uint64.(of_int 0xffffffff |> logand uint |> to_uint32), - {uinteger = Uint64.(shift_right uint 32 |> to_uint32); - has_uint32 = true; - s = s'} + let next_uint32 t = + match Common.next_uint32 ~next:next t.s t.ustore with + | u, s, ustore -> u, {s; ustore} - let next_double t = match next_uint64 t with - | u, t' -> - Uint64.(shift_right u 11 |> to_int) |> Float.of_int |> ( *. ) (1.0 /. 9007199254740992.0), t' + let next_double t = Common.next_double ~nextu64:next_uint64 t - let advance delta {s = {state; increment}; has_uint32; uinteger} = + let advance delta {s = {state; increment}; ustore} = let open Uint128 in let rec lcg d am ap cm cp = (* advance state using LCG method *) match d = zero, logand d one = one with @@ -89,8 +66,7 @@ end = struct | false, true -> lcg (shift_right d 1) (am * cm) (ap * cm + cp) (cm * cm) (cp * (cm + one)) | false, false -> lcg (shift_right d 1) am ap (cm * cm) (cp * (cm + one)) in - {s = {state = lcg (Uint128.of_int128 delta) one zero multiplier increment; increment}; - has_uint32; uinteger} + {s = {state = lcg (Uint128.of_int128 delta) one zero multiplier increment; increment}; ustore} let set_seed (shi, slo, ihi, ilo) = @@ -112,7 +88,5 @@ end = struct let initialize seed = let state = Seed.SeedSequence.generate_64bit_state 4 seed in - {s = set_seed (state.(0), state.(1), state.(2), state.(3)); - uinteger = Uint32.zero; - has_uint32 = false} + {s = set_seed (state.(0), state.(1), state.(2), state.(3)); ustore = Common.Empty} end diff --git a/lib/philox.ml b/lib/philox.ml index 1e3a52c..46186fb 100644 --- a/lib/philox.ml +++ b/lib/philox.ml @@ -34,23 +34,7 @@ module Philox : sig |> List.map Philox64.initialize ]} *) - type t - (** [t] is the state of the Philox64 bitgenerator *) - - val next_uint64 : t -> uint64 * t - (** Generate a random unsigned 64-bit integer and return a state of the - generator advanced by one step forward *) - - val next_uint32 : t -> uint32 * t - (** Generate a random unsigned 32-bit integer and return a state of the - generator advanced by one step forward *) - - val next_double : t -> float * t - (** Generate a random 64 bit float and return a state of the - generator advanced by one step forward *) - - val initialize : Seed.SeedSequence.t -> t - (** Get the initial state of the generator using a {!SeedSequence} type as input *) + include Common.BITGEN val initialize_ctr : counter:uint64 * uint64 * uint64 * uint64 -> Seed.SeedSequence.t -> t (** Get the initial state of the generator using a 4-element unsigned 64-bit tuple as @@ -59,17 +43,15 @@ module Philox : sig val jump : t -> t (** [jump t] is equivalent to {m 2^{128}} calls to {!Philox64.next_uint64}. *) - end = struct type t = { - ctr : counter; key: key; - buffer_pos : int; - buffer : uint64 * uint64 * uint64 * uint64; - has_uint32 : bool; - uinteger : uint32} + ctr : counter; + buffer : buffer; + ustore : uint32 Common.store} and counter = uint64 * uint64 * uint64 * uint64 and key = uint64 * uint64 + and buffer = Buffer of int * counter let bumpk0 = Uint64.of_string "0x9E3779B97F4A7C15" @@ -113,27 +95,20 @@ end = struct | 0 -> b0 | 1 -> b1 | 2 -> b2 | _ -> b3 - let next_uint64 t = match t.buffer_pos with - | i when i < 4 -> index t.buffer i, {t with buffer_pos = i + 1} + let next_uint64 t = match t.buffer with + | Buffer (i, buf) when i < 4 -> index buf i, {t with buffer = Buffer (i + 1, buf)} | _ -> let ctr' = next t.ctr in let buf = ten_rounds ctr' t.key in - index buf 0, {t with ctr = ctr'; buffer = buf; buffer_pos = 1} + index buf 0, {t with ctr = ctr'; buffer = Buffer (1, buf)} let next_uint32 t = - match t.has_uint32 with - | true -> t.uinteger, {t with has_uint32 = false} - | false -> - let uint, t' = next_uint64 t in - Uint64.(of_int 0xffffffff |> logand uint |> to_uint32), (* low 32 bits *) - {t' with has_uint32 = true; - uinteger = Uint64.(shift_right uint 32 |> to_uint32)} (* high 32 bits *) + match Common.next_uint32 ~next:next_uint64 t t.ustore with + | u, s, ustore -> u, {s with ustore = ustore} - let next_double t = match next_uint64 t with - | u, t' -> Uint64.(shift_right u 11 |> to_int) - |> Float.of_int |> ( *. ) (1.0 /. 9007199254740992.0), t' + let next_double t = Common.next_double ~nextu64:next_uint64 t let jump t = @@ -150,11 +125,9 @@ end = struct let initialize_ctr ~counter seed = let istate = Seed.SeedSequence.generate_64bit_state 2 seed in {ctr = counter; - buffer = zeros; - key = (istate.(0), istate.(1)); - uinteger = Uint32.zero; - has_uint32 = false; - buffer_pos = 4} + ustore = Common.Empty; + buffer = Buffer (4, zeros); + key = (istate.(0), istate.(1))} let initialize seed = initialize_ctr ~counter:zeros seed diff --git a/lib/sfc.ml b/lib/sfc.ml index 28554c5..1eb722d 100644 --- a/lib/sfc.ml +++ b/lib/sfc.ml @@ -16,56 +16,28 @@ module SFC64 : sig each iteration. The input seed is processed by {!SeedSequence} to generate the first 3 values, then the algorithm is iterated a small number of times to mix.*) - type t - (** [t] is the state of the SFC64 bitgenerator *) - - val next_uint64 : t -> uint64 * t - (** Generate a random unsigned 64-bit integer and return a state of the - generator advanced by one step forward *) - - val next_uint32 : t -> uint32 * t - (** Generate a random unsigned 32-bit integer and return a state of the - generator advanced by one step forward *) - - val next_double : t -> float * t - (** Generate a random 64 bit float and return a state of the - generator advanced by one step forward *) - - val initialize : Seed.SeedSequence.t -> t - (** Get the initial state of the generator using a {!SeedSequence} type as input *) - + include Common.BITGEN end = struct (* last uint64 value is the counter *) - type t = {s : state; has_uint32 : bool; uinteger : uint32} + type t = {s : state; ustore : uint32 Common.store} and state = uint64 * uint64 * uint64 * uint64 - let next (w, x, y, z) = - let uint = Uint64.(w + x + z) - and y' = Uint64.(logor (shift_left y 24) (shift_right y 40)) in - uint, Uint64.(shift_right x 11 |> logxor x, y + shift_left y 3, y' + uint, z + one) + let next (w, x, y, z) = match Uint64.(w + x + z) with + | u -> u, Uint64.(shift_right x 11 |> logxor x, y + shift_left y 3, + logor (shift_left y 24) (shift_right y 40) + u, z + one) - let next_uint64 t = - let uint, s' = next t.s in - uint, {t with s = s'} + let next_uint64 t = match next t.s with + | u, s' -> u, {t with s = s'} let next_uint32 t = - match t.has_uint32 with - | true -> t.uinteger, {t with has_uint32 = false} - | false -> let uint, s' = next t.s in - Uint64.(of_int 0xffffffff |> logand uint |> to_uint32), { - uinteger = Uint64.(shift_right uint 32 |> to_uint32); - has_uint32 = true; - s = s' - } + match Common.next_uint32 ~next:next t.s t.ustore with + | u, s, ustore -> u, {s; ustore} - let next_double t = - let uint, t' = next_uint64 t in - let rnd = Uint64.(shift_right uint 11 |> to_int) |> Float.of_int in - rnd *. (1.0 /. 9007199254740992.0), t' + let next_double t = Common.next_double ~nextu64:next_uint64 t let set_seed (w, x, y) = @@ -78,6 +50,5 @@ end = struct let initialize seed = let istate = Seed.SeedSequence.generate_64bit_state 3 seed in - {s = set_seed (istate.(0), istate.(1), istate.(2)); - has_uint32 = false; uinteger = Uint32.zero} + {s = set_seed (istate.(0), istate.(1), istate.(2)); ustore = Common.Empty} end diff --git a/lib/xoshiro.ml b/lib/xoshiro.ml index 60a3ba6..7f467e9 100644 --- a/lib/xoshiro.ml +++ b/lib/xoshiro.ml @@ -20,68 +20,37 @@ module Xoshiro256StarStar : sig random numbers have been generated. This allows the original sequence to be split so that distinct segments can be used in each worker process.*) - type t - (** [t] is the state of the Xoshiro256** bitgenerator *) - - val next_uint64 : t -> uint64 * t - (** Generate a random unsigned 64-bit integer and return a state of the - generator advanced by one step forward *) - - val next_uint32 : t -> uint32 * t - (** Generate a random unsigned 32-bit integer and return a state of the - generator advanced by one step forward *) - - val next_double : t -> float * t - (** Generate a random 64 bit float and return a state of the - generator advanced by one step forward *) - - val initialize : Seed.SeedSequence.t -> t - (** Get the initial state of the generator using a {!SeedSequence} type as input *) + include Common.BITGEN val jump : t -> t (** [jump t] is equivalent to {m 2^{128}} calls to {!Xoshiro256.next_uint64}; it can be used to generate {m 2^{128}} non-overlapping subsequences for parallel computations. *) end = struct - - type t = {s : state; has_uint32 : bool; uinteger : uint32} + type t = {s : state; ustore : uint32 Common.store} and state = uint64 * uint64 * uint64 * uint64 let rotl x k = - 64 - k |> Uint64.shift_right x |> Uint64.logor (Uint64.shift_left x k) + Uint64.shift_right x (64 - k) |> Uint64.(logor (shift_left x k)) - let next ((w, x, y, z) : state) : uint64 * state = - let open Uint64 in - let y' = logxor y w in - let z' = logxor z x in - of_int 9 * rotl (x * of_int 5) 7, - (logxor w z', logxor x y', logxor y' (shift_left x 17), rotl z' 45) + let next (w, x, y, z) = + let open Uint64 in match logxor y w, logxor z x with + | y', z' -> of_int 9 * rotl (x * of_int 5) 7, + (logxor w z', logxor x y', logxor y' (shift_left x 17), rotl z' 45) - let next_uint64 t = - let u, s' = next t.s in - u, {t with s = s'} + let next_uint64 t = match next t.s with + | u, s' -> u, {t with s = s'} let next_uint32 t = - match t.has_uint32 with - | true -> t.uinteger, {t with has_uint32 = false} - | false -> - let uint, s' = next t.s in - Uint64.(of_int 0xffffffff |> logand uint |> to_uint32), - {uinteger = Uint64.(shift_right uint 32 |> to_uint32); - has_uint32 = true; - s = s'} + match Common.next_uint32 ~next:next t.s t.ustore with + | u, s, ustore -> u, {s; ustore} - let next_double t = - let uint, t' = next_uint64 t in - Uint64.shift_right uint 11 - |> Uint64.to_string - |> Float.of_string - |> Float.mul (1.0 /. 9007199254740992.0), t' + let next_double t = Common.next_double ~nextu64:next_uint64 t let jump = Uint64.( @@ -103,7 +72,5 @@ end = struct let initialize seed = let istate = Seed.SeedSequence.generate_64bit_state 4 seed in - {s = (istate.(0), istate.(1), istate.(2), istate.(3)); - has_uint32 = false; - uinteger = Uint32.zero} + {s = (istate.(0), istate.(1), istate.(2), istate.(3)); ustore = Common.Empty} end