Besseres Random() - eure Vorschläge

[Frickler]

Haste die Generatoren einlaufen lassen? Also erstmal jeweils ein paar tausend Werte abgefragt, bevor Du ihn benutzt?

[Stevie]

Solang eh nur 32bit Zahlen rauskommen isses eh Mumpitz, damit kannst nicht mal nen Skatspiel (!32) vernünftig mischen.

[DieDolly]

Was beim Mersenne übrigens scheinbar kein bisschen zu funktionierern scheint ist es Integer einfach in Int64 umzuschreiben
Hier mal die komplette Unit wenn die jemanden interessiert

Delphi-Quellcode:

unit Shared.Math.Random.MT19937;

{$R-} {range checking off}
{$Warnings OFF}
// {$Q-} {overflow checking off}
{----------------------------------------------------------------------
 Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
 Pseudo-Random Number Generator.

 What is Mersenne Twister?
 Mersenne Twister(MT) is a pseudorandom number generator developped by
 Makoto Matsumoto and Takuji Nishimura (alphabetical order) during
 1996-1997. MT has the following merits:
 It is designed with consideration on the flaws of various existing
 generators.
 Far longer period and far higher order of equidistribution than any
 other implemented generators. (It is proved that the period is 2^19937-1,
 and 623-dimensional equidistribution property is assured.)
 Fast generation. (Although it depends on the system, it is reported that
 MT is sometimes faster than the standard ANSI-C library in a system
 with pipeline and cache memory.)
 Efficient use of the memory. (The implemented C-code mt19937.c
 consumes only 624 words of working area.)

 home page
 http://www.math.keio.ac.jp/~matumoto/emt.html
 original c source
 http://www.math.keio.ac.jp/~nisimura/random/int/mt19937int.c

 Coded by Takuji Nishimura, considering the suggestions by
 Topher Cooper and Marc Rieffel in July-Aug. 1997.

 This library is free software; you can redistribute it and/or
 modify it under the terms of the GNU Library General Public
 License as published by the Free Software Foundation; either
 version 2 of the License, or (at your option) any later
 version.
 This library is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 See the GNU Library General Public License for more details.
 You should have received a copy of the GNU Library General
 Public License along with this library; if not, write to the
 Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 02111-1307  USA

 Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura.
 When you use this, send an email to: [e-mail]
 with an appropriate reference to your work.

 REFERENCE
 M. Matsumoto and T. Nishimura,
 "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
 Pseudo-Random Number Generator",
 ACM Transactions on Modeling and Computer Simulation,
 Vol. 8, No. 1, January 1998, pp 3--30.


 Translated to OP and Delphi interface added by Roman Krejci (6.12.1999)

 http://www.rksolution.cz/delphi/tips.htm

 Revised 21.6.2000: Bug in the function RandInt_MT19937 fixed
 ----------------------------------------------------------------------}

interface

{Period parameter}
const
 MT19937N = 624;

type
 TMT19937StateArray = array [0 .. MT19937N - 1] of Int64; // the array for the state vector

procedure sgenrand_MT19937(seed: Int64); // Initialization by seed
procedure lsgenrand_MT19937(const seed_array: TMT19937StateArray); // Initialization by array of seeds
procedure randomize_MT19937; // randomization
function randInt_MT19937(Range: Int64): Int64; overload; // Int64 Random with positive range
function randInt_MT19937(Min, Max: Int64): Int64; overload; // Int64 Random with positive range between min and max
function genrand_MT19937: Int64; // random Int64 (full range);
// function randFloat_MT19937: Double; // float RANDOM on 0..1 interval

implementation

{Period parameters}
const
 MT19937M = 397;
 MT19937MATRIX_A = $9908B0DF; // constant vector a
 MT19937UPPER_MASK = $80000000; // most significant w-r bits
 MT19937LOWER_MASK = $7FFFFFFF; // least significant r bits

 // Tempering parameters
 TEMPERING_MASK_B = $9D2C5680;
 TEMPERING_MASK_C = $EFC60000;

