<commit_before>%%
%% %CopyrightBegin%
%%
%% Copyright Ericsson AB 2000-2017. All Rights Reserved.
%%
%% Licensed under the Apache License, Version 2.0 (the "License");
%% you may not use this file except in compliance with the License.
%% You may obtain a copy of the License at
%%
%%     http://www.apache.org/licenses/LICENSE-2.0
%%
%% Unless required by applicable law or agreed to in writing, software
%% distributed under the License is distributed on an "AS IS" BASIS,
%% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
%% See the License for the specific language governing permissions and
%% limitations under the License.
%%
%% %CopyrightEnd%

-module(rand_SUITE).
-compile({nowarn_deprecated_function,[{random,seed,1},
                                      {random,uniform_s,1},
                                      {random,uniform_s,2}]}).

-export([all/0, suite/0, groups/0, group/1]).

-export([interval_int/1, interval_float/1, seed/1,
         api_eq/1, reference/1,
	 basic_stats_uniform_1/1, basic_stats_uniform_2/1,
	 basic_stats_standard_normal/1,
	 basic_stats_normal/1,
         uniform_real_conv/1,
	 plugin/1, measure/1,
	 reference_jump_state/1, reference_jump_procdict/1]).

-export([test/0, gen/1]).

-export([uniform_real_gen/1, uniform_gen/2]).

-include_lib("common_test/include/ct.hrl").

-define(LOOP, 1000000).

suite() ->
    [{ct_hooks,[ts_install_cth]},
     {timetrap,{minutes,3}}].

all() ->
    [seed, interval_int, interval_float,
     api_eq,
     reference,
     {group, basic_stats}, uniform_real_conv,
     plugin, measure,
     {group, reference_jump}
    ].

groups() ->
    [{basic_stats, [parallel],
      [basic_stats_uniform_1, basic_stats_uniform_2,
       basic_stats_standard_normal]},
     {reference_jump, [parallel],
      [reference_jump_state, reference_jump_procdict]}].

group(basic_stats) ->
    %% valgrind needs a lot of time
    [{timetrap,{minutes,10}}];
group(reference_jump) ->
    %% valgrind needs a lot of time
    [{timetrap,{minutes,10}}].

%% A simple helper to test without test_server during dev
test() ->
    Tests = all(),
    lists:foreach(
      fun (Test) ->
              try
                  ok = ?MODULE:Test([]),
                  io:format("~p: ok~n", [Test])
              catch _:Reason ->
                      io:format("Failed: ~p: ~p ~p~n",
                                [Test, Reason, erlang:get_stacktrace()])
              end
      end, Tests).

algs() ->
    [exrop, exsp, exs1024s, exs64, exsplus, exs1024].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Test that seed and seed_s and export_seed/0 is working.
seed(Config) when is_list(Config) ->
    Algs = algs(),
    Test = fun(Alg) ->
		   try seed_1(Alg)
		   catch _:Reason ->
			   ct:fail({Alg, Reason, erlang:get_stacktrace()})
		   end
	   end,
    [Test(Alg) || Alg <- Algs],
    ok.

seed_1(Alg) ->
    %% Check that uniform seeds automatically,
    _ = rand:uniform(),
    S00 = get(rand_seed),
    erase(),
    _ = rand:uniform_real(),
    false = S00 =:= get(rand_seed), %% hopefully

    %% Choosing algo and seed
    S0 = rand:seed(Alg, {0, 0, 0}),
    %% Check that (documented?) process_dict variable is correct
    S0 = get(rand_seed),
    S0 = rand:seed_s(Alg, {0, 0, 0}),
    %% Check that process_dict should not be used for seed_s functionality
    _ = rand:seed_s(Alg, {1, 0, 0}),
    S0 = get(rand_seed),
    %% Test export
    ES0 = rand:export_seed(),
    ES0 = rand:export_seed_s(S0),
    S0 = rand:seed(ES0),
    S0 = rand:seed_s(ES0),
    %% seed/1 calls should be unique
    S1 = rand:seed(Alg),
    false = (S1 =:= rand:seed_s(Alg)),
    %% Negative integers works
    _ = rand:seed_s(Alg, {-1,-1,-1}),
    %% Check that export_seed/1 returns 'undefined' if there is no seed
    erase(rand_seed),
    undefined = rand:export_seed(),

    %% Other term do not work
    {'EXIT', _} = (catch rand:seed_s(foobar, os:timestamp())),
    {'EXIT', _} = (catch rand:seed_s(Alg, {asd, 1, 1})),
    {'EXIT', _} = (catch rand:seed_s(Alg, {0, 234.1234, 1})),
    {'EXIT', _} = (catch rand:seed_s(Alg, {0, 234, [1, 123, 123]})),
    ok.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check that both APIs are consistent with each other.
api_eq(_Config) ->
    Algs = algs(),
    Small = fun(Alg) ->
		    Seed = rand:seed(Alg),
		    io:format("Seed ~p~n",[rand:export_seed_s(Seed)]),
		    api_eq_1(Seed)
	    end,
    _ = [Small(Alg) || Alg <- Algs],
    ok.

api_eq_1(S00) ->
    Check = fun(_, Seed) ->
		    {V0, S0} = rand:uniform_s(Seed),
		    V0 = rand:uniform(),
		    {V1, S1} = rand:uniform_s(1000000, S0),
		    V1 = rand:uniform(1000000),
		    {V2, S2} = rand:normal_s(S1),
		    V2 = rand:normal(),
		    S2
	    end,
    S1 = lists:foldl(Check, S00, lists:seq(1, 200)),
    S1 = get(rand_seed),
    {V0, S2} = rand:uniform_s(S1),
    V0 = rand:uniform(),
    S2 = get(rand_seed),

    Exported = rand:export_seed(),
    Exported = rand:export_seed_s(S2),

    S3 = lists:foldl(Check, S2, lists:seq(1, 200)),
    S3 = get(rand_seed),

    S4 = lists:foldl(Check, S3, lists:seq(1, 200)),
    S4 = get(rand_seed),
    %% Verify that we do not have loops
    false = S1 =:= S2,
    false = S2 =:= S3,
    false = S3 =:= S4,

    S2 = rand:seed(Exported),
    S3 = lists:foldl(Check, S2, lists:seq(1, 200)),
    ok.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check that uniform/1 returns values within the proper interval.
interval_int(Config) when is_list(Config) ->
    Algs = algs(),
    Small = fun(Alg) ->
		    Seed = rand:seed(Alg),
		    io:format("Seed ~p~n",[rand:export_seed_s(Seed)]),
		    Max = interval_int_1(100000, 7, 0),
		    Max =:= 7 orelse exit({7, Alg, Max})
	    end,
    _ = [Small(Alg) || Alg <- Algs],
    %% Test large integers
    Large = fun(Alg) ->
		    Seed = rand:seed(Alg),
		    io:format("Seed ~p~n",[rand:export_seed_s(Seed)]),
		    Max = interval_int_1(100000, 1 bsl 128, 0),
		    Max > 1 bsl 64 orelse exit({large, Alg, Max})
	    end,
    [Large(Alg) || Alg <- Algs],
    ok.

interval_int_1(0, _, Max) -> Max;
interval_int_1(N, Top, Max) ->
    X = rand:uniform(Top),
    if
	0 < X, X =< Top ->
	    ok;
	true ->
	    io:format("X=~p Top=~p 0<~p<~p~n", [X,Top,X,Top]),
	    exit({X, rand:export_seed()})
    end,
    interval_int_1(N-1, Top, max(X, Max)).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check that uniform/0 returns values within the proper interval.
interval_float(Config) when is_list(Config) ->
    Algs = algs(),
    Test = fun(Alg) ->
		   _ = rand:seed(Alg),
		   interval_float_1(100000)
	   end,
    [Test(Alg) || Alg <- Algs],
    ok.

