Просмотр исходного кода

* fixed distribution of random(int64) based on patch by Pangea
(mantis #35878)

git-svn-id: trunk@42508 -

Jonas Maebe 6 лет назад
Родитель
Сommit
789f288771
4 измененных файлов с 73 добавлено и 9 удалено
  1. 2 0
      .gitattributes
  2. 27 9
      rtl/inc/system.inc
  3. 30 0
      tests/webtbs/tw35878.pp
  4. 14 0
      tests/webtbs/tw35878a.pp

+ 2 - 0
.gitattributes

@@ -17770,6 +17770,8 @@ tests/webtbs/tw3578.pp svneol=native#text/plain
 tests/webtbs/tw3579.pp svneol=native#text/plain
 tests/webtbs/tw3583.pp svneol=native#text/plain
 tests/webtbs/tw35862.pp svneol=native#text/pascal
+tests/webtbs/tw35878.pp svneol=native#text/plain
+tests/webtbs/tw35878a.pp svneol=native#text/plain
 tests/webtbs/tw35886.pp svneol=native#text/plain
 tests/webtbs/tw3589.pp svneol=native#text/plain
 tests/webtbs/tw3594.pp svneol=native#text/plain

+ 27 - 9
rtl/inc/system.inc

@@ -682,16 +682,34 @@ begin
 end;
 
 function random(l:int64): int64;
+var
+ a, b, c, d: cardinal;
+ q, bd, ad, bc, ac: qword;
+ carry: qword;
 begin
-  { always call random, so the random generator cycles (TP-compatible) (JM) }
-  { also do it in two separate steps, so the order in which the two calls
-    are performed is predictable (JM) }
-  random:=mtwist_u32rand;
-  random:=random or ((qword(mtwist_u32rand) shl 32) and high(int64));
-  if (l<>0) then
-    random := random mod l
-  else
-    random := 0;
+  if (l < 0) then
+    inc(l);
+ q:=qword(l);
+ if q>qword(low(int64)) then
+   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;
 
 {$ifndef FPUNONE}

+ 30 - 0
tests/webtbs/tw35878.pp

@@ -0,0 +1,30 @@
+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 18 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;
+  if abs(0.5-(s/n))>0.1 then
+    halt(1);
+end.

+ 14 - 0
tests/webtbs/tw35878a.pp

@@ -0,0 +1,14 @@
+var
+  i, rand, val: int64;
+begin
+  rand:=-5;
+  for i:=1 to 10000 do
+    begin
+      val:=random(rand);
+      if (val<-5) or (val>0) then
+        begin
+          writeln('error: ', val);
+          halt(1);
+        end;
+    end;
+end.