Bigint random numbers skewed?

General
1

I wonder whether bigint random numbers are
evenly generated. Here is my test case, the
expected value of generating random numbers

in the interval [-A,A+1) should be zero right?
Thats the query I am using, I am measuring
the bits of the “zero”, since it doesn’t reach zero,

thats only the theoretical expected value:

?- L is -(1<<100), H is (1<<100)+1, K is (1<<20), 
    between(1,10,_),  aggregate_all(sum(X),
         (between(1,K,_), random(L, H, X)), S), 
    M is msb(abs(S)>>20), write(M), nl, fail; true.

If I run this query with SWI-Prolog I get:

89
87
89
85
87
84
88
88
87
89

If I run this query with a JDK 19 based Prolog I get:

82
84
83
84
81
84
83
84
81
82

The second result looks more “zero” to me.

Edit 10.10.2022
I saw the skewing when fuzzying the float/1 function
on bigint, by generating 100 bit large bigints. By means
of this simple Prolog code:

fuzzer :-
  L is -(1<<100),
  H is (1<<100)+1,
  between(1, 100, N),
  random(L, H, X),
  Y is float(X),
  write(case(N, X, Y)), write('.'), nl,
  fail.
fuzzer.
2

Here is a better fuzzer, it has a higher mean lsb/1 value,
and therefore mostlikely tests more HALF_EVEN:

?- aggregate_all(sum(N), (case(_,X,_), N is lsb(abs(X))), S), A is S/100.
S = 94,
A = 0.94.

?- aggregate_all(sum(N), (case2(_,X,_), N is lsb(abs(X))), S), A is S/100.
S = 5050,
A = 50.5.

The fuzzer reads as follows:

fuzzer2 :-
  between(1, 100, N),
  L is -(1<<N),
  H is (1<<N)+1,
  random(L, H, X),
  Z is X<<(100-N),
  Y is float(Z),
  write(case2(N, Z, Y)), write('.'), nl,
  fail.
fuzzer2.
3

Any idea why the bigint random number is skewed in SWI-Prolog?
Nevermind. Here is yet another generator, it generates more
cases where SWI-Prolog fails. Since I found a case2/3 problem at

57, was using this value to generate more case2/3 problem similar
problems, now for facts case3/3:

fuzzer3 :-
   L is -(1<<57),
   H is (1<<57)+1,
   between(1, 100, N),
   random(L, H, X),
   Z is X<<43,
   Y is float(Z),
   write(case3(N, Z, Y)), write('.'), nl,
   fail.
fuzzer3.

The whole test cases, now case/3, case2/3 and case3/3 are here:

cases.p.log (21,0 KB)

4

All I know is that in the end it uses mpz_urandomm() from GMP and the initialization happens with gmp_randinit_mt() when available and gmp_randinit_default() otherwise.

5

Thanks. All these now work with bugfixed version:

?- case(N, X, Y), Y =\= float(X).
false.

?- case2(N, X, Y), Y =\= float(X).
false.

?- case3(N, X, Y), Y =\= float(X).
false.

I’ll generate a PR in the next day or so.

6

I guess I should do some more testing, is random/3 bootstrapped
from random/2? Yes it seems so, its the same what I am using on
my side in the other Prolog system now:

random(L, U, R) :-
    integer(L), integer(U),
    !,
    R is L+random(U-L).

https://www.swi-prolog.org/pldoc/doc/SWI/library/random.pl?show=src#random/3

If there is a skew in random/3, should also see it for random/2.

7

The low level operations are the arithmetic functions random(Max) and random_float. Given GMP, the integer version uses mpz_urandomm() and the float version mpf_urandomb(); The good old library(random) is built on top of these.

11

The library(random) has a very long history, dating back to the famous dec10 library. It has been extended by the YAP people and got SICStus compatible extensions. Originally it used a very simple algorithm. I added GMP based random numbers to arithmetic and reimplemented library(random) to use these random generators. AFAIK, I never changed the protocol.

Our dream … a standard library …

12

I totally agree!!

15

The error in Ciao Prolog is much larger than in SWI-Prolog,
i.e. the distance to the expected value. Was using this query,
taking into account the different random/3 protocoll, and

manually compensating for the missing aggregate_all/3
and msb/1 evaluable function:

$ ciaosh
Ciao 1.22.0 (2022-09-28 21:20:35 +0200) [LINUXx86_64]
?- use_module(library(between)).

?- ['aggregate.p'].

?- E is (1<<100), M is (1<<101), K is (1<<20),
between(1,10,_), aggregate_all(sum(X),
(between(1,K,_), random(0,M,X)), S),
H is abs(S-E*K), msb(H,J), N is J-20, write(N), nl, fail; true.
100.0
100.0
100.0
100.0
100.0
100.0
100.0
100.0
100.0
100.0