type
 Int64Rec = packed record
  case Integer of
   0:
    (Lo, Hi: Cardinal);
   1:
    (Cardinals: array [0 .. 1] of Cardinal);
   2:
    (Words: array [0 .. 3] of Word);
   3:
    (Bytes: array [0 .. 7] of Byte);
 end;

var
 mt: TMT19937StateArray;
 mti: Int64 = MT19937N + 1; // mti=MT19937N+1 means mt[] is not initialized

 // Initializing the array with a seed
procedure sgenrand_MT19937(seed: Int64);
var
 i: Int64;
begin
 for i := 0 to MT19937N - 1 do
  begin
   mt[i] := seed and $FFFF0000;
   seed := 69069 * seed + 1;
   mt[i] := mt[i] or ((seed and $FFFF0000) shr 16);
   seed := 69069 * seed + 1;
  end;

 mti := MT19937N;
end;

{
 Initialization by "sgenrand_MT19937()" is an example. Theoretically,
 there are 2^19937-1 possible states as an intial state.
 This function (lsgenrand_MT19937) allows to choose any of 2^19937-1 ones.
 Essential bits in "seed_array[]" is following 19937 bits:
 (seed_array[0]&MT19937UPPER_MASK), seed_array[1], ..., seed_array[MT19937-1].
 (seed_array[0]&MT19937LOWER_MASK) is discarded.
 Theoretically,
 (seed_array[0]&MT19937UPPER_MASK), seed_array[1], ..., seed_array[MT19937N-1]
 can take any values except all zeros.
}

procedure lsgenrand_MT19937(const seed_array: TMT19937StateArray);
var
 i: Int64;
begin
 for i := 0 to MT19937N - 1 do
  mt[i] := seed_array[i];

 mti := MT19937N;
end;

function genrand_MT19937: Int64;
const
 mag01: array [0 .. 1] of Int64 = (0, MT19937MATRIX_A);
var
 y, kk: Int64;
begin
 if mti >= MT19937N then // generate MT19937N Int64 at one time
  begin
   if mti = (MT19937N + 1) then // if sgenrand_MT19937() has not been called,
    sgenrand_MT19937(4357); // default initial seed is used

   for kk := 0 to MT19937N - MT19937M - 1 do
    begin
     y := (mt[kk] and MT19937UPPER_MASK) or (mt[kk + 1] and MT19937LOWER_MASK);
     mt[kk] := mt[kk + MT19937M] xor (y shr 1) xor mag01[y and $00000001];
    end;

   for kk := MT19937N - MT19937M to MT19937N - 2 do
    begin
     y := (mt[kk] and MT19937UPPER_MASK) or (mt[kk + 1] and MT19937LOWER_MASK);
     mt[kk] := mt[kk + (MT19937M - MT19937N)] xor (y shr 1) xor mag01[y and $00000001];
    end;

   y := (mt[MT19937N - 1] and MT19937UPPER_MASK) or (mt[0] and MT19937LOWER_MASK);
   mt[MT19937N - 1] := mt[MT19937M - 1] xor (y shr 1) xor mag01[y and $00000001];
   mti := 0;
  end;

 y := mt[mti];
 Inc(mti);
 y := y xor (y shr 11);
 y := y xor (y shl 7) and TEMPERING_MASK_B;
 y := y xor (y shl 15) and TEMPERING_MASK_C;
 y := y xor (y shr 18);

 Result := y;
end;

// Delphi interface

procedure randomize_MT19937;
var
 OldRandSeed: Integer;
begin
 OldRandSeed := System.RandSeed; // save system RandSeed value
 System.Randomize; // randseed value based on system time is generated
 sgenrand_MT19937(System.RandSeed); // initialize generator state array
 System.RandSeed := OldRandSeed; // restore system RandSeed
end;

// rewritten to delphi code due to 64bit compatibility
// Uwe Raabe -  delphipraxis.com
function randInt_MT19937(Range: Int64): Int64;
begin
 Result := Int64Rec(Int64(Range) * Cardinal(genrand_MT19937)).Hi;