interval_float_1(0) -> ok;
interval_float_1(N) ->
    X = rand:uniform(),
    Y = rand:uniform_real(),
    if
	0.0 =< X, X < 1.0, 0.0 < Y, Y < 1.0 ->
	    ok;
	true ->
	    io:format("X=~p 0.0=<~p<1.0~n", [X,X]),
	    io:format("Y=~p 0.0<~p<1.0~n", [Y,Y]),
	    exit({X, rand:export_seed()})
    end,
    interval_float_1(N-1).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check if each algorithm generates the proper sequence.
reference(Config) when is_list(Config) ->
    [reference_1(Alg) || Alg <- algs()],
    ok.

reference_1(Alg) ->
    Refval  = reference_val(Alg),
    Testval = gen(Alg),
    case Refval =:= Testval of
        true -> ok;
        false when Refval =:= not_implemented ->
            exit({not_implemented,Alg});
        false ->
	    io:format("Failed: ~p~n",[Alg]),
	    io:format("Length ~p ~p~n",[length(Refval), length(Testval)]),
	    io:format("Head ~p ~p~n",[hd(Refval), hd(Testval)]),
	    exit(wrong_value)
    end.

gen(Algo) ->
    State =
        case Algo of
            exs64 -> %% Printed with orig 'C' code and this seed
                rand:seed_s({exs64, 12345678});
            _ when Algo =:= exsplus; Algo =:= exsp; Algo =:= exrop ->
                %% Printed with orig 'C' code and this seed
                rand:seed_s({Algo, [12345678|12345678]});
            _ when Algo =:= exs1024; Algo =:= exs1024s ->
                %% Printed with orig 'C' code and this seed
                rand:seed_s({Algo, {lists:duplicate(16, 12345678), []}});
            _ ->
                rand:seed(Algo, {100, 200, 300})
        end,
    Max = range(State),
    gen(?LOOP, State, Max, []).

gen(N, State0, Max, Acc) when N > 0 ->
    {Random, State} = rand:uniform_s(Max, State0),
    case N rem (?LOOP div 100) of
	0 -> gen(N-1, State, Max, [Random|Acc]);
	_ -> gen(N-1, State, Max, Acc)
    end;
gen(_, _, _, Acc) -> lists:reverse(Acc).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This just tests the basics so we have not made any serious errors
%% when making the conversion from the original algorithms.
%% The algorithms must have good properties to begin with
%%

%% Check that the algorithms generate sound values.

basic_stats_uniform_1(Config) when is_list(Config) ->
    ct:timetrap({minutes,15}), %% valgrind needs a lot of time
    [basic_uniform_1(?LOOP, rand:seed_s(Alg), 0.0, array:new([{default, 0}]))
     || Alg <- algs()],
    ok.

basic_stats_uniform_2(Config) when is_list(Config) ->
    ct:timetrap({minutes,15}), %% valgrind needs a lot of time
    [basic_uniform_2(?LOOP, rand:seed_s(Alg), 0, array:new([{default, 0}]))
     || Alg <- algs()],
    ok.

basic_stats_standard_normal(Config) when is_list(Config) ->
    ct:timetrap({minutes,6}), %% valgrind needs a lot of time
    io:format("Testing standard normal~n",[]),
    IntendedMean = 0,
    IntendedVariance = 1,
    [basic_normal_1(?LOOP, IntendedMean, IntendedVariance,
                    rand:seed_s(Alg), 0, 0)
     || Alg <- algs()],
    ok.

basic_stats_normal(Config) when is_list(Config) ->
    IntendedMeans = [-1.0e6, -50, -math:pi(), -math:exp(-1),
                     0.12345678, math:exp(1), 100, 1.0e6],
    IntendedVariances = [1.0e-6, math:exp(-1), 1, math:pi(), 1.0e6],
    IntendedMeanVariancePairs =
        [{Mean, Variance} || Mean <- IntendedMeans,
                             Variance <- IntendedVariances],

    ct:timetrap({minutes, 6 * length(IntendedMeanVariancePairs)}), %% valgrind needs a lot of time
    lists:foreach(
      fun ({IntendedMean, IntendedVariance}) ->
              ct:pal(
                "Testing normal(~.2f, ~.2f)~n",
                [float(IntendedMean), float(IntendedVariance)]),
              [basic_normal_1(?LOOP, IntendedMean, IntendedVariance,
                              rand:seed_s(Alg), 0, 0)
               || Alg <- algs()]
      end,
      IntendedMeanVariancePairs).

basic_uniform_1(N, S0, Sum, A0) when N > 0 ->
    {X,S} =
        case N band 1 of
            0 ->
                rand:uniform_s(S0);
            1 ->
                rand:uniform_real_s(S0)
        end,
    I = trunc(X*100),
    A = array:set(I, 1+array:get(I,A0), A0),
    basic_uniform_1(N-1, S, Sum+X, A);