And this helper code (file aggregate.p):

:- use_module(library(numlists)).
:- use_module(library(random)).

aggregate_all(sum(X), G, S) :-
   findall(X, G, L),
   sum_list(L, S).

msb(X, Y) :- Y is log(X)/log(2).

Edit 13.10.2022
Also a little strange that it is exactly 100.0 in all runs?
Have to dig deeper whats going on.

17

For those who tried to understand what this is testing, I created an annotated version:

run :-
    run(100, 20).

run(Bits, LgTimes) :-
    E is 1<<Bits,                       % Expected avg + 0.5
    M is E*2,                           % M is twice as big
    K is 1<<LgTimes,                    % 1M
    (   between(1,10,_),                % 10 runs
        aggregate_all(sum(X),           % sum K random [0..M)
                      (between(1,K,_), random(0,M,X)),
                      S),
        H is abs(S-E*K),                % Expect H = K*0.5
        J is msb(H),                    % Expect J = LgTimes-1
        N is J-(LgTimes-1),             % Expect N = 0
        writeln(N),
        fail
    ;   true
    ).

This is a little different, subtracting 19 from J rather than 20. I think that is what we need to expect 0. That is a bit too simple though, we must expect this not to be precise. My statistics are a but too rusty to know what exactly we must expect from this.

For further reference, SWI-Prolog uses the GMP Mersenne Twister implementation. See Mersenne Twister - Wikipedia.

18

There is a small glitch, this here:

% sum K random [0..M)
random(0,M,X)

Should read:

/* SWI-Prolog Variant */
% sum K random [0..M1)
M1 is M+1, X is random(M1)

/* Ciao-Prolog Variant, it doesn't have evaluable so we go random/3 */
% sum K random [0..M]
random(0, M, X)

Edit 13.10.2022
But you could also use this code, on SWI-Prolog, where the interval can be anything:

/**
 * stats(M, D):
 * Compute the absolute difference from the average of 2^20 samples
 * in the interval [0..M) to the expected value, (M-1)/2 of an uniformly
 * distributed discrete random variable in the same interval.
 */
% stats(+Integer, -Float)
stats(M, D) :-
   K is (1<<20),
   aggregate_all(sum(X), (between(1,K,_), X is random(M)), S),
%   D is abs(S/K-(M-1)/2).
   D is abs((2*S-K*(M-1))/(2*K)).
19

The difference between the two systems gets smaller
with more bits. With M=2^51 the systems look quite different:

/* SWI-Prolog 8.5.18 */
?- M is (1<<51), between(1,10,_), stats(M, D),
N is log(D)/log(2), write(N), nl, fail; true.
38.030490564740845
36.64608133068678
39.90266908214614
39.40783704682056
31.502181228865066
38.46811325797751
37.23074324742856
36.82002406398337
34.9708016050052
38.75827431893174

/* Jekejeke Prolog 1.5.5 */
?- M is (1<<51), between(1,10,_), stats(M, D),
N is log(D)/log(2), write(N), nl, fail; true.
36.25036911731235
36.80713135437714
36.06365511352839
33.96122918542425
36.03296123667612
35.301308660385246
35.177861119830645
34.1454022656268
35.7980788891358
35.49236391391306

With M=2^201 the systems look rather similar:

?- M is (1<<201), between(1,10,_), stats(M, D),
N is log(D)/log(2), write(N), nl, fail; true.
189.53672010636757
188.95400489833128
187.3383068117556
188.6908371443576
188.22921137411072
186.30323985809352
188.6864997125651
189.79251475115166
189.2815717081842
190.26286276270744

/* Jekejeke Prolog 1.5.5 */
?- M is (1<<201), between(1,10,_), stats(M, D),
N is log(D)/log(2), write(N), nl, fail; true.
185.14260912573846
188.59895742142805
185.31313349159757
185.11803349540932
183.72143323996613
188.9781991830158
188.25081311791666
189.98099006923877
187.77570104217577
187.64532463688968

What does this mean? I guess the system with a smaller value is
possible the less good random generator, because some periodicity was hit.
Java uses a Lehmer algorithm which has a very short period and

uses very few bits in its state. On the other hand the Twister uses
quite some bits I guess, and has a much larger periodicity.
Twister would be therefore the better random number generator.

Edit 13.10.2022
Have to be careful to not get carried away, and should maybe introduce
some other testing on the random number generator. Main question is, can
I do a fuzzer on the SWI-Prolog platform so that it finds some float/1 test cases?

Actually this seems ok, so don’t worry! This was already affirmed in some other thread:

20

I like that :slight_smile: There are probably cases where such a skew is a problem, but I do not see it for random testing. Pseudo random numbers are complicated. From the Wikipedia page it seems that the twister is still recognized as good and widely used random generator.