anchor.diff:  %% uniform/1 returns a random integer X where 1 =< X =< N,
   %% updating the state in the process dictionary.
   
  --spec uniform(N :: pos_integer()) -> X::pos_integer().
  +-spec uniform(N :: pos_integer()) -> X :: pos_integer().
   uniform(N) ->
       {X, Seed} = uniform_s(N, seed_get()),
       _ = seed_put(Seed),
       X.
   
   %% uniform_s/1: given a state, uniform_s/1
  -%% returns a random float X where 0.0 < X < 1.0,
  +%% returns a random float X where 0.0 =< X < 1.0,
   %% and a new state.
   
  --spec uniform_s(state()) -> {X::float(), NewS :: state()}.
  +-spec uniform_s(State :: state()) -> {X :: float(), NewState :: state()}.
   uniform_s(State = {#{uniform:=Uniform}, _}) ->
  -    Uniform(State).
  +    Uniform(State);
  +uniform_s({#{bits:=Bits, next:=Next} = Alg, R0}) ->
  +    {V, R1} = Next(R0),
  +    %% Produce floats on the form N * 2^(-53)
  +    {(V bsr (Bits - 53)) * ?TWO_POW_MINUS53, {Alg, R1}};
  +uniform_s({#{max:=Max, next:=Next} = Alg, R0}) ->
  +    {V, R1} = Next(R0),
  +    %% Old broken algorithm with non-uniform density
  +    {V / (Max + 1), {Alg, R1}}.
  +
   
   %% uniform_s/2: given an integer N >= 1 and a state, uniform_s/2
   %% uniform_s/2 returns a random integer X where 1 =< X =< N,
   %% and a new state.
   
  --spec uniform_s(N::pos_integer(), state()) -> {X::pos_integer(), NewS::state()}.
  -uniform_s(N, State = {#{uniform_n:=Uniform, max:=Max}, _})
  -  when 0 < N, N =< Max ->
  -    Uniform(N, State);
  -uniform_s(N, State0 = {#{uniform:=Uniform}, _})
  -  when is_integer(N), 0 < N ->
  -    {F, State} = Uniform(State0),
  -    {trunc(F * N) + 1, State}.
  +-spec uniform_s(N :: pos_integer(), State :: state()) ->
  +                       {X :: pos_integer(), NewState :: state()}.
  +uniform_s(N, State = {#{uniform_n:=UniformN}, _})
  +  when is_integer(N), 1 =< N ->
  +    UniformN(N, State);
  +uniform_s(N, {#{bits:=Bits, next:=Next} = Alg, R0})
  +  when is_integer(N), 1 =< N ->
  +    {V, R1} = Next(R0),
  +    MaxMinusN = ?BIT(Bits) - N,
  +    ?uniform_range(N, Alg, R1, V, MaxMinusN, I);
  +uniform_s(N, {#{max:=Max, next:=Next} = Alg, R0})
  +  when is_integer(N), 1 =< N ->
  +    %% Old broken algorithm with skewed probability
  +    %% and gap in ranges > Max
  +    {V, R1} = Next(R0),  
  +    if
  +        N =< Max ->
  +            {(V rem N) + 1, {Alg, R1}};
  +        true ->
  +            F = V / (Max + 1),
  +            {trunc(F * N) + 1, {Alg, R1}}
  +    end.
   
   %% jump/1: given a state, jump/1
   %% returns a new state which is equivalent to that
   %% after a large number of call defined for each algorithm.
   %% The large number is algorithm dependent.
   
  --spec jump(state()) -> NewS :: state().
  +-spec jump(state()) -> NewState :: state().
   jump(State = {#{jump:=Jump}, _}) ->
  -    Jump(State).
  +    Jump(State);
  +jump({#{}, _}) ->
  +    erlang:error(not_implemented).
  +
   
   %% jump/0: read the internal state and
   %% apply the jump function for the state as in jump/1
   %% and write back the new value to the internal state,
   %% then returns the new value.
   
  --spec jump() -> NewS :: state().
  -
  +-spec jump() -> NewState :: state().
   jump() ->
       seed_put(jump(seed_get())).
   
  
positive.code:     NL = [H|lists:reverse(RL)],
      exs1024_next({NL, []}).
  
  exs1024_uniform({Alg, R0}) ->
      {V, R1} = exs1024_next(R0),
      {V / 18446744073709551616, {Alg, R1}}.
  
  exs1024_uniform(Max, {Alg, R}) ->
      {V, R1} = exs1024_next(R),
      {(V rem Max) + 1, {Alg, R1}}.
  
  %% This is the jump function for the exs1024 generator, equivalent
  %% to 2^512 calls to next(); it can be used to generate 2^512
  
negative.code:     _ = seed_put(Seed),
      X.
  
  %% normal/2: returns a random float with N(μ, σ²) normal distribution
  %% updating the state in the process dictionary.
  
  -spec normal(Mean :: number(), Variance :: number()) -> float().
  normal(Mean, Variance) ->
      Mean + (math:sqrt(Variance) * normal()).
  
  %% normal_s/1: returns a random float with standard normal distribution
  %% The Ziggurat Method for generating random variables - Marsaglia and Tsang
  %% Paper and reference code: http://www.jstatsoft.org/v05/i08/
  
  -spec normal_s(State :: state()) -> {float(), NewState :: state()}.
  normal_s(State0) ->
      {Sign, R, State} = get_52(State0),
      Idx = ?MASK(8, R),
      Idx1 = Idx+1,
      {Ki, Wi} = normal_kiwi(Idx1),
      X = R * Wi,
  