Output of Random(SomeInteger) appears to be biased lazarus / 64-FPC development ---------------------------------------------------------------------- Pangea · 2019-07-24 20:22 As suggested i open a forum thread for the problem discussed in https://bugs.freepascal.org/view.php?id=35878. If this is the wrong subforum, please move the thread. Background information: I noticed that the output of the 64-bit random function appears to be biased towards smaller numbers. I originally thought, it would be enough to adapt the same trick the 32 bit implementation uses to the 64 bit function and opened the bug tracker entry. But this was a mistake, as this also seems to be insufficient. It is possible to construct two really simple programs to show that both (the 32 and the 64-bit version) do not produce unbiased output. I apologize for this confusion. This is the program to reproduce this bug for the 64-bit version ("function Random(l:int64):int64;") Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i program random_test;{$mode objfpc}{$H+}{$APPTYPE CONSOLE}const   l: UInt64 = 6148914691236517205;var   s,n: UInt64;   i,j: UInt64;begin  WriteLn('Experiment:', LineEnding);  WriteLn(' Draw a natural number r from the intervall [0,l-1] and');  WriteLn(' increment a counter s when r < l div 2 is satisfied.');  WriteLn(' Repeat this step n times and calculate the ratio s/n.', LineEnding);   WriteLn(' Expected ratio: ', (l div 2)/l:30, LineEnding);   WriteLn('Input size n':16, 'Observed ratio s/n':30);  l := 6148914691236517205;  for j := 4 to 63 do  begin    n := (Int64(1) shl j);    s := 0;    for i := 0 to n-1 do      if Random(Int64(l)) < l div 2 then        Inc(s);    WriteLn( (UInt64(1) shl j):16, s/n:30);  end;end. And this is the code for the 32 bit version ("function Random(l:longint):longint;"): Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i program random_test;{$mode objfpc}{$H+}{$APPTYPE CONSOLE}var   s,n: longint;   i,j: longint;   l: LongInt;begin  l := 1227133513;   WriteLn('Experiment:', LineEnding);  WriteLn(' Draw a natural number r from the intervall [0,l-1] and');  WriteLn(' increment a counter s when r mod 2 = 0 is satisfied.');  WriteLn(' Repeat this step n times and calculate the ratio s/n.', LineEnding);   WriteLn(' Expected ratio: ', (l div 2)/l:30, LineEnding);   WriteLn('Input size n':16, 'Observed ratio s/n':30);  for j := 4 to 31 do  begin    n := (LongInt(1) shl j);    s := 0;    for i := 0 to n-1 do      if Random(LongInt(l)) mod 2 = 0 then        Inc(s);    WriteLn( (LongInt(1) shl j):16, s/n:30);  end;end. Both programs have an expectation value of 1/2, but the current random implementation does not converge to this. My take on this matter: This is not a problem of the underlying random number generator, but the way the output range of the PRNG is mapped onto the smaller interval. Therefore it is completly possible for the output of this functions to be biased, while the PRNG reproduces the values of the reference implementation. The bias is not introduced in the generator, but in the map. I also think it is necessary to adapt rejection sampling to guarantee an unbiased result. Feel free to tear this apart. : ) Edit: I tested the programs against version 3.0.4. ---------------------------------------------------------------------- Martin_fr · 2019-07-24 23:08 Actually you appear to be right. From system unit: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i function random(l:longint): longint;begin  { otherwise we can return values = l (JM) }  if (l < 0) then    inc(l);  random := longint((int64(cardinal(genrand_MT19937))*l) shr 32);end; function random(l:int64): int64;begin  { always call random, so the random generator cycles (TP-compatible) (JM) }  random := int64((qword(cardinal(genrand_MT19937)) or ((qword(cardinal(genrand_MT19937)) shl 32))) and $7fffffffffffffff);  if (l<>0) then    random := random mod l  else    random := 0;end;  genrand_MT19937 returns a random 32 bit value. Lets assume that part to be as correct as can be. There is a generic problem. If you do "Random(100)" you need to map the 2^32 values to 0..99. But 2^31 is not int dividable by 100. So some of 0..99 will be more likely. Once you evenly distributed ‭4,294,967,200 out of the ‭4,294,967,296‬ values, you have 96 values left, to evenly distribute. How to you map 96 values to 0..99 so that each of 0..99 is equally likely chosen? This problem then manifests in both of the implementations 64bit:     random := random mod l means that  96 to 99 are slightly less likely. (If it was QWord) With random (2^63 + n), the last 2*n values are only half as likely to occur. 32bit: The numbers with different likelihood are equally spread across the requested range. -- The issue is that any stable function to map the result of genrand_MT19937 to a given range will have the same problem. You would need additional inputs to map them different (e.g counting how often "random" was called.) But such additional input may be guess-able. And if it is part of the random result, then an attacker may be able to make predictions based on that. Of course if an attacker knows your range, and that you use the fpc implementation, said attacker can make guesses too. -- Since using "Mod" the bias on the numbers depends on how big the range is (relative to the max of genrand_MT19937 ), you can get a bigger random number, by calling genrand_MT19937  several times. It still will be biased, but statistically less likely. ---------------------------------------------------------------------- Thaddy · 2019-07-25 05:12 It can be right, it can likely be solved by scaling back and forth by the forrmulea -slightly adjusted - I gave and its inverse. That will give a best fit but is slow. The algorithm itself, however is a standard. Still: we can't change that, because that is how it worked since TP days.  Existing data will become unusable and that is not acceptable. (I have quite a lot of such data dating as far back as its introduction in FreePascal. There are billions of such data sets, because the algorithm is so widely used in multiple languages. My suggestion is, and the focus of a solution I will propose on the mailing list, is to overload and where not possible add differently named functionality. Caveat is a platform without an FPU must support soft float and that is not always the case so that is frankly not a solution. I will investigate if it is possible to implement the first suggestion, but only for ranged randoms: these are not part of the actual standard. The good news is I have the 64 bit MT working and I am now testing if the output is equal to the original test data at http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt19937-64.out.txt. If that works I will work on the problem we discuss here so it can be solved for both 32 an 64 bit. Reason is that I want the 64 bit version standards compliant first, then add it to my compendium and then propose a patch. Non- Xplatform solutions will not be considered by me. ---------------------------------------------------------------------- ASerge · 2019-07-25 06:01 For those who want a quick workaround, I can suggest adding this to the beginning of the test: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i function Random(X: SizeInt): SizeInt;begin  Result := Round(System.Random * X);end; ---------------------------------------------------------------------- Pangea · 2019-07-25 09:15 Quote from: ASerge on July 25, 2019, 08:01:37 amFor those who want a quick workaround, I can suggest adding this to the beginning of the test: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i function Random(X: SizeInt): SizeInt;begin  Result := Round(System.Random * X);end; This still fails when SizeInt has 64 bit and L is large. As a test please consider this adapted version of the 64-bit program (i am using the 64-bit compiler, so sizeint has 64-bit): Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i program random_test; {$mode objfpc}{$H+}{$APPTYPE CONSOLE} function RoundRandom(x: SizeInt): SizeInt;begin  Result := round(x*Random);end; var   s,n: Int64;   i,j: Int64;   l: Int64;begin  l := 6148914691236517205;   WriteLn('Experiment:', LineEnding);  WriteLn(' Draw a natural number r from the intervall [0,l-1] and');  WriteLn(' increment a counter s when r mod 2 = 0 is satisfied.');  WriteLn(' Repeat this step n times and calculate the ratio s/n.', LineEnding);   WriteLn(' Expected ratio: ', (l div 2)/l:30, LineEnding);   WriteLn('Input size n':16, 'Observed ratio s/n':30);  for j := 4 to 31 do  begin    n := (Int64(1) shl j);    s := 0;    for i := 0 to n-1 do      if RoundRandom(Int64(l)) mod 2 = 0 then        Inc(s);    WriteLn( (Int64(1) shl j):16, s/n:30);  end;end.     So here is an variant using rejection sampling (no floating point arithmetic, just reject the samples you cannot map correctly): Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i   function RandomRejectionSampling(L: UInt64): Uint64;  var    thr, rem: UInt64;  begin    if L = 0 then      Result := 0    else    begin      // We need to calculate High(UInt64)+1 mod L without using 128 bit arithmetic.      rem := (High(UInt64) mod L) + 1;      if rem = L then        rem := 0;      // Calculate the threshold      thr := (High(UInt64) - rem);      if thr < L-1 then        thr := L-1;      // Perform rejection sampling      repeat        Result := ((UInt64(GenRand32) shl 32) or UInt64(GenRand32));      until Result <= thr;      Result := Result mod L;    end;  end; It's just the very basic idea, so there might be room for optimization. Edit: Fixed a bug in the code. ---------------------------------------------------------------------- ASerge · 2019-07-25 09:25 Quote from: Pangea on July 25, 2019, 11:15:04 amThis still fails when SizeInt has 64 bit and L is large. As a test please consider this adapted version of the 64-bit program (i am using the 64-bit compiler, so sizeint has 64-bit): And where fails? I've added my function to both your implementations. Everything was compiled and started without problems. And the distribution became about 0.5. ---------------------------------------------------------------------- Pangea · 2019-07-25 09:41 It only produces even numbers for large L. You can test it against the 64 bit program when you change the condition to random(l) mod 2 = 0. ---------------------------------------------------------------------- Martin_fr · 2019-07-25 10:27 QuoteCode: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i      // Perform rejection sampling      repeat        Result := ((UInt64(GenRand32) shl 32) or UInt64(GenRand32));      until Result <= thr; Depending on the underlying random, you may have to reject a lot of numbers... Also this makes the random function, take different amount of time depending on where in the underlaying random sequence it is. Are you sure this will not open the door for timing attacks? ---------------------------------------------------------------------- Pangea · 2019-07-25 11:38 I attached a plot showing the acceptance probability over the choice L. The acceptance should never be smaller than 1/2 (unless there is a bug). My personal opinion on the timing part is: The underlying PRNG is not cryptographically secure in any way, so securing the random number generation against side channel attacks does not really help. ---------------------------------------------------------------------- ASerge · 2019-07-25 11:47 Quote from: Pangea on July 25, 2019, 11:41:58 amIt only produces even numbers for large L. You can test it against the 64 bit program when you change the condition to random(l) mod 2 = 0. That's expected. The accuracy of fraction part for double is 52, for large numbers only the exponent is used. There is no such problem for x32, of course. ---------------------------------------------------------------------- Thaddy · 2019-07-25 14:03 Quote from: Pangea on July 25, 2019, 01:38:45 pmMy personal opinion on the timing part is: The underlying PRNG is not cryptographically secure in any way, so securing the random number generation against side channel attacks does not really help. It is not secure, but it is compiled in B+ mode, so Boolean shortcuts have no timing influence because every expression is fully evaluated.. This was done on my request. [edit] This is not true for random, but for the hashes. But timing attacks can be resolved by B+. ---------------------------------------------------------------------- Pangea · 2019-07-25 16:53 I don't think it is necessary to harden random against timing attacks. It is not designed to be used as cryptographic primitive. I eliminated an unnecessary comparison and this is, what i would propose: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i   function Random(L: UInt32): UInt32; overload;  var    thr, rem: UInt32;  begin    if L = 0 then      Result := 0    else    begin      rem := (UInt64(High(UInt32)) + 1) mod L;      thr := (High(UInt32) - rem);    repeat        Result := genrand_MT19937;      until Result <= thr;    Result := Result mod L;    end;  end;        function Random(L: UInt64): UInt64; overload;  var    thr, rem: UInt64;  begin    if L = 0 then      Result := 0    else    begin      // We need to calculate High(UInt64)+1 mod L without using 128 bit arithmetic.      rem := (High(UInt64) mod L) + 1;      if rem = L then        rem := 0;      // Calculate the threshold      thr := (High(UInt64) - rem);      repeat        Result := ((UInt64(genrand_MT19937) shl 32) or UInt64(genrand_MT19937));      until Result <= thr;      Result := Result mod L;    end;  end;  However this does NOT behave like the old implementation. This is especially relevant for former negative input. If negative input is desired, maybe this is more appropriate (Also the code is simpler): Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i function Random(L: Int32): Int32; overload;const  UpperBound: UInt32 = (High(UInt32) shr 1); var   thr, rem: UInt32; begin   if L <= 0 then   Result := 0   else   begin     rem := UInt32(UpperBound+1) mod UInt32(L);     thr := (UpperBound - rem);     repeat       Result := genrand_MT19937 and UpperBound;     until UInt32(Result) <= thr;     Result := Result mod L;   end; end; function Random(L: Int64): Int64; overload;const  UpperBound: UInt64 = (High(UInt64) shr 1); var   thr, rem: UInt64; begin   if L <= 0 then     Result := 0   else   begin     rem := (UpperBound +1) mod UInt64(L);     thr := (UpperBound - rem);     repeat       Result :=  ((UInt64(genrand_MT19937) shl 32) or UInt64(genrand_MT19937)) and UpperBound;     until UInt64(Result) <= thr;     Result := Result mod L;   end; end;     However this still behaves differently for neagtive L in the 64 bit case. Maybe it should be clarified in the documentation what happens when L is smaller or euqal to zero. Edit: I read the old code wrong. If it is negative it should return a negative number in the respective range. So this is still insufficient. Edit: It is too hot. I can't think clearly and i do not trust myself. Currently random(Int32(-5)) seems to return numbers in the interval [-4; 0], but random(Int64(-5)) returns numbers in the interval [0, 4]??? Can someone clarify this please? ---------------------------------------------------------------------- PascalDragon · 2019-07-26 07:23 Quote from: Pangea on July 25, 2019, 06:53:04 pmI don't think it is necessary to harden random against timing attacks. It is not designed to be used as cryptographic primitive.Correct. That's not its purpose. ---------------------------------------------------------------------- devEric69 · 2019-07-26 08:31 If it's a question of having a real random number generator (a long sequence of numbers, stored and studied in a spreadsheet program, to see if it's roughly (no computer generator is perfect: it's a question of numbers generated, before it becomes biased. So, the real question is: does the count, of random numbers that I will use, remain more or less uniformly distributed in the considered range of random numbers that interests me, with this generator?) formally distributed according to a bell curve, over a given range of numbers, centered, and according to a normal or other statistical law), then it's preferable to use a more advanced algorithm like Tore, or another one. My 2 cents. ---------------------------------------------------------------------- Pangea · 2019-07-27 19:53 Quote from: devEric69 on July 26, 2019, 10:31:15 amit's preferable to use a more advanced algorithm like Tore I googled this algorithm and apparently this is a way to construct low discrepancy sequences by taking the fractional part of multiplies of irrational numbers (e.g. the square roots of primes). Is this correct? I also found [1] and [2], both giving an overview over different sampling methods. I went ahead and implemented them together with a simple check, wether their output is really unbiased (a.k.a. i did not implement a bug). In my opinion the one using the intrinsic (no division) and the one with only a single division are really interesting. [1]  https://arxiv.org/pdf/1805.10941 [2]  http://www.pcg-random.org/posts/bounded-rands.html Edit: Removed source. ---------------------------------------------------------------------- Thaddy · 2019-07-27 20:30 Did you test with at least BigCrush? ---------------------------------------------------------------------- Pangea · 2019-07-27 21:48 I checked TestU01, but i seems like it only accepts powers of 2 as the output range of your PRNG. So i was not able to run the samplers against it. Edit: If the upper range of the interval is divides the upper bound of the generator, no bias for any sampler is expected. Edit2: Unless, of course, the underlying PRNG is also biased. Edit3: Actually, I can convert the output of the sampler to some floating point number and assume some precision. Then i can feed to the testsuit. Maybe i try this. Edit4: Ok, i tested the samplers against SmallCrush. Find the details and results below. However i don't think i will run this against BigCrush, because completing this takes forever and i cannot run my PC overnight. Also i have read somwhere, that the Mersenne Twister fails some test in big crush? Would this mean it will also fail for the composite number, even when the rejecting maps are working correctly? Results: As expected both naive non rejecting methods fail several tests. To be precise: Code: [Select]========= Summary results of SmallCrush =========  Version:          TestU01 1.2.3  Generator:        TMulSampler  Number of statistics:  15  Total CPU time:   00:00:41.60  The following tests gave p-values outside [0.001, 0.9990]:  (eps  means a value < 1.0e-300):  (eps1 means a value < 1.0e-15):        Test                          p-value  ----------------------------------------------   3  Gap                              eps    4  SimpPoker                        eps    5  CouponCollector                  eps    7  WeightDistrib                    eps    9  HammingIndep                     eps   10  RandomWalk1 H                    eps   10  RandomWalk1 M                    eps   10  RandomWalk1 J                    eps   10  RandomWalk1 C                   2.5e-5  ----------------------------------------------  All other tests were passed ========= Summary results of SmallCrush =========  Version:          TestU01 1.2.3  Generator:        TModSampler  Number of statistics:  15  Total CPU time:   00:00:41.26  The following tests gave p-values outside [0.001, 0.9990]:  (eps  means a value < 1.0e-300):  (eps1 means a value < 1.0e-15):        Test                          p-value  ----------------------------------------------   2  Collision                        eps    6  MaxOft                           eps    6  MaxOft AD                      1 - eps1   7  WeightDistrib                    eps    9  HammingIndep                   5.6e-16  10  RandomWalk1 H                    eps   10  RandomWalk1 M                    eps   10  RandomWalk1 J                    eps   10  RandomWalk1 R                    eps   10  RandomWalk1 C                    eps   ----------------------------------------------  All other tests were passedTMulSampler refers to the multiply and shift method and TModSampler refers to the modulo maping, both (at least in a somewhat similar form) currently implemented in the RTL.  All other implemented rejecting methods pass SmallCrush without failing any test. I attached the full output log to this post. Used method: To feed the output of the samplers into testsuite i construct a floating point number iterative (via Random(L)/(L-1) + Random(L)/(L-1)^2 + ...) until the number has absorbed at least 52 random bits. For the test i used a generator width of 31 bit (backed by the mersenne twister generator in the RTL) and target output interval size of L= 1610612736. I attached the source code to this post (The code is a bit ugly, but well...). To build it you need fpc, gcc, testu01 installed and all header files and libraries in the respective search path. First call "gcc -c bat.c" and after this build the pascal application. If you do not want to run against testu01, remove the testu01 define in the pascal program and it should compile on any fpc installation. Since all samplers are tested against testu01 directly after each other it is helpful to pipe the output in some file and just watch it with tail -f or something. Feel free to run it against big or small crush yourself. I will also remove the attached code on the older forums post, to avoid having two different sources online. Edit5: Actually it should be Random(L)/L + Random(L)/L^2 + ... . I will rerun the test with the correct formula .... Edit6: Ok, it yields basically the same results. I updated  the attachments. ---------------------------------------------------------------------- Pangea · 2019-07-31 02:03 I have seen, that the issue has been marked as resolved. However i still see two minor problems with this: - The generated output interval for negative input appears to be one to small? Or is this intended? - The output is still biased. What is the stance on this? Is the bias in this version of the function deemed acceptable (since it has the same bias the 32 bit version has)? This is the current code in trunk: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i function random(l:int64): int64;var a, b, c, d: cardinal; q, bd, ad, bc, ac: qword; carry: qword;begin  if (l < 0) then    begin      inc(l);      q:=qword(-l)    end  else    q:=qword(l);   a:=mtwist_u32rand; b:=mtwist_u32rand;  c:=q shr 32; d:=cardinal(q);  bd:=b*d; ad:=a*d; bc:=b*c; ac:=a*c;  // We only need the carry bit carry:=((bd shr 32)+cardinal(ad)+cardinal(bc)) shr 32;  // Calculate the final result result:=int64(carry+(ad shr 32)+(bc shr 32)+ac); if l<0 then   result:=-result;end;