end;
// bug fixed 21.6.2000.
// asm
// PUSH EAX // Parameter 1 = Range auf Stack sichern
// CALL genrand_MT19937 // Random Int64 aufrufen, Ergebnis steht in EAX (keine Parameter)
// POP EDX // Range vom Stack in EDX holen
// MUL EDX // EAX mit EDX multiplizieren, Ergebnis steht in EDX,EAX
// MOV EAX,EDX // Die h�heren 32-Bits in EDX als Result in EAX zur�ckgeben
// end;

function randInt_MT19937(Min, Max: Int64): Int64;
begin
 Result := randInt_MT19937(Max - Min) + Min;
end;

// function randFloat_MT19937: Double;
// const
// Minus32: Double = -32.0;
// asm
// CALL    genrand_MT19937
// PUSH    0
// PUSH    EAX
// FLD     Minus32
// FILD    qword ptr [ESP]
// ADD     ESP,8
// FSCALE
// FSTP    ST(1)
// end;

initialization

// added 25.11.2019
randomize_MT19937;

end.

[Michael II]

Zitat von Stevie:

Solang eh nur 32bit Zahlen rauskommen isses eh Mumpitz, damit kannst nicht mal nen Skatspiel (!32) vernünftig mischen.

Ja wenn Z die Menge aller Zufallszahlen und V die Menge aller Kartenverteilungen ist und du jedem z aus Z eine Kartenverteilung v aus V zuordnen willst, dann ist die Funktion F:Z->V natürlich nicht surjektiv. (V mächtiger als Z.)

Das ist aber kein Problem, da du eine beliebige Anzahl von Zufallszahlen erzeugen kannst und so immer eine Funktion F:ZxZ..xZ -> V findest, welche alle Elemente von V trifft.

Es reicht also bereits ein 1 Bit Zufallsgenerator.

Problematisch ist u.U. die Periodenlänge: Bei vielen P-Zufallsgeneratoren werden ausgehend von einem Startwert s (von Startwerten s[i]) n Bit P-Zufallszahlen berechnet. Irgendwann wird dabei wieder der Startwert erreicht und alles beginnt von vorne. Dabei werden oft nicht alle möglichen 2^n Werte durchlaufen.

[Michael II]

Kennt jemand eine Delphi "Übersetzung" von David Zuckermans und Eshan Chattopadhyays Arbeit?

Infos:

https://www.chip.de/news/Diese-Nachr..._97383944.html

https://news.utexas.edu/2016/05/16/c...cybersecurity/

Paper:

https://eccc.weizmann.ac.il//report/2015/119/

[markus5766h]

Bei Ansicht dieses Artikels sind mir meine alten
Spielereien zur Zahl Pi eingefallen . . . hier
die Monte-Carlo-Methode ...

Beim Screenshot im Anhang wurde der blaue Graph mit der Methode aus Post #30 erzeugt,
beim roten Graph mit Randomize ... x:=Random, y := Random

Die Unterschiede sind mit steigender Zahl der Durchläufe eher als vernachlässigbar zu bezeichnen . . .

Die Graphen wurden mittels 2500000 Durchläufen zur Berechnung von Pi erstellt.

[Michael II]

Post #30 zeigt wie Delphi (in 10.3.3) random() berechnet.

Es ist also nicht erstaunlich, wenn random und random zu einem ähnlichen (schlechten Zufalls-)Resultat führen. Bei gleichem Startwert randseed sollten sogar die gleichen Zahlenreihen entstehen.

Interessant wäre wie bereits erwähnt eine Übersetzung der Arbeit von Zuckerman & Chattopadhyay (Link weiter oben) ins Delphische.

[himitsu]

Inzwischen sind die Funktionen "auch" als Pure-Pascal verfügbar (für andere Platformen wo dieser Assembler nicht läuft)
und was ich grad bemerkte und als besonders praktisch empfand, wurden dort irgendwann zwei Variablen eingebaut, so dass man da ein eigenes Random (IntRandom) und Randomize registrieren und es somit auch "überall" im Delphi austauschen kann.

