1:- module(r_random,
    2[
    3    rg/1,
    4    rg/3,
    5    random_generator/1,
    6    random_generator/3,
    7    random/3,
    8    grandom//1,
    9    random_bounded/4,
   10    grandom_bounded//2,
   11    random_between/5,
   12    grandom_between//3,
   13    random_list/4,
   14    grandom_list//2,
   15    random_boundedlist/5,
   16    grandom_boundedlist//3,
   17    random_betweenlist/6,
   18    grandom_betweenlist//4,
   19    probability/4,
   20    gprobability//2,
   21    random_permutation/4,
   22    grandom_permutation//2
   23]).   24
   25:- use_module(library(reif)).   26:- use_module(library(clpfd)).   27:- set_prolog_flag(clpfd_monotonic, true).   28
   29
   30% reference Minimal C Implementation: https://www.pcg-random.org/download.html
   31% PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation
   32% Author: "Melissa E. O'Neill"
   33% based on https://github.com/imneme/pcg-c-basic/blob/master/pcg_basic.c
   34
   35% typedef struct { uint64_t state;  uint64_t inc; } pcg32_random_t;
   36% pcg32randomt(State,Inc).
   37
   38% -------------------- public predicates ---------------------------:
   39
   40% A default random number generator.
   41random_generator( pcg32randomt(0x853c49e6748fea9b, 0xda3e39cb94b95bdb) ).
   42
   43random_generator( InitState, InitSeq, Pcg32randomt ) :-
   44    pcg32_srandom_r( InitState, InitSeq, Pcg32randomt ).
   45
   46% abbreviated functor name for random number generator.
   47rg( pcg32randomt( 0x853c49e6748fea9b, 0xda3e39cb94b95bdb ) ).
   48rg( InitState, InitSeq, Pcg32randomt ) :-
   49    pcg32_srandom_r( InitState, InitSeq, Pcg32randomt ).
   50
   51% Random integer in [0,2^32-1].
   52random( RandomInt, Pcg32randomt, NextPcg32randomt ) :-
   53    pcg32_random_r( RandomInt, Pcg32randomt, NextPcg32randomt ).
   54
   55grandom(RandomInt) --> state( Pcg32randomt, NextPcg32randomt ),
   56    { pcg32_random_r( RandomInt, Pcg32randomt, NextPcg32randomt ) }.
   57
   58% Random integer in [0,Max-1].
   59random_bounded( Max, RandomInt, Pcg32randomt, NextPcg32randomt ) :-
   60    pcg32_boundedrand_r( Max, RandomInt, Pcg32randomt, NextPcg32randomt ).
   61
   62grandom_bounded( Max, RandomInt ) --> state( Pcg32randomt, NextPcg32randomt ),
   63    { pcg32_boundedrand_r( Max, RandomInt, Pcg32randomt, NextPcg32randomt ) }.
   64
   65% Random integer in [Min,Max].
   66random_between( Min, Max, RandomInt, Pcg32randomt, NextPcg32randomt ) :-
   67    #(Max) - #(Min) #< 2^32, #(Min) #< #(Max),
   68    #(Bound) #= #(Max) - #(Min),
   69    pcg32_boundedrand_r( Bound, R, Pcg32randomt, NextPcg32randomt ),
   70    #(RandomInt) #= #(R) + #(Min).
   71
   72grandom_between(Min, Max, RandomInt) --> state( Pcg32randomt, NextPcg32randomt ),
   73    { random_between( Min, Max, RandomInt, Pcg32randomt, NextPcg32randomt ) }.
   74
   75% Relationship between a natural number N and a list of random numbers L such that length(L,N).
   76
   77random_list( N, RL, Pcg32randomt, FinalPcg32randomt ) :-
   78    length(RL,N), foldl( random, RL, Pcg32randomt, FinalPcg32randomt ).
   79
   80grandom_list(N, RL) --> state( Pcg32randomt, NextPcg32randomt ),
   81    { random_list( N, RL, Pcg32randomt, NextPcg32randomt) }.
   82
   83% Relationship between a natural number N and a list of random bounded numbers L such that length(L,N).
   84random_boundedlist( Bound, N, RL, Pcg32randomt, FinalPcg32randomt ) :-
   85    length(RL,N), foldl( random_bounded(Bound), RL, Pcg32randomt, FinalPcg32randomt ).
   86
   87grandom_boundedlist( Bound, N, RL ) --> state( Pcg32randomt, NextPcg32randomt ),
   88    { random_boundedlist( Bound, N, RL, Pcg32randomt, NextPcg32randomt ) }.
   89
   90% Relationship between a natural number N and a list of random numbers between Min and Max) L such that length(L,N).
   91random_betweenlist( Min, Max, N, RL, Pcg32randomt, FinalPcg32randomt ) :-
   92    length(RL,N), foldl( random_between(Min,Max), RL, Pcg32randomt, FinalPcg32randomt ).
   93    
   94grandom_betweenlist( Min, Max, N, RL ) --> state( Pcg32randomt, NextPcg32randomt ),
   95    { random_betweenlist( Min, Max, N, RL, Pcg32randomt, NextPcg32randomt ) }.
   96
   97% T is reified - if an event with probability Probability/Precision happens, T = 1, else T = 0.
   98probability( Probability/Precision, T, Pcg32randomt, NextPcg32randomt ) :-
   99    pcg32_boundedrand_r( Precision, RandomNumber, Pcg32randomt, NextPcg32randomt ),
  100    #(RandomNumber) #=< #(Probability) #<==> #(T).
  101
  102gprobability( Probability/Precision, T ) --> state( Pcg32randomt, NextPcg32randomt ),
  103    { probability( Probability/Precision, T, Pcg32randomt, NextPcg32randomt ) }.
  104
  105random_permutation(List, RandomPermutation, Pcg32randomt, NextPcg32randomt) :-
  106    key_random(List, Keyed, Pcg32randomt, NextPcg32randomt),
  107    keysort(Keyed, Sorted),
  108    pairs_values(Sorted, RandomPermutation).
  109
  110grandom_permutation(List1, List2) --> state( Pcg32randomt, NextPcg32randomt ),
  111    { random_permutation(List1, List2, Pcg32randomt, NextPcg32randomt) }.
  112
  113% -------------------- private predicates ---------------------------:
  114
  115% random_permutation predicates:
  116
  117key_random([], [], Pcg32randomt, Pcg32randomt).
  118key_random([H|T0], [K-H|T], Pcg32randomt, FinalPcg32randomt) :-
  119    random(K, Pcg32randomt, NextPcg32randomt),
  120    key_random(T0, T, NextPcg32randomt, FinalPcg32randomt).
  121
  122% Auxiliary dcgs:
  123
  124state(S0, S), [S] --> [S0]. % implicitly pass states.
  125
  126% Random Number Generator
  127pcg32_srandom_r( InitState, InitSeq, FinalPcg32randomt ) :-
  128    Sup is 2^64 - 1, InitState in 0..Sup, InitSeq in 0..Sup,
  129    1 #= #(InitSeq) /\ 1,
  130    #(Inc) #= ( ( #(InitSeq) << 1) \/ 1 ) mod Sup,
  131    pcg32_random_r( _, pcg32randomt(0,Inc), pcg32randomt(S2,Inc) ),
  132    #(S3) #= ( #(S2) + #(InitState) ) mod Sup,
  133    pcg32_random_r( _, pcg32randomt(S3,Inc), FinalPcg32randomt ).
  134
  135% Uniformly distributed 32-bit random number:
  136pcg32_random_r( RandomValue, pcg32randomt(OldState,Inc), pcg32randomt(NextState,Inc) ) :-
  137    #(NextState) #= ( ( #(OldState) * 6364136223846793005) + (#(Inc) \/ 1) ) mod (2^64),
  138    #(XorShifted) #=( ( ( #(OldState) >> 18) xor #(OldState)) >> 27) mod (2^32), % here and in the line below mod 2^32 might be unnecessary, but i'm using it for precaution.
  139    #(Rot) #= ( #(OldState) >> 59) mod (2^32),
  140    #(RandomValue) #= ( (#(XorShifted) >> #(Rot)) \/ (#(XorShifted) << (- #(Rot) /\ 31)) ) mod (2^32).
  141
  142% Bounded random number using rejection sampling
  143pcg32_boundedrand_r( Bound, Result, Pcg32randomt, NextPcg32randomt ) :-
  144    Sup is 2^32, Bound in 1..Sup,
  145    #(Threshold) #= 2^32 mod #(Bound),
  146    bounded_sample( Bound, Threshold, Result, Pcg32randomt, NextPcg32randomt ).
  147
  148% Rejection sampling:
  149bounded_sample( Bound, Threshold, Result, Pcg32randomt, FinalPcg32randomt ) :-
  150    pcg32_random_r( R, Pcg32randomt, NextPcg32randomt ),
  151    (   #(R)