basic_uniform_1(0, {#{type:=Alg}, _}, Sum, A) ->
    AverN = Sum / ?LOOP,
    io:format("~.12w: Average: ~.4f~n", [Alg, AverN]),
    Counters = array:to_list(A),
    Min = lists:min(Counters),
    Max = lists:max(Counters),
    io:format("~.12w: Min: ~p Max: ~p~n", [Alg, Min, Max]),

    %% Verify that the basic statistics are ok
    %% be gentle we don't want to see to many failing tests
    abs(0.5 - AverN) < 0.005 orelse ct:fail({average, Alg, AverN}),
    abs(?LOOP div 100 - Min) < 1000 orelse ct:fail({min, Alg, Min}),
    abs(?LOOP div 100 - Max) < 1000 orelse ct:fail({max, Alg, Max}),
    ok.

basic_uniform_2(N, S0, Sum, A0) when N > 0 ->
    {X,S} = rand:uniform_s(100, S0),
    A = array:set(X-1, 1+array:get(X-1,A0), A0),
    basic_uniform_2(N-1, S, Sum+X, A);
basic_uniform_2(0, {#{type:=Alg}, _}, Sum, A) ->
    AverN = Sum / ?LOOP,
    io:format("~.12w: Average: ~.4f~n", [Alg, AverN]),
    Counters = tl(array:to_list(A)),
    Min = lists:min(Counters),
    Max = lists:max(Counters),
    io:format("~.12w: Min: ~p Max: ~p~n", [Alg, Min, Max]),

    %% Verify that the basic statistics are ok
    %% be gentle we don't want to see to many failing tests
    abs(50.5 - AverN) < 0.5 orelse ct:fail({average, Alg, AverN}),
    abs(?LOOP div 100 - Min) < 1000 orelse ct:fail({min, Alg, Min}),
    abs(?LOOP div 100 - Max) < 1000 orelse ct:fail({max, Alg, Max}),
    ok.

basic_normal_1(N, IntendedMean, IntendedVariance, S0, StandardSum, StandardSq) when N > 0 ->
    {X,S} = normal_s(IntendedMean, IntendedVariance, S0),
    % We now shape X into a standard normal distribution (in case it wasn't already)
    % in order to minimise the accumulated error on Sum / SumSq;
    % otherwise said error would prevent us of making a fair judgment on
    % the overall distribution when targeting large means and variances.
    StandardX = (X - IntendedMean) / math:sqrt(IntendedVariance),
    basic_normal_1(N-1, IntendedMean, IntendedVariance, S,
                   StandardX+StandardSum, StandardX*StandardX+StandardSq);
basic_normal_1(0, _IntendedMean, _IntendedVariance, {#{type:=Alg}, _}, StandardSum, StandardSumSq) ->
    StandardMean = StandardSum / ?LOOP,
    StandardVariance = (StandardSumSq - (StandardSum*StandardSum/?LOOP))/(?LOOP - 1),
    StandardStdDev =  math:sqrt(StandardVariance),
    io:format("~.12w: Standardised Average: ~7.4f, Standardised StdDev ~6.4f~n",
              [Alg, StandardMean, StandardStdDev]),
    %% Verify that the basic statistics are ok
    %% be gentle we don't want to see to many failing tests
    abs(StandardMean) < 0.005 orelse ct:fail({average, Alg, StandardMean}),
    abs(StandardStdDev - 1.0) < 0.005 orelse ct:fail({stddev, Alg, StandardStdDev}),
    ok.

normal_s(Mean, Variance, State0) when Mean == 0, Variance == 1 ->
    % Make sure we're also testing the standard normal interface
    rand:normal_s(State0);
normal_s(Mean, Variance, State0) ->
    rand:normal_s(Mean, Variance, State0).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% White box test of the conversion to float

uniform_real_conv(Config) when is_list(Config) ->
    [begin
%%         ct:pal("~13.16.0bx~3.16.0b: ~p~n", [M,E,Gen]),
         uniform_real_conv_check(M, E, Gen)
     end || {M, E, Gen} <- uniform_real_conv_data()],
    uniform_real_scan(0),
    uniform_real_scan(3).

uniform_real_conv_data() ->
    [{16#fffffffffffff,  -1, [16#3ffffffffffffff]},
     {16#fffffffffffff,  -1, [16#3ffffffffffffe0]},
     {16#ffffffffffffe,  -1, [16#3ffffffffffffdf]},
     %%
     {16#0000000000000,  -1, [16#200000000000000]},
     {16#fffffffffffff,  -2, [16#1ffffffffffffff]},
     {16#fffffffffffff,  -2, [16#1fffffffffffff0]},
     {16#ffffffffffffe,  -2, [16#1ffffffffffffef]},
     %%
     {16#0000000000000,  -2, [16#100000000000000]},
     {16#fffffffffffff,  -3, [16#0ffffffffffffff]},
     {16#fffffffffffff,  -3, [16#0fffffffffffff8]},
     {16#ffffffffffffe,  -3, [16#0fffffffffffff7]},
     %%
     {16#0000000000000,  -3, [16#080000000000000]},
     {16#fffffffffffff,  -4, [16#07fffffffffffff]},
     {16#fffffffffffff,  -4, [16#07ffffffffffffc]},
     {16#ffffffffffffe,  -4, [16#07ffffffffffffb]},
     %%
     {16#0000000000000,  -4, [16#040000000000000]},
     {16#fffffffffffff,  -5, [16#03fffffffffffff,16#3ffffffffffffff]},
     {16#fffffffffffff,  -5, [16#03ffffffffffffe,16#200000000000000]},
     {16#ffffffffffffe,  -5, [16#03fffffffffffff,16#1ffffffffffffff]},
     {16#ffffffffffffe,  -5, [16#03fffffffffffff,16#100000000000000]},
     %%
     {16#0000000000001, -56, [16#000000000000007,16#00000000000007f]},
     {16#0000000000001, -56, [16#000000000000004,16#000000000000040]},
     {16#0000000000000, -57, [16#000000000000003,16#20000000000001f]},
     {16#0000000000000, -57, [16#000000000000000,16#200000000000000]},
     {16#fffffffffffff, -58, [16#000000000000003,16#1ffffffffffffff]},
     {16#fffffffffffff, -58, [16#000000000000000,16#1fffffffffffff0]},
     {16#ffffffffffffe, -58, [16#000000000000000,16#1ffffffffffffef]},
     {16#ffffffffffffe, -58, [16#000000000000000,16#1ffffffffffffe0]},
     %%
     {16#0000000000000, -58, [16#000000000000000,16#10000000000000f]},
     {16#0000000000000, -58, [16#000000000000000,16#100000000000000]},
     {2#11001100000000000000000000000000000000000011000000011, % 53 bits
      -1022,
      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, % 18 zeros
       2#1100110000000000000000000000000000000000001 bsl 2, % 43 bits
       2#1000000011 bsl (56-10+2)]}, % 10 bits
     {0, -1, % 0.5 after retry
      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, % 18 zeros
       2#111111111111111111111111111111111111111111 bsl 2, % 42 bits - retry
       16#200000000000003]}]. % 0.5

-define(UNIFORM_REAL_SCAN_PATTERN, (16#19000000000009)). % 53 bits
-define(UNIFORM_REAL_SCAN_NUMBER, (1021)).

uniform_real_scan_template(K) ->
    <<0:?UNIFORM_REAL_SCAN_NUMBER,
      ?UNIFORM_REAL_SCAN_PATTERN:53,K:2,0:1>>.

uniform_real_scan(K) ->
    Templ = uniform_real_scan_template(K),
    N = ?UNIFORM_REAL_SCAN_NUMBER,
    uniform_real_scan(Templ, N, K).

uniform_real_scan(Templ, N, K) when 0 =< N ->
    <<_:N/bits,T/bits>> = Templ,
    Data = uniform_real_scan_data(T, K),
    uniform_real_conv_check(
      ?UNIFORM_REAL_SCAN_PATTERN, N - 1 - ?UNIFORM_REAL_SCAN_NUMBER, Data),
    uniform_real_scan(Templ, N - 1, K);
uniform_real_scan(_, _, _) ->
    ok.

uniform_real_scan_data(Templ, K) ->
    case Templ of
        <<X:56, T/bits>> ->
            B = rand:bc64(X),
            [(X bsl 2) bor K |
             if
                 53 =< B ->
                     [];
                 true ->
                     uniform_real_scan_data(T, K)
             end];
        _ ->
            <<X:56, _/bits>> = <<Templ/bits, 0:56>>,
            [(X bsl 2) bor K]
    end.

uniform_real_conv_check(M, E, Gen) ->
    <<F/float>> = <<0:1, (E + 16#3ff):11, M:52>>,
    try uniform_real_gen(Gen) of
        F -> F;
        FF ->
            ct:pal(
              "~s =/= ~s: ~s~n",
              [rand:float2str(FF), rand:float2str(F),
               [["16#",integer_to_list(G,16),$\s]||G<-Gen]]),
            ct:fail({neq, FF, F})
    catch
        Error:Reason ->
            ct:pal(
              "~w:~p ~s: ~s~n",
              [Error, Reason, rand:float2str(F),
               [["16#",integer_to_list(G,16),$\s]||G<-Gen]]),
            ct:fail({Error, Reason, F, erlang:get_stacktrace()})
    end.


uniform_real_gen(Gen) ->
    State = rand_state(Gen),
    {F, {#{type := rand_SUITE_list},[]}} = rand:uniform_real_s(State),
    F.

uniform_gen(Range, Gen) ->
    State = rand_state(Gen),
    {N, {#{type := rand_SUITE_list},[]}} = rand:uniform_s(Range, State),
    N.

%% Loaded dice for white box tests
rand_state(Gen) ->
    {#{type => rand_SUITE_list, bits => 58, weak_low_bits => 1,
       next => fun ([H|T]) -> {H, T} end},
     Gen}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Test that the user can write algorithms.
plugin(Config) when is_list(Config) ->
    try crypto:strong_rand_bytes(1) of
        <<_>> ->
            _ = lists:foldl(
                  fun(_, S0) ->
                          {V1, S1} = rand:uniform_s(10000, S0),
                          true = is_integer(V1),
                          {V2, S2} = rand:uniform_s(S1),
                          true = is_float(V2),
                          S2
                  end, crypto64_seed(), lists:seq(1, 200)),
            ok
    catch
        error:low_entropy ->
            {skip,low_entropy};
        error:undef ->
            {skip,no_crypto}
    end.

%% Test implementation
crypto64_seed() ->
    {#{type=>crypto64,
       bits=>64,
       next=>fun crypto64_next/1,
       uniform=>fun crypto64_uniform/1,
       uniform_n=>fun crypto64_uniform_n/2},
     <<>>}.

%% Be fair and create bignums i.e. 64bits otherwise use 58bits
crypto64_next(<<Num:64, Bin/binary>>) ->
    {Num, Bin};
crypto64_next(_) ->
    crypto64_next(crypto:strong_rand_bytes((64 div 8)*100)).

crypto64_uniform({Api, Data0}) ->
    {Int, Data} = crypto64_next(Data0),
    {Int / (1 bsl 64), {Api, Data}}.

crypto64_uniform_n(N, {Api, Data0}) when N < (1 bsl 64) ->
    {Int, Data} = crypto64_next(Data0),
    {(Int rem N)+1, {Api, Data}};
crypto64_uniform_n(N, State0) ->
    {F,State} = crypto64_uniform(State0),
    {trunc(F * N) + 1, State}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Not a test but measures the time characteristics of the different algorithms
measure(Config) ->
    ct:timetrap({minutes,60}), %% valgrind needs a lot of time
    case ct:get_timetrap_info() of
        {_,{_,1}} -> % No scaling
            do_measure(Config);
        {_,{_,Scale}} ->
            {skip,{will_not_run_in_scaled_time,Scale}}
    end.

-define(CHECK_UNIFORM_RANGE(Gen, Range, X, St),
        case (Gen) of
            {(X), (St)} when is_integer(X), 1 =< (X), (X) =< (Range) ->
                St
        end).
-define(CHECK_UNIFORM(Gen, X, St),
        case (Gen) of
            {(X), (St)} when is_float(X), 0.0 =< (X), (X) < 1.0 ->
                St
        end).
-define(CHECK_UNIFORM_NZ(Gen, X, St),
        case (Gen) of
            {(X), (St)} when is_float(X), 0.0 < (X), (X) =< 1.0 ->
                St
        end).
-define(CHECK_NORMAL(Gen, X, St),
        case (Gen) of
            {(X), (St)} when is_float(X) ->
                St
        end).

do_measure(_Config) ->
    Algs =
        algs() ++
        try crypto:strong_rand_bytes(1) of
            <<_>> -> [crypto64, crypto_cache, crypto]
        catch
            error:low_entropy -> [];
            error:undef -> []
        end,
    %%
    ct:pal("~nRNG uniform integer range 10000 performance~n",[]),
    _ =
        measure_1(
          fun (_) -> 10000 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer 32 bit performance~n",[]),
    _ =
        measure_1(
          fun (_) -> 1 bsl 32 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer half range performance~n",[]),
    _ =
        measure_1(
          fun (State) -> half_range(State) end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer half range + 1 performance~n",[]),
    _ =
        measure_1(
          fun (State) -> half_range(State) + 1 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer full range - 1 performance~n",[]),
    _ =
        measure_1(
          fun (State) -> (half_range(State) bsl 1) - 1 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer full range performance~n",[]),
    _ =
        measure_1(
          fun (State) -> half_range(State) bsl 1 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer full range + 1 performance~n",[]),
    _ =
        measure_1(
          fun (State) -> (half_range(State) bsl 1) + 1 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer double range performance~n",[]),
    _ =
        measure_1(
          fun (State) ->
                  half_range(State) bsl 2
          end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer double range + 1  performance~n",[]),
    _ =
        measure_1(
          fun (State) ->
                  (half_range(State) bsl 2) + 1
          end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer 64 bit performance~n",[]),
    _ =
        measure_1(
          fun (_) -> 1 bsl 64 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform float performance~n",[]),
    _ =
        measure_1(
          fun (_) -> 0 end,
          fun (State, _, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM(Mod:uniform_s(St0), X, St)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform_real float performance~n",[]),
    _ =
        measure_1(
          fun (_) -> 0 end,
          fun (State, _, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM(Mod:uniform_real_s(St0), X, St)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG normal float performance~n",[]),
    [TMarkNormalFloat|_] =
        measure_1(
          fun (_) -> 0 end,
          fun (State, _, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_NORMAL(Mod:normal_s(St0), X, St1)
                    end,
                    State)
          end,
          Algs),
    %% Just for fun try an implementation of the Box-Muller
    %% transformation for creating normal distribution floats
    %% to compare with our Ziggurat implementation.
    %% Generates two numbers per call that we add so they
    %% will not be optimized away.  Hence the benchmark time
    %% is twice what it should be.
    TwoPi = 2 * math:pi(),
    _ =
        measure_1(
          fun (_) -> 0 end,
          fun (State, _, Mod) ->
                  measure_loop(
                    fun (State0) ->
                            {U1, State1} = Mod:uniform_real_s(State0),
                            {U2, State2} = Mod:uniform_s(State1),
                            R = math:sqrt(-2.0 * math:log(U1)),
                            T = TwoPi * U2,
                            Z0 = R * math:cos(T),
                            Z1 = R * math:sin(T),
                            ?CHECK_NORMAL({Z0 + Z1, State2}, X, State3)
                    end,
                    State)
          end,
          exrop, TMarkNormalFloat),
    ok.

measure_loop(Fun, State) ->
    measure_loop(Fun, State, ?LOOP).
%%
measure_loop(Fun, State, N) when 0 < N ->
    measure_loop(Fun, Fun(State), N-1);
measure_loop(_, _, _) ->
    ok.

measure_1(RangeFun, Fun, Algs) ->
    TMark = measure_1(RangeFun, Fun, hd(Algs), undefined),
    [TMark] ++
        [measure_1(RangeFun, Fun, Alg, TMark) || Alg <- tl(Algs)].

measure_1(RangeFun, Fun, Alg, TMark) ->
    Parent = self(),
    {Mod, State} =
        case Alg of
            crypto64 ->
                {rand, crypto64_seed()};
            crypto_cache ->
                {rand, crypto:rand_seed_alg(crypto_cache)};
            crypto ->
                {rand, crypto:rand_seed_s()};
            random ->
                {random, random:seed(os:timestamp()), get(random_seed)};
            _ ->
                {rand, rand:seed_s(Alg)}
        end,
    Range = RangeFun(State),
    Pid = spawn_link(
            fun() ->
                    {Time, ok} = timer:tc(fun () -> Fun(State, Range, Mod) end),
                    Percent =
                        case TMark of
                            undefined -> 100;
                            _ -> (Time * 100 + 50) div TMark
                        end,
                    io:format(
                      "~.12w: ~p ns ~p% [16#~.16b]~n",
                      [Alg, (Time * 1000 + 500) div ?LOOP, Percent, Range]),
                    Parent ! {self(), Time},
                    normal
            end),
    receive
	{Pid, Msg} -> Msg
    end.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The jump sequence tests has two parts
%% for those with the functional API (jump/1)
%% and for those with the internal state
%% in process dictionary (jump/0).

-define(LOOP_JUMP, (?LOOP div 1000)).

%% Check if each algorithm generates the proper jump sequence
%% with the functional API.
reference_jump_state(Config) when is_list(Config) ->
    [reference_jump_1(Alg) || Alg <- algs()],
    ok.

reference_jump_1(Alg) ->
    Refval  = reference_jump_val(Alg),
    Testval = gen_jump_1(Alg),
    case Refval =:= Testval of
        true -> ok;
        false -><commit_msg>Scale down measure rounds<commit_after>%%
%% %CopyrightBegin%
%%
%% Copyright Ericsson AB 2000-2017. All Rights Reserved.
%%
%% Licensed under the Apache License, Version 2.0 (the "License");
%% you may not use this file except in compliance with the License.
%% You may obtain a copy of the License at
%%
%%     http://www.apache.org/licenses/LICENSE-2.0
%%
%% Unless required by applicable law or agreed to in writing, software
%% distributed under the License is distributed on an "AS IS" BASIS,
%% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
%% See the License for the specific language governing permissions and
%% limitations under the License.
%%
%% %CopyrightEnd%

-module(rand_SUITE).
-compile({nowarn_deprecated_function,[{random,seed,1},
                                      {random,uniform_s,1},
                                      {random,uniform_s,2}]}).

-export([all/0, suite/0, groups/0, group/1]).

-export([interval_int/1, interval_float/1, seed/1,
         api_eq/1, reference/1,
	 basic_stats_uniform_1/1, basic_stats_uniform_2/1,
	 basic_stats_standard_normal/1,
	 basic_stats_normal/1,
         uniform_real_conv/1,
	 plugin/1, measure/1,
	 reference_jump_state/1, reference_jump_procdict/1]).

-export([test/0, gen/1]).

-export([uniform_real_gen/1, uniform_gen/2]).

-include_lib("common_test/include/ct.hrl").

-define(LOOP, 1000000).

suite() ->
    [{ct_hooks,[ts_install_cth]},
     {timetrap,{minutes,3}}].

all() ->
    [seed, interval_int, interval_float,
     api_eq,
     reference,
     {group, basic_stats}, uniform_real_conv,
     plugin, measure,
     {group, reference_jump}
    ].

groups() ->
    [{basic_stats, [parallel],
      [basic_stats_uniform_1, basic_stats_uniform_2,
       basic_stats_standard_normal]},
     {reference_jump, [parallel],
      [reference_jump_state, reference_jump_procdict]}].

group(basic_stats) ->
    %% valgrind needs a lot of time
    [{timetrap,{minutes,10}}];
group(reference_jump) ->
    %% valgrind needs a lot of time
    [{timetrap,{minutes,10}}].

%% A simple helper to test without test_server during dev
test() ->
    Tests = all(),
    lists:foreach(
      fun (Test) ->
              try
                  ok = ?MODULE:Test([]),
                  io:format("~p: ok~n", [Test])
              catch _:Reason ->
                      io:format("Failed: ~p: ~p ~p~n",
                                [Test, Reason, erlang:get_stacktrace()])
              end
      end, Tests).

algs() ->
    [exrop, exsp, exs1024s, exs64, exsplus, exs1024].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Test that seed and seed_s and export_seed/0 is working.
seed(Config) when is_list(Config) ->
    Algs = algs(),
    Test = fun(Alg) ->
		   try seed_1(Alg)
		   catch _:Reason ->
			   ct:fail({Alg, Reason, erlang:get_stacktrace()})
		   end
	   end,
    [Test(Alg) || Alg <- Algs],
    ok.

seed_1(Alg) ->
    %% Check that uniform seeds automatically,
    _ = rand:uniform(),
    S00 = get(rand_seed),
    erase(),
    _ = rand:uniform_real(),
    false = S00 =:= get(rand_seed), %% hopefully

    %% Choosing algo and seed
    S0 = rand:seed(Alg, {0, 0, 0}),
    %% Check that (documented?) process_dict variable is correct
    S0 = get(rand_seed),
    S0 = rand:seed_s(Alg, {0, 0, 0}),
    %% Check that process_dict should not be used for seed_s functionality
    _ = rand:seed_s(Alg, {1, 0, 0}),
    S0 = get(rand_seed),
    %% Test export
    ES0 = rand:export_seed(),
    ES0 = rand:export_seed_s(S0),
    S0 = rand:seed(ES0),
    S0 = rand:seed_s(ES0),
    %% seed/1 calls should be unique
    S1 = rand:seed(Alg),
    false = (S1 =:= rand:seed_s(Alg)),
    %% Negative integers works
    _ = rand:seed_s(Alg, {-1,-1,-1}),
    %% Check that export_seed/1 returns 'undefined' if there is no seed
    erase(rand_seed),
    undefined = rand:export_seed(),

    %% Other term do not work
    {'EXIT', _} = (catch rand:seed_s(foobar, os:timestamp())),
    {'EXIT', _} = (catch rand:seed_s(Alg, {asd, 1, 1})),
    {'EXIT', _} = (catch rand:seed_s(Alg, {0, 234.1234, 1})),
    {'EXIT', _} = (catch rand:seed_s(Alg, {0, 234, [1, 123, 123]})),
    ok.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check that both APIs are consistent with each other.
api_eq(_Config) ->
    Algs = algs(),
    Small = fun(Alg) ->
		    Seed = rand:seed(Alg),
		    io:format("Seed ~p~n",[rand:export_seed_s(Seed)]),
		    api_eq_1(Seed)
	    end,
    _ = [Small(Alg) || Alg <- Algs],
    ok.

api_eq_1(S00) ->
    Check = fun(_, Seed) ->
		    {V0, S0} = rand:uniform_s(Seed),
		    V0 = rand:uniform(),
		    {V1, S1} = rand:uniform_s(1000000, S0),
		    V1 = rand:uniform(1000000),
		    {V2, S2} = rand:normal_s(S1),
		    V2 = rand:normal(),
		    S2
	    end,
    S1 = lists:foldl(Check, S00, lists:seq(1, 200)),
    S1 = get(rand_seed),
    {V0, S2} = rand:uniform_s(S1),
    V0 = rand:uniform(),
    S2 = get(rand_seed),

    Exported = rand:export_seed(),
    Exported = rand:export_seed_s(S2),

    S3 = lists:foldl(Check, S2, lists:seq(1, 200)),
    S3 = get(rand_seed),

    S4 = lists:foldl(Check, S3, lists:seq(1, 200)),
    S4 = get(rand_seed),
    %% Verify that we do not have loops
    false = S1 =:= S2,
    false = S2 =:= S3,
    false = S3 =:= S4,

    S2 = rand:seed(Exported),
    S3 = lists:foldl(Check, S2, lists:seq(1, 200)),
    ok.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check that uniform/1 returns values within the proper interval.
interval_int(Config) when is_list(Config) ->
    Algs = algs(),
    Small = fun(Alg) ->
		    Seed = rand:seed(Alg),
		    io:format("Seed ~p~n",[rand:export_seed_s(Seed)]),
		    Max = interval_int_1(100000, 7, 0),
		    Max =:= 7 orelse exit({7, Alg, Max})
	    end,
    _ = [Small(Alg) || Alg <- Algs],
    %% Test large integers
    Large = fun(Alg) ->
		    Seed = rand:seed(Alg),
		    io:format("Seed ~p~n",[rand:export_seed_s(Seed)]),
		    Max = interval_int_1(100000, 1 bsl 128, 0),
		    Max > 1 bsl 64 orelse exit({large, Alg, Max})
	    end,
    [Large(Alg) || Alg <- Algs],
    ok.

interval_int_1(0, _, Max) -> Max;
interval_int_1(N, Top, Max) ->
    X = rand:uniform(Top),
    if
	0 < X, X =< Top ->
	    ok;
	true ->
	    io:format("X=~p Top=~p 0<~p<~p~n", [X,Top,X,Top]),
	    exit({X, rand:export_seed()})
    end,
    interval_int_1(N-1, Top, max(X, Max)).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check that uniform/0 returns values within the proper interval.
interval_float(Config) when is_list(Config) ->
    Algs = algs(),
    Test = fun(Alg) ->
		   _ = rand:seed(Alg),
		   interval_float_1(100000)
	   end,
    [Test(Alg) || Alg <- Algs],
    ok.

interval_float_1(0) -> ok;
interval_float_1(N) ->
    X = rand:uniform(),
    Y = rand:uniform_real(),
    if
	0.0 =< X, X < 1.0, 0.0 < Y, Y < 1.0 ->
	    ok;
	true ->
	    io:format("X=~p 0.0=<~p<1.0~n", [X,X]),
	    io:format("Y=~p 0.0<~p<1.0~n", [Y,Y]),
	    exit({X, rand:export_seed()})
    end,
    interval_float_1(N-1).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check if each algorithm generates the proper sequence.
reference(Config) when is_list(Config) ->
    [reference_1(Alg) || Alg <- algs()],
    ok.

reference_1(Alg) ->
    Refval  = reference_val(Alg),
    Testval = gen(Alg),
    case Refval =:= Testval of
        true -> ok;
        false when Refval =:= not_implemented ->
            exit({not_implemented,Alg});
        false ->
	    io:format("Failed: ~p~n",[Alg]),
	    io:format("Length ~p ~p~n",[length(Refval), length(Testval)]),
	    io:format("Head ~p ~p~n",[hd(Refval), hd(Testval)]),
	    exit(wrong_value)
    end.

gen(Algo) ->
    State =
        case Algo of
            exs64 -> %% Printed with orig 'C' code and this seed
                rand:seed_s({exs64, 12345678});
            _ when Algo =:= exsplus; Algo =:= exsp; Algo =:= exrop ->
                %% Printed with orig 'C' code and this seed
                rand:seed_s({Algo, [12345678|12345678]});
            _ when Algo =:= exs1024; Algo =:= exs1024s ->
                %% Printed with orig 'C' code and this seed
                rand:seed_s({Algo, {lists:duplicate(16, 12345678), []}});
            _ ->
                rand:seed(Algo, {100, 200, 300})
        end,
    Max = range(State),
    gen(?LOOP, State, Max, []).

gen(N, State0, Max, Acc) when N > 0 ->
    {Random, State} = rand:uniform_s(Max, State0),
    case N rem (?LOOP div 100) of
	0 -> gen(N-1, State, Max, [Random|Acc]);
	_ -> gen(N-1, State, Max, Acc)
    end;
gen(_, _, _, Acc) -> lists:reverse(Acc).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This just tests the basics so we have not made any serious errors
%% when making the conversion from the original algorithms.
%% The algorithms must have good properties to begin with
%%

%% Check that the algorithms generate sound values.

basic_stats_uniform_1(Config) when is_list(Config) ->
    ct:timetrap({minutes,15}), %% valgrind needs a lot of time
    [basic_uniform_1(?LOOP, rand:seed_s(Alg), 0.0, array:new([{default, 0}]))
     || Alg <- algs()],
    ok.

basic_stats_uniform_2(Config) when is_list(Config) ->
    ct:timetrap({minutes,15}), %% valgrind needs a lot of time
    [basic_uniform_2(?LOOP, rand:seed_s(Alg), 0, array:new([{default, 0}]))
     || Alg <- algs()],
    ok.

basic_stats_standard_normal(Config) when is_list(Config) ->
    ct:timetrap({minutes,6}), %% valgrind needs a lot of time
    io:format("Testing standard normal~n",[]),
    IntendedMean = 0,
    IntendedVariance = 1,
    [basic_normal_1(?LOOP, IntendedMean, IntendedVariance,
                    rand:seed_s(Alg), 0, 0)
     || Alg <- algs()],
    ok.

basic_stats_normal(Config) when is_list(Config) ->
    IntendedMeans = [-1.0e6, -50, -math:pi(), -math:exp(-1),
                     0.12345678, math:exp(1), 100, 1.0e6],
    IntendedVariances = [1.0e-6, math:exp(-1), 1, math:pi(), 1.0e6],
    IntendedMeanVariancePairs =
        [{Mean, Variance} || Mean <- IntendedMeans,
                             Variance <- IntendedVariances],

    ct:timetrap({minutes, 6 * length(IntendedMeanVariancePairs)}), %% valgrind needs a lot of time
    lists:foreach(
      fun ({IntendedMean, IntendedVariance}) ->
              ct:pal(
                "Testing normal(~.2f, ~.2f)~n",
                [float(IntendedMean), float(IntendedVariance)]),
              [basic_normal_1(?LOOP, IntendedMean, IntendedVariance,
                              rand:seed_s(Alg), 0, 0)
               || Alg <- algs()]
      end,
      IntendedMeanVariancePairs).

basic_uniform_1(N, S0, Sum, A0) when N > 0 ->
    {X,S} =
        case N band 1 of
            0 ->
                rand:uniform_s(S0);
            1 ->
                rand:uniform_real_s(S0)
        end,
    I = trunc(X*100),
    A = array:set(I, 1+array:get(I,A0), A0),
    basic_uniform_1(N-1, S, Sum+X, A);
basic_uniform_1(0, {#{type:=Alg}, _}, Sum, A) ->
    AverN = Sum / ?LOOP,
    io:format("~.12w: Average: ~.4f~n", [Alg, AverN]),
    Counters = array:to_list(A),
    Min = lists:min(Counters),
    Max = lists:max(Counters),
    io:format("~.12w: Min: ~p Max: ~p~n", [Alg, Min, Max]),

    %% Verify that the basic statistics are ok
    %% be gentle we don't want to see to many failing tests
    abs(0.5 - AverN) < 0.005 orelse ct:fail({average, Alg, AverN}),
    abs(?LOOP div 100 - Min) < 1000 orelse ct:fail({min, Alg, Min}),
    abs(?LOOP div 100 - Max) < 1000 orelse ct:fail({max, Alg, Max}),
    ok.

basic_uniform_2(N, S0, Sum, A0) when N > 0 ->
    {X,S} = rand:uniform_s(100, S0),
    A = array:set(X-1, 1+array:get(X-1,A0), A0),
    basic_uniform_2(N-1, S, Sum+X, A);
basic_uniform_2(0, {#{type:=Alg}, _}, Sum, A) ->
    AverN = Sum / ?LOOP,
    io:format("~.12w: Average: ~.4f~n", [Alg, AverN]),
    Counters = tl(array:to_list(A)),
    Min = lists:min(Counters),
    Max = lists:max(Counters),
    io:format("~.12w: Min: ~p Max: ~p~n", [Alg, Min, Max]),

    %% Verify that the basic statistics are ok
    %% be gentle we don't want to see to many failing tests
    abs(50.5 - AverN) < 0.5 orelse ct:fail({average, Alg, AverN}),
    abs(?LOOP div 100 - Min) < 1000 orelse ct:fail({min, Alg, Min}),
    abs(?LOOP div 100 - Max) < 1000 orelse ct:fail({max, Alg, Max}),
    ok.

basic_normal_1(N, IntendedMean, IntendedVariance, S0, StandardSum, StandardSq) when N > 0 ->
    {X,S} = normal_s(IntendedMean, IntendedVariance, S0),
    % We now shape X into a standard normal distribution (in case it wasn't already)
    % in order to minimise the accumulated error on Sum / SumSq;
    % otherwise said error would prevent us of making a fair judgment on
    % the overall distribution when targeting large means and variances.
    StandardX = (X - IntendedMean) / math:sqrt(IntendedVariance),
    basic_normal_1(N-1, IntendedMean, IntendedVariance, S,
                   StandardX+StandardSum, StandardX*StandardX+StandardSq);
basic_normal_1(0, _IntendedMean, _IntendedVariance, {#{type:=Alg}, _}, StandardSum, StandardSumSq) ->
    StandardMean = StandardSum / ?LOOP,
    StandardVariance = (StandardSumSq - (StandardSum*StandardSum/?LOOP))/(?LOOP - 1),
    StandardStdDev =  math:sqrt(StandardVariance),
    io:format("~.12w: Standardised Average: ~7.4f, Standardised StdDev ~6.4f~n",
              [Alg, StandardMean, StandardStdDev]),
    %% Verify that the basic statistics are ok
    %% be gentle we don't want to see to many failing tests
    abs(StandardMean) < 0.005 orelse ct:fail({average, Alg, StandardMean}),
    abs(StandardStdDev - 1.0) < 0.005 orelse ct:fail({stddev, Alg, StandardStdDev}),
    ok.

normal_s(Mean, Variance, State0) when Mean == 0, Variance == 1 ->
    % Make sure we're also testing the standard normal interface
    rand:normal_s(State0);
normal_s(Mean, Variance, State0) ->
    rand:normal_s(Mean, Variance, State0).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% White box test of the conversion to float

uniform_real_conv(Config) when is_list(Config) ->
    [begin
%%         ct:pal("~13.16.0bx~3.16.0b: ~p~n", [M,E,Gen]),
         uniform_real_conv_check(M, E, Gen)
     end || {M, E, Gen} <- uniform_real_conv_data()],
    uniform_real_scan(0),
    uniform_real_scan(3).

uniform_real_conv_data() ->
    [{16#fffffffffffff,  -1, [16#3ffffffffffffff]},
     {16#fffffffffffff,  -1, [16#3ffffffffffffe0]},
     {16#ffffffffffffe,  -1, [16#3ffffffffffffdf]},
     %%
     {16#0000000000000,  -1, [16#200000000000000]},
     {16#fffffffffffff,  -2, [16#1ffffffffffffff]},
     {16#fffffffffffff,  -2, [16#1fffffffffffff0]},
     {16#ffffffffffffe,  -2, [16#1ffffffffffffef]},
     %%
     {16#0000000000000,  -2, [16#100000000000000]},
     {16#fffffffffffff,  -3, [16#0ffffffffffffff]},
     {16#fffffffffffff,  -3, [16#0fffffffffffff8]},
     {16#ffffffffffffe,  -3, [16#0fffffffffffff7]},
     %%
     {16#0000000000000,  -3, [16#080000000000000]},
     {16#fffffffffffff,  -4, [16#07fffffffffffff]},
     {16#fffffffffffff,  -4, [16#07ffffffffffffc]},
     {16#ffffffffffffe,  -4, [16#07ffffffffffffb]},
     %%
     {16#0000000000000,  -4, [16#040000000000000]},
     {16#fffffffffffff,  -5, [16#03fffffffffffff,16#3ffffffffffffff]},
     {16#fffffffffffff,  -5, [16#03ffffffffffffe,16#200000000000000]},
     {16#ffffffffffffe,  -5, [16#03fffffffffffff,16#1ffffffffffffff]},
     {16#ffffffffffffe,  -5, [16#03fffffffffffff,16#100000000000000]},
     %%
     {16#0000000000001, -56, [16#000000000000007,16#00000000000007f]},
     {16#0000000000001, -56, [16#000000000000004,16#000000000000040]},
     {16#0000000000000, -57, [16#000000000000003,16#20000000000001f]},
     {16#0000000000000, -57, [16#000000000000000,16#200000000000000]},
     {16#fffffffffffff, -58, [16#000000000000003,16#1ffffffffffffff]},
     {16#fffffffffffff, -58, [16#000000000000000,16#1fffffffffffff0]},
     {16#ffffffffffffe, -58, [16#000000000000000,16#1ffffffffffffef]},
     {16#ffffffffffffe, -58, [16#000000000000000,16#1ffffffffffffe0]},
     %%
     {16#0000000000000, -58, [16#000000000000000,16#10000000000000f]},
     {16#0000000000000, -58, [16#000000000000000,16#100000000000000]},
     {2#11001100000000000000000000000000000000000011000000011, % 53 bits
      -1022,
      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, % 18 zeros
       2#1100110000000000000000000000000000000000001 bsl 2, % 43 bits
       2#1000000011 bsl (56-10+2)]}, % 10 bits
     {0, -1, % 0.5 after retry
      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, % 18 zeros
       2#111111111111111111111111111111111111111111 bsl 2, % 42 bits - retry
       16#200000000000003]}]. % 0.5

-define(UNIFORM_REAL_SCAN_PATTERN, (16#19000000000009)). % 53 bits
-define(UNIFORM_REAL_SCAN_NUMBER, (1021)).

uniform_real_scan_template(K) ->
    <<0:?UNIFORM_REAL_SCAN_NUMBER,
      ?UNIFORM_REAL_SCAN_PATTERN:53,K:2,0:1>>.

uniform_real_scan(K) ->
    Templ = uniform_real_scan_template(K),
    N = ?UNIFORM_REAL_SCAN_NUMBER,
    uniform_real_scan(Templ, N, K).

uniform_real_scan(Templ, N, K) when 0 =< N ->
    <<_:N/bits,T/bits>> = Templ,
    Data = uniform_real_scan_data(T, K),
    uniform_real_conv_check(
      ?UNIFORM_REAL_SCAN_PATTERN, N - 1 - ?UNIFORM_REAL_SCAN_NUMBER, Data),
    uniform_real_scan(Templ, N - 1, K);
uniform_real_scan(_, _, _) ->
    ok.

uniform_real_scan_data(Templ, K) ->
    case Templ of
        <<X:56, T/bits>> ->
            B = rand:bc64(X),
            [(X bsl 2) bor K |
             if
                 53 =< B ->
                     [];
                 true ->
                     uniform_real_scan_data(T, K)
             end];
        _ ->
            <<X:56, _/bits>> = <<Templ/bits, 0:56>>,
            [(X bsl 2) bor K]
    end.

uniform_real_conv_check(M, E, Gen) ->
    <<F/float>> = <<0:1, (E + 16#3ff):11, M:52>>,
    try uniform_real_gen(Gen) of
        F -> F;
        FF ->
            ct:pal(
              "~s =/= ~s: ~s~n",
              [rand:float2str(FF), rand:float2str(F),
               [["16#",integer_to_list(G,16),$\s]||G<-Gen]]),
            ct:fail({neq, FF, F})
    catch
        Error:Reason ->
            ct:pal(
              "~w:~p ~s: ~s~n",
              [Error, Reason, rand:float2str(F),
               [["16#",integer_to_list(G,16),$\s]||G<-Gen]]),
            ct:fail({Error, Reason, F, erlang:get_stacktrace()})
    end.


uniform_real_gen(Gen) ->
    State = rand_state(Gen),
    {F, {#{type := rand_SUITE_list},[]}} = rand:uniform_real_s(State),
    F.

uniform_gen(Range, Gen) ->
    State = rand_state(Gen),
    {N, {#{type := rand_SUITE_list},[]}} = rand:uniform_s(Range, State),
    N.

%% Loaded dice for white box tests
rand_state(Gen) ->
    {#{type => rand_SUITE_list, bits => 58, weak_low_bits => 1,
       next => fun ([H|T]) -> {H, T} end},
     Gen}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Test that the user can write algorithms.
plugin(Config) when is_list(Config) ->
    try crypto:strong_rand_bytes(1) of
        <<_>> ->
            _ = lists:foldl(
                  fun(_, S0) ->
                          {V1, S1} = rand:uniform_s(10000, S0),
                          true = is_integer(V1),
                          {V2, S2} = rand:uniform_s(S1),
                          true = is_float(V2),
                          S2
                  end, crypto64_seed(), lists:seq(1, 200)),
            ok
    catch
        error:low_entropy ->
            {skip,low_entropy};
        error:undef ->
            {skip,no_crypto}
    end.

%% Test implementation
crypto64_seed() ->
    {#{type=>crypto64,
       bits=>64,
       next=>fun crypto64_next/1,
       uniform=>fun crypto64_uniform/1,
       uniform_n=>fun crypto64_uniform_n/2},
     <<>>}.

%% Be fair and create bignums i.e. 64bits otherwise use 58bits
crypto64_next(<<Num:64, Bin/binary>>) ->
    {Num, Bin};
crypto64_next(_) ->
    crypto64_next(crypto:strong_rand_bytes((64 div 8)*100)).

crypto64_uniform({Api, Data0}) ->
    {Int, Data} = crypto64_next(Data0),
    {Int / (1 bsl 64), {Api, Data}}.

crypto64_uniform_n(N, {Api, Data0}) when N < (1 bsl 64) ->
    {Int, Data} = crypto64_next(Data0),
    {(Int rem N)+1, {Api, Data}};
crypto64_uniform_n(N, State0) ->
    {F,State} = crypto64_uniform(State0),
    {trunc(F * N) + 1, State}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Not a test but measures the time characteristics of the different algorithms
measure(Config) ->
    ct:timetrap({minutes,60}), %% valgrind needs a lot of time
    case ct:get_timetrap_info() of
        {_,{_,1}} -> % No scaling
            do_measure(Config);
        {_,{_,Scale}} ->
            {skip,{will_not_run_in_scaled_time,Scale}}
    end.

-define(CHECK_UNIFORM_RANGE(Gen, Range, X, St),
        case (Gen) of
            {(X), (St)} when is_integer(X), 1 =< (X), (X) =< (Range) ->
                St
        end).
-define(CHECK_UNIFORM(Gen, X, St),
        case (Gen) of
            {(X), (St)} when is_float(X), 0.0 =< (X), (X) < 1.0 ->
                St
        end).
-define(CHECK_UNIFORM_NZ(Gen, X, St),
        case (Gen) of
            {(X), (St)} when is_float(X), 0.0 < (X), (X) =< 1.0 ->
                St
        end).
-define(CHECK_NORMAL(Gen, X, St),
        case (Gen) of
            {(X), (St)} when is_float(X) ->
                St
        end).

do_measure(_Config) ->
    Algs =
        algs() ++
        try crypto:strong_rand_bytes(1) of
            <<_>> -> [crypto64, crypto_cache, crypto]
        catch
            error:low_entropy -> [];
            error:undef -> []
        end,
    %%
    ct:pal("~nRNG uniform integer range 10000 performance~n",[]),
    _ =
        measure_1(
          fun (_) -> 10000 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer 32 bit performance~n",[]),
    _ =
        measure_1(
          fun (_) -> 1 bsl 32 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer half range performance~n",[]),
    _ =
        measure_1(
          fun (State) -> half_range(State) end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer half range + 1 performance~n",[]),
    _ =
        measure_1(
          fun (State) -> half_range(State) + 1 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer full range - 1 performance~n",[]),
    _ =
        measure_1(
          fun (State) -> (half_range(State) bsl 1) - 1 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer full range performance~n",[]),
    _ =
        measure_1(
          fun (State) -> half_range(State) bsl 1 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer full range + 1 performance~n",[]),
    _ =
        measure_1(
          fun (State) -> (half_range(State) bsl 1) + 1 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer double range performance~n",[]),
    _ =
        measure_1(
          fun (State) ->
                  half_range(State) bsl 2
          end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer double range + 1  performance~n",[]),
    _ =
        measure_1(
          fun (State) ->
                  (half_range(State) bsl 2) + 1
          end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform integer 64 bit performance~n",[]),
    _ =
        measure_1(
          fun (_) -> 1 bsl 64 end,
          fun (State, Range, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM_RANGE(
                               Mod:uniform_s(Range, St0), Range,
                               X, St1)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform float performance~n",[]),
    _ =
        measure_1(
          fun (_) -> 0 end,
          fun (State, _, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM(Mod:uniform_s(St0), X, St)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG uniform_real float performance~n",[]),
    _ =
        measure_1(
          fun (_) -> 0 end,
          fun (State, _, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_UNIFORM(Mod:uniform_real_s(St0), X, St)
                    end,
                    State)
          end,
          Algs),
    %%
    ct:pal("~nRNG normal float performance~n",[]),
    [TMarkNormalFloat|_] =
        measure_1(
          fun (_) -> 0 end,
          fun (State, _, Mod) ->
                  measure_loop(
                    fun (St0) ->
                            ?CHECK_NORMAL(Mod:normal_s(St0), X, St1)
                    end,
                    State)
          end,
          Algs),
    %% Just for fun try an implementation of the Box-Muller
    %% transformation for creating normal distribution floats
    %% to compare with our Ziggurat implementation.
    %% Generates two numbers per call that we add so they
    %% will not be optimized away.  Hence the benchmark time
    %% is twice what it should be.
    TwoPi = 2 * math:pi(),
    _ =
        measure_1(
          fun (_) -> 0 end,
          fun (State, _, Mod) ->
                  measure_loop(
                    fun (State0) ->
                            {U1, State1} = Mod:uniform_real_s(State0),
                            {U2, State2} = Mod:uniform_s(State1),
                            R = math:sqrt(-2.0 * math:log(U1)),
                            T = TwoPi * U2,
                            Z0 = R * math:cos(T),
                            Z1 = R * math:sin(T),
                            ?CHECK_NORMAL({Z0 + Z1, State2}, X, State3)
                    end,
                    State)
          end,
          exrop, TMarkNormalFloat),
    ok.

-define(LOOP_MEASURE, (?LOOP div 5)).

measure_loop(Fun, State) ->
    measure_loop(Fun, State, ?LOOP_MEASURE).
%%
measure_loop(Fun, State, N) when 0 < N ->
    measure_loop(Fun, Fun(State), N-1);
measure_loop(_, _, _) ->
    ok.

measure_1(RangeFun, Fun, Algs) ->
    TMark = measure_1(RangeFun, Fun, hd(Algs), undefined),
    [TMark] ++
        [measure_1(RangeFun, Fun, Alg, TMark) || Alg <- tl(Algs)].

measure_1(RangeFun, Fun, Alg, TMark) ->
    Parent = self(),
    {Mod, State} =
        case Alg of
            crypto64 ->
                {rand, crypto64_seed()};
            crypto_cache ->
                {rand, crypto:rand_seed_alg(crypto_cache)};
            crypto ->
                {rand, crypto:rand_seed_s()};
            random ->
                {random, random:seed(os:timestamp()), get(random_seed)};
            _ ->
                {rand, rand:seed_s(Alg)}
        end,
    Range = RangeFun(State),
    Pid = spawn_link(
            fun() ->
                    {Time, ok} = timer:tc(fun () -> Fun(State, Range, Mod) end),
                    Percent =
                        case TMark of
                            undefined -> 100;
                            _ -> (Time * 100 + 50) div TMark
                        end,
                    io:format(
                      "~.12w: ~p ns ~p% [16#~.16b]~n",
                      [Alg, (Time * 1000 + 500) div ?LOOP_MEASURE,
                       Percent, Range]),
                    Parent ! {self(), Time},
                    normal
            end),
    receive
	{Pid, Msg} -> Msg
    end.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The jump sequence tests has two parts
%% for those with the functional API (jump/1)
%% and for those with the internal state
%% in process dictionary (jump/0).

-define(LOOP_JUMP, (?LOOP div 1000)).

%% Check if each algorithm generates the proper jump sequence
%% with the functional API.
reference_jump_state(Config) when is_list(Config) ->
    [reference_jump_1(Alg) || Alg <- algs()],
    ok.

reference_jump_1(Alg) ->
    Refval  = reference_jump_val(Alg),
    Testval = gen_jump_1(Alg),
    case Refval =:= Testval of
        true -> ok;
        false ->
{'fraction_non_alphanumeric': 0.16160696873171684, 'fraction_numerical': 0.056052135474224794, 'mean_word_length': 1.9896511514915947, 'pattern_counts': {'":': 0, '<': 141, '<?xml version=': 0, '>': 397, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 162, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6790493', 'n_tokens_mistral': 23627, 'n_tokens_neox': 20670, 'n_words': 4677}