PS: Das aus #30 noch ein bissl zusammengefasst, dann kommt das bei raus.

Delphi-Quellcode:

function ra(aRange: Integer): Integer;
begin
  RandSeed := RandSeed * 134775813 + 1; // 3*17*131*20173
  Result := (Int64(aRange) * RandSeed) shr 32;  // aka Result := MulDiv(aRange, RandSeed, $100000000);
end;

Aber die Grundlagen des "deterministic linear congruential generator with 134775813 as a and 1 as c" haben sich nicht geändert.
Das ist seit über 20 Jahren konstant. Nur wurde vor längerer Zeit mal das GetTickCount gegen QueryPerformanceCounter im Randomize ersetzt, womit ein "falsch" benutztes und zu oft/schnell aufgerufenes Randomize keine großen Nachteile bringt.

Delphi-Quellcode:

type
  TRandom32Proc = function: UInt32;
  TRandomizeProc = procedure(NewSeed: UInt64);

function DefaultRandom32: UInt32;
procedure DefaultRandomize(NewSeed: UInt64);

var
  Random32Proc: TRandom32Proc = DefaultRandom32;
  RandomizeProc: TRandomizeProc = DefaultRandomize;

procedure Randomize;

function Random(const ARange: Integer): Integer; overload;
function Random: Extended; overload;



function DefaultRandom32: UInt32;
{$IFDEF PUREPASCAL}
begin
  Result := UInt32(RandSeed) * $08088405 + 1;
  RandSeed := Result
end;
{$ELSE !PUREPASCAL}
asm
{     <-EAX     Result  }
{$IFDEF PIC}
        PUSH    EBX
        CALL    GetGOT
        MOV     EBX,EAX
        MOV     ECX,[EBX].RandSeed
        IMUL    EAX,[ECX],08088405H
        INC     EAX
        MOV     [ECX],EAX
        POP     EBX
{$ELSE !PIC}
        IMUL    EAX,RandSeed,08088405H
        INC     EAX
        MOV     RandSeed,EAX
{$ENDIF !PIC}
end;
{$ENDIF !PUREPASCAL}

function Random(const ARange: Integer): Integer;
{$IFDEF PUREPASCAL}
var
  Temp: UInt32;
begin
  Temp := Random32Proc;
  Result := (UInt64(UInt32(ARange)) * UInt64(Temp)) shr 32;
end;
{$ELSE !PUREPASCAL}
asm
{     ->EAX     Range   }
{     <-EAX     Result  }
{$IFDEF PIC}
        PUSH    EBX
        PUSH    EAX
...

[Michael II]

Es ist mir bewusst, dass der Code aus #30 zusammengefasst werden kann. Es ging ja darum den ASM Code Zeile für Zeile in Delphi zu schreiben.

Deine Zusammenfassung von #30 liefert bei mir bei einem Startwert von randseed=1 bereits beim zweiten Aufruf von random(1000000) nicht mehr den korrekten Wert zurück. Es liegt wahrscheinlich an meinen Compilerschaltern (?). Wenn du's so geschrieben hättest...

Delphi-Quellcode:

function random(aRange: Integer): Integer;
begin
{$Q-}
{$R-}
  RandSeed := RandSeed * 134775813 + 1; // 3*17*131*20173
  Result := (aRange * Cardinal(RandSeed)) shr 32; // aka Result := MulDiv(aRange, RandSeed, $100000000);
end;

...dann wäre auch mein australisches Delphi zufrieden.

Infos zum Verfahren:
https://de.wikipedia.org/wiki/Kongruenzgenerator

[Rollo62]

@himitsu
QueryPerformanceCounter ist aber auch nur für Windows relevant ...

Ich habe jetzt nicht die hohen Ansprüche daran herauszufinden ob der Zufall in unserem Universum überhaupt existiert und gleichverteilt geblieben ist über die letzten 13 Mrd. Jahre,
aber gut dass die Funktionen jetzt im Source vorhanden sind, ist eine gute Sache.