From 437555fd6c495915773b0f9ade7aad3fd0a73a1b Mon Sep 17 00:00:00 2001 From: Raimo Niskanen Date: Tue, 21 Mar 2017 16:36:33 +0100 Subject: Implement Xoroshiro116+ and improve statisticals Implement Xoroshiro116+ as 'exrop' with fixes. Deprecate all old algorithms but reincarnate 'exs1024' as 'exs1024s' and 'exsplus' as 'exsp' with fixes. Fixes: * Avoid skew for uniform integers caused by using a simple 'rem' operation for range confinement. Correctness requires retry with new random value for an unfortunate first value. * Implement a correct algorithm that collects enough random bits for ranges larger than the generator's precision. * Fix uniform density for floats by acquiring 53 bits then multiplying with 2.0^(-53) which produces floats on the form N * 2.0^(-53). --- lib/stdlib/test/rand_SUITE.erl | 544 ++++++++++++++++++++++++++++++++--------- 1 file changed, 425 insertions(+), 119 deletions(-) (limited to 'lib/stdlib/test/rand_SUITE.erl') diff --git a/lib/stdlib/test/rand_SUITE.erl b/lib/stdlib/test/rand_SUITE.erl index 098eefeb61..51bb03f572 100644 --- a/lib/stdlib/test/rand_SUITE.erl +++ b/lib/stdlib/test/rand_SUITE.erl @@ -1,7 +1,7 @@ %% %% %CopyrightBegin% %% -%% Copyright Ericsson AB 2000-2016. All Rights Reserved. +%% 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. @@ -66,18 +66,19 @@ group(reference_jump) -> %% 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). + 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() -> - [exs64, exsplus, exs1024]. + [exs64, exsplus, exsp, exrop, exs1024, exs1024s]. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -226,10 +227,10 @@ interval_float_1(0) -> ok; interval_float_1(N) -> X = rand:uniform(), if - 0.0 < X, X < 1.0 -> + 0.0 =< X, X < 1.0 -> ok; true -> - io:format("X=~p 0<~p<1.0~n", [X,X]), + io:format("X=~p 0=<~p<1.0~n", [X,X]), exit({X, rand:export_seed()}) end, interval_float_1(N-1). @@ -246,6 +247,8 @@ reference_1(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)]), @@ -254,25 +257,29 @@ reference_1(Alg) -> end. gen(Algo) -> - Seed = case Algo of - exsplus -> %% Printed with orig 'C' code and this seed - rand:seed_s({exsplus, [12345678|12345678]}); - exs64 -> %% Printed with orig 'C' code and this seed - rand:seed_s({exs64, 12345678}); - exs1024 -> %% Printed with orig 'C' code and this seed - rand:seed_s({exs1024, {lists:duplicate(16, 12345678), []}}); - _ -> - rand:seed(Algo, {100, 200, 300}) - end, - gen(?LOOP, Seed, []). - -gen(N, State0 = {#{max:=Max}, _}, Acc) when N > 0 -> + 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, [Random|Acc]); - _ -> gen(N-1, State, Acc) + 0 -> gen(N-1, State, Max, [Random|Acc]); + _ -> gen(N-1, State, Max, Acc) end; -gen(_, _, Acc) -> lists:reverse(Acc). +gen(_, _, _, Acc) -> lists:reverse(Acc). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% This just tests the basics so we have not made any serious errors @@ -307,11 +314,11 @@ basic_uniform_1(N, S0, Sum, A0) when N > 0 -> basic_uniform_1(N-1, S, Sum+X, A); basic_uniform_1(0, {#{type:=Alg}, _}, Sum, A) -> AverN = Sum / ?LOOP, - io:format("~.10w: Average: ~.4f~n", [Alg, AverN]), + io:format("~.12w: Average: ~.4f~n", [Alg, AverN]), Counters = array:to_list(A), Min = lists:min(Counters), Max = lists:max(Counters), - io:format("~.10w: Min: ~p Max: ~p~n", [Alg, Min, Max]), + 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 @@ -326,11 +333,11 @@ basic_uniform_2(N, S0, Sum, A0) when N > 0 -> basic_uniform_2(N-1, S, Sum+X, A); basic_uniform_2(0, {#{type:=Alg}, _}, Sum, A) -> AverN = Sum / ?LOOP, - io:format("~.10w: Average: ~.4f~n", [Alg, AverN]), + io:format("~.12w: Average: ~.4f~n", [Alg, AverN]), Counters = tl(array:to_list(A)), Min = lists:min(Counters), Max = lists:max(Counters), - io:format("~.10w: Min: ~p Max: ~p~n", [Alg, Min, Max]), + 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 @@ -345,7 +352,7 @@ basic_normal_1(N, S0, Sum, Sq) when N > 0 -> basic_normal_1(0, {#{type:=Alg}, _}, Sum, SumSq) -> Mean = Sum / ?LOOP, StdDev = math:sqrt((SumSq - (Sum*Sum/?LOOP))/(?LOOP - 1)), - io:format("~.10w: Average: ~7.4f StdDev ~6.4f~n", [Alg, Mean, StdDev]), + io:format("~.12w: Average: ~7.4f StdDev ~6.4f~n", [Alg, Mean, StdDev]), %% Verify that the basic statistics are ok %% be gentle we don't want to see to many failing tests abs(Mean) < 0.005 orelse ct:fail({average, Alg, Mean}), @@ -365,7 +372,7 @@ plugin(Config) when is_list(Config) -> {V2, S2} = rand:uniform_s(S1), true = is_float(V2), S2 - end, crypto_seed(), lists:seq(1, 200)), + end, crypto64_seed(), lists:seq(1, 200)), ok catch error:low_entropy -> @@ -375,86 +382,220 @@ plugin(Config) when is_list(Config) -> end. %% Test implementation -crypto_seed() -> - {#{type=>crypto, - max=>(1 bsl 64)-1, - next=>fun crypto_next/1, - uniform=>fun crypto_uniform/1, - uniform_n=>fun crypto_uniform_n/2}, +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 -crypto_next(<>) -> +crypto64_next(<>) -> {Num, Bin}; -crypto_next(_) -> - crypto_next(crypto:strong_rand_bytes((64 div 8)*100)). +crypto64_next(_) -> + crypto64_next(crypto:strong_rand_bytes((64 div 8)*100)). -crypto_uniform({Api, Data0}) -> - {Int, Data} = crypto_next(Data0), +crypto64_uniform({Api, Data0}) -> + {Int, Data} = crypto64_next(Data0), {Int / (1 bsl 64), {Api, Data}}. -crypto_uniform_n(N, {Api, Data0}) when N < (1 bsl 64) -> - {Int, Data} = crypto_next(Data0), +crypto64_uniform_n(N, {Api, Data0}) when N < (1 bsl 64) -> + {Int, Data} = crypto64_next(Data0), {(Int rem N)+1, {Api, Data}}; -crypto_uniform_n(N, State0) -> - {F,State} = crypto_uniform(State0), +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(Suite) when is_atom(Suite) -> []; -measure(_Config) -> - ct:timetrap({minutes,15}), %% valgrind needs a lot of time +measure(Config) -> + ct:timetrap({minutes,30}), %% 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. + +do_measure(_Config) -> Algos = try crypto:strong_rand_bytes(1) of - <<_>> -> [crypto64] + <<_>> -> [crypto64, crypto] catch error:low_entropy -> []; error:undef -> [] end ++ algs(), - io:format("RNG uniform integer performance~n",[]), - _ = measure_1(random, fun(State) -> {int, random:uniform_s(10000, State)} end), - _ = [measure_1(Algo, fun(State) -> {int, rand:uniform_s(10000, State)} end) || Algo <- Algos], - io:format("RNG uniform float performance~n",[]), - _ = measure_1(random, fun(State) -> {uniform, random:uniform_s(State)} end), - _ = [measure_1(Algo, fun(State) -> {uniform, rand:uniform_s(State)} end) || Algo <- Algos], - io:format("RNG normal float performance~n",[]), - io:format("~.10w: not implemented (too few bits)~n", [random]), - _ = [measure_1(Algo, fun(State) -> {normal, rand:normal_s(State)} end) || Algo <- Algos], + %% + ct:pal("RNG uniform integer performance~n",[]), + TMark1 = + measure_1( + random, + fun (_) -> 10000 end, + undefined, + fun (Range, State) -> + {int, random:uniform_s(Range, State)} + end), + _ = + [measure_1( + Algo, + fun (_) -> 10000 end, + TMark1, + fun (Range, State) -> + {int, rand:uniform_s(Range, State)} + end) || Algo <- Algos], + %% + ct:pal("~nRNG uniform integer 2^(N-1) performance~n",[]), + RangeTwoPowFun = fun (State) -> quart_range(State) bsl 1 end, + TMark2 = + measure_1( + random, + RangeTwoPowFun, + undefined, + fun (Range, State) -> + {int, random:uniform_s(Range, State)} + end), + _ = + [measure_1( + Algo, + RangeTwoPowFun, + TMark2, + fun (Range, State) -> + {int, rand:uniform_s(Range, State)} + end) || Algo <- Algos], + %% + ct:pal("~nRNG uniform integer 3*2^(N-2)+1 performance~n",[]), + RangeLargeFun = fun (State) -> 3 * quart_range(State) + 1 end, + TMark3 = + measure_1( + random, + RangeLargeFun, + undefined, + fun (Range, State) -> + {int, random:uniform_s(Range, State)} + end), + _ = + [measure_1( + Algo, + RangeLargeFun, + TMark3, + fun (Range, State) -> + {int, rand:uniform_s(Range, State)} + end) || Algo <- Algos], + %% + ct:pal("~nRNG uniform integer 2^128 performance~n",[]), + TMark4 = + measure_1( + random, + fun (_) -> 1 bsl 128 end, + undefined, + fun (Range, State) -> + {int, random:uniform_s(Range, State)} + end), + _ = + [measure_1( + Algo, + fun (_) -> 1 bsl 128 end, + TMark4, + fun (Range, State) -> + {int, rand:uniform_s(Range, State)} + end) || Algo <- Algos], + %% + ct:pal("~nRNG uniform integer 2^128 + 1 performance~n",[]), + TMark5 = + measure_1( + random, + fun (_) -> (1 bsl 128) + 1 end, + undefined, + fun (Range, State) -> + {int, random:uniform_s(Range, State)} + end), + _ = + [measure_1( + Algo, + fun (_) -> (1 bsl 128) + 1 end, + TMark5, + fun (Range, State) -> + {int, rand:uniform_s(Range, State)} + end) || Algo <- Algos], + %% + ct:pal("~nRNG uniform float performance~n",[]), + TMark6 = + measure_1( + random, + fun (_) -> 0 end, + undefined, + fun (_, State) -> + {uniform, random:uniform_s(State)} + end), + _ = + [measure_1( + Algo, + fun (_) -> 0 end, + TMark6, + fun (_, State) -> + {uniform, rand:uniform_s(State)} + end) || Algo <- Algos], + %% + ct:pal("~nRNG normal float performance~n",[]), + io:format("~.12w: not implemented (too few bits)~n", [random]), + _ = [measure_1( + Algo, + fun (_) -> 0 end, + TMark6, + fun (_, State) -> + {normal, rand:normal_s(State)} + end) || Algo <- Algos], ok. -measure_1(Algo, Gen) -> +measure_1(Algo, RangeFun, TMark, Gen) -> Parent = self(), - Seed = fun(crypto64) -> crypto_seed(); - (random) -> random:seed(os:timestamp()), get(random_seed); - (Alg) -> rand:seed_s(Alg) - end, - - Pid = spawn_link(fun() -> - Fun = fun() -> measure_2(?LOOP, Seed(Algo), Gen) end, - {Time, ok} = timer:tc(Fun), - io:format("~.10w: ~pµs~n", [Algo, Time]), - Parent ! {self(), ok}, - normal - end), + Seed = + case Algo of + crypto64 -> + crypto64_seed(); + crypto -> + crypto:rand_seed_s(); + random -> + random:seed(os:timestamp()), get(random_seed); + _ -> + rand:seed_s(Algo) + end, + Range = RangeFun(Seed), + Pid = spawn_link( + fun() -> + Fun = fun() -> measure_2(?LOOP, Range, Seed, Gen) end, + {Time, ok} = timer:tc(Fun), + Percent = + case TMark of + undefined -> 100; + _ -> (Time * 100 + 50) div TMark + end, + io:format( + "~.12w: ~p ns ~p% [16#~.16b]~n", + [Algo, (Time * 1000 + 500) div ?LOOP, Percent, Range]), + Parent ! {self(), Time}, + normal + end), receive {Pid, Msg} -> Msg end. -measure_2(N, State0, Fun) when N > 0 -> - case Fun(State0) of +measure_2(N, Range, State0, Fun) when N > 0 -> + case Fun(Range, State0) of {int, {Random, State}} - when is_integer(Random), Random >= 1, Random =< 100000 -> - measure_2(N-1, State, Fun); - {uniform, {Random, State}} when is_float(Random), Random > 0, Random < 1 -> - measure_2(N-1, State, Fun); + when is_integer(Random), Random >= 1, Random =< Range -> + measure_2(N-1, Range, State, Fun); + {uniform, {Random, State}} + when is_float(Random), 0.0 =< Random, Random < 1.0 -> + measure_2(N-1, Range, State, Fun); {normal, {Random, State}} when is_float(Random) -> - measure_2(N-1, State, Fun); + measure_2(N-1, Range, State, Fun); Res -> exit({error, Res, State0}) end; -measure_2(0, _, _) -> ok. +measure_2(0, _, _, _) -> ok. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% The jump sequence tests has two parts @@ -479,36 +620,43 @@ reference_jump_1(Alg) -> 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)]), + io:format("Vals ~p ~p~n",[Refval, Testval]), exit(wrong_value) end. gen_jump_1(Algo) -> - Seed = case Algo of - exsplus -> %% Printed with orig 'C' code and this seed - rand:seed_s({exsplus, [12345678|12345678]}); - exs1024 -> %% Printed with orig 'C' code and this seed - rand:seed_s({exs1024, {lists:duplicate(16, 12345678), []}}); - exs64 -> %% Test exception of not_implemented notice - try rand:jump(rand:seed_s(exs64)) - catch - error:not_implemented -> not_implemented - end; - _ -> % unimplemented - not_implemented - end, - case Seed of + State = + case Algo of + exs64 -> %% Test exception of not_implemented notice + try rand:jump(rand:seed_s(exs64)) + catch + error:not_implemented -> not_implemented + end; + _ 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), []}}); + _ -> % unimplemented + not_implemented + end, + case State of not_implemented -> [not_implemented]; - S -> gen_jump_1(?LOOP_JUMP, S, []) + _ -> + Max = range(State), + gen_jump_1(?LOOP_JUMP, State, Max, []) end. -gen_jump_1(N, State0 = {#{max:=Max}, _}, Acc) when N > 0 -> +gen_jump_1(N, State0, Max, Acc) when N > 0 -> {_, State1} = rand:uniform_s(Max, State0), {Random, State2} = rand:uniform_s(Max, rand:jump(State1)), case N rem (?LOOP_JUMP div 100) of - 0 -> gen_jump_1(N-1, State2, [Random|Acc]); - _ -> gen_jump_1(N-1, State2, Acc) + 0 -> gen_jump_1(N-1, State2, Max, [Random|Acc]); + _ -> gen_jump_1(N-1, State2, Max, Acc) end; -gen_jump_1(_, _, Acc) -> lists:reverse(Acc). +gen_jump_1(_, _, _, Acc) -> lists:reverse(Acc). + %% Check if each algorithm generates the proper jump sequence %% with the internal state in the process dictionary. @@ -530,25 +678,26 @@ reference_jump_0(Alg) -> gen_jump_0(Algo) -> Seed = case Algo of - exsplus -> %% Printed with orig 'C' code and this seed - rand:seed({exsplus, [12345678|12345678]}); - exs1024 -> %% Printed with orig 'C' code and this seed - rand:seed({exs1024, {lists:duplicate(16, 12345678), []}}); exs64 -> %% Test exception of not_implemented notice - try - _ = rand:seed(exs64), - rand:jump() - catch - error:not_implemented -> not_implemented - end; + try + _ = rand:seed(exs64), + rand:jump() + catch + error:not_implemented -> not_implemented + end; + _ when Algo =:= exsplus; Algo =:= exsp; Algo =:= exrop -> + %% Printed with orig 'C' code and this seed + rand:seed({Algo, [12345678|12345678]}); + _ when Algo =:= exs1024; Algo =:= exs1024s -> + %% Printed with orig 'C' code and this seed + rand:seed({Algo, {lists:duplicate(16, 12345678), []}}); _ -> % unimplemented not_implemented end, case Seed of not_implemented -> [not_implemented]; - S -> - {Seedmap=#{}, _} = S, - Max = maps:get(max, Seedmap), + _ -> + Max = range(Seed), gen_jump_0(?LOOP_JUMP, Max, []) end. @@ -643,9 +792,77 @@ reference_val(exsplus) -> 16#6c6145ffa1169d,16#18ec2c393d45359,16#1f1a5f256e7130c,16#131cc2f49b8004f, 16#36f715a249f4ec2,16#1c27629826c50d3,16#914d9a6648726a,16#27f5bf5ce2301e8, 16#3dd493b8012970f,16#be13bed1e00e5c,16#ceef033b74ae10,16#3da38c6a50abe03, - 16#15cbd1a421c7a8c,16#22794e3ec6ef3b1,16#26154d26e7ea99f,16#3a66681359a6ab6]. + 16#15cbd1a421c7a8c,16#22794e3ec6ef3b1,16#26154d26e7ea99f,16#3a66681359a6ab6]; + +reference_val(exsp) -> + reference_val(exsplus); +reference_val(exs1024s) -> + reference_val(exs1024); +reference_val(exrop) -> +%% #include +%% #include +%% +%% uint64_t s[2]; +%% uint64_t next(void); +%% /* Xoroshiro116+ PRNG here */ +%% +%% int main(char *argv[]) { +%% int n; +%% uint64_t r; +%% s[0] = 12345678; +%% s[1] = 12345678; +%% +%% for (n = 1000000; n > 0; n--) { +%% r = next(); +%% if ((n % 10000) == 0) { +%% printf("%llu,", (unsigned long long) (r + 1)); +%% } +%% } +%% printf("\n"); +%% } + [24691357,29089185972758626,135434857127264790, + 277209758236304485,101045429972817342, + 241950202080388093,283018380268425711,268233672110762489, + 173241488791227202,245038518481669421, + 253627577363613736,234979870724373477,115607127954560275, + 96445882796968228,166106849348423677, + 83614184550774836,109634510785746957,68415533259662436, + 12078288820568786,246413981014863011, + 96953486962147513,138629231038332640,206078430370986460, + 11002780552565714,238837272913629203, + 60272901610411077,148828243883348685,203140738399788939, + 131001610760610046,30717739120305678, + 262903815608472425,31891125663924935,107252017522511256, + 241577109487224033,263801934853180827, + 155517416581881714,223609336630639997,112175917931581716, + 16523497284706825,201453767973653420, + 35912153101632769,211525452750005043,96678037860996922, + 70962216125870068,107383886372877124, + 223441708670831233,247351119445661499,233235283318278995, + 280646255087307741,232948506631162445, + %% + 117394974124526779,55395923845250321,274512622756597759, + 31754154862553492,222645458401498438, + 161643932692872858,11771755227312868,93933211280589745, + 92242631276348831,197206910466548143, + 150370169849735808,229903773212075765,264650708561842793, + 30318996509793571,158249985447105184, + 220423733894955738,62892844479829080,112941952955911674, + 203157000073363030,54175707830615686, + 50121351829191185,115891831802446962,62298417197154985, + 6569598473421167,69822368618978464, + 176271134892968134,160793729023716344,271997399244980560, + 59100661824817999,150500611720118722, + 23707133151561128,25156834940231911,257788052162304719, + 176517852966055005,247173855600850875, + 83440973524473396,94711136045581604,154881198769946042, + 236537934330658377,152283781345006019, + 250789092615679985,78848633178610658,72059442721196128, + 98223942961505519,191144652663779840, + 102425686803727694,89058927716079076,80721467542933080, + 8462479817391645,2774921106204163]. -%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% reference_jump_val(exsplus) -> [82445318862816932, 145810727464480743, 16514517716894509, 247642377064868650, @@ -701,4 +918,93 @@ reference_jump_val(exs1024) -> 17936751184378118743, 4224632875737239207, 15888641556987476199, 9586888813112229805, 9476861567287505094, 14909536929239540332, 17996844556292992842, 2699310519182298856]; -reference_jump_val(exs64) -> [not_implemented]. +reference_jump_val(exsp) -> + reference_jump_val(exsplus); +reference_jump_val(exs1024s) -> + reference_jump_val(exs1024); +reference_jump_val(exs64) -> [not_implemented]; +reference_jump_val(exrop) -> +%% #include +%% #include +%% +%% uint64_t s[2]; +%% uint64_t next(void); +%% /* Xoroshiro116+ PRNG here */ +%% +%% int main(char *argv[]) { +%% int n; +%% uint64_t r; +%% s[0] = 12345678; +%% s[1] = 12345678; + +%% for (n = 1000; n > 0; n--) { +%% next(); +%% jump(); +%% r = next(); +%% if ((n % 10) == 0) { +%% printf("%llu,", (unsigned long long) (r + 1)); +%% } +%% } +%% printf("\n"); +%% } + [60301713907476001,135397949584721850,4148159712710727, + 110297784509908316,18753463199438866, + 106699913259182846,2414728156662676,237591345910610406, + 48519427605486503,38071665570452612, + 235484041375354592,45428997361037927,112352324717959775, + 226084403445232507,270797890380258829, + 160587966336947922,80453153271416820,222758573634013699, + 195715386237881435,240975253876429810, + 93387593470886224,23845439014202236,235376123357642262, + 22286175195310374,239068556844083490, + 120126027410954482,250690865061862527,113265144383673111, + 57986825640269127,206087920253971490, + 265971029949338955,40654558754415167,185972161822891882, + 72224917962819036,116613804322063968, + 129103518989198416,236110607653724474,98446977363728314, + 122264213760984600,55635665885245081, + 42625530794327559,288031254029912894,81654312180555835, + 261800844953573559,144734008151358432, + 77095621402920587,286730580569820386,274596992060316466, + 97977034409404188,5517946553518132, + %% + 56460292644964432,252118572460428657,38694442746260303, + 165653145330192194,136968555571402812, + 64905200201714082,257386366768713186,22702362175273017, + 208480936480037395,152926769756967697, + 256751159334239189,130982960476845557,21613531985982870, + 87016962652282927,130446710536726404, + 188769410109327420,282891129440391928,251807515151187951, + 262029034126352975,30694713572208714, + 46430187445005589,176983177204884508,144190360369444480, + 14245137612606100,126045457407279122, + 169277107135012393,42599413368851184,130940158341360014, + 113412693367677211,119353175256553456, + 96339829771832349,17378172025472134,110141940813943768, + 253735613682893347,234964721082540068, + 85668779779185140,164542570671430062,18205512302089755, + 282380693509970845,190996054681051049, + 250227633882474729,171181147785250210,55437891969696407, + 241227318715885854,77323084015890802, + 1663590009695191,234064400749487599,222983191707424780, + 254956809144783896,203898972156838252]. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% The old algorithms used a range 2^N - 1 for their reference val +%% tests, which was incorrect but works as long as you do not draw +%% the value 2^N, which is very unlikely. It was not possible +%% to simply correct the range to 2^N due to another incorrectness +%% in that the old algorithms changed to using the broken +%% (multiply a float approach with too few bits) approach for +%% ranges >= 2^N. This function digs out the range to use +%% for the reference tests for old and new algorithms. +range({#{bits:=Bits}, _}) -> 1 bsl Bits; +range({#{max:=Max}, _}) -> Max; %% Old incorrect range +range({_, _, _}) -> 51. % random + + +quart_range({#{bits:=Bits}, _}) -> 1 bsl (Bits - 2); +quart_range({#{max:=Max}, _}) -> (Max bsr 2) + 1; +quart_range({#{}, _}) -> 1 bsl 62; % crypto +quart_range({_, _, _}) -> 1 bsl 49. % random -- cgit v1.2.3