Entwickler-Ecke

Algorithmen, Optimierung und Assembler - Zufallsgenerator mit 64-Bit Seed


Spaceguide - Fr 09.09.05 18:45
Titel: Zufallsgenerator mit 64-Bit Seed
Hat jemand einen Zufallsgenerator, der Int64 als Seed-Variable frisst und auch möglichst lange Perioden hat?


DBR - Di 13.09.05 20:06

in der Unit Math gibst das:


Delphi-Quelltext
1:
function RandG(Mean, StdDev: Extended): Extended;                    


Schau halt mal in die Hilfe.

Gruß DBR


Spaceguide - Di 13.09.05 20:21

user profile iconDBR hat folgendes geschrieben:
in der Unit Math gibst das:
Schau halt mal in die Hilfe.

Sieht so aus, als ob du nicht in die Hilfe geschaut hast. Nur weil da Häufungspunkt und Abweichung im Format extended vorliegen, heisst das noch lange nicht, dass der Seed (wird wohl die globale RandSeed Variable sein, die ist ein Longint!) 64-bit lang ist.


alzaimar - Di 13.09.05 20:55

Blöde Idee (Sollte aber ordendlich funktionieren):

Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
Function Random64 : Int64;
Type 
  TMap64 = Record 
    i0, i1 : Integer 
    End;
Begin
  TMap64(Result).i0 := Random (maxint);
  TMap64(Result).i1 := Random (maxint);
End;


Liefert eine Zufallszahl im Bereich 0... maxint64. Ich denke, mit einer schön langen Periode.

Moderiert von user profile iconraziel: Delphi-Tags hinzugefügt.


Spaceguide - Di 13.09.05 21:01

Ich brauche doch keine 64-Bit Ausgabe, sondern einen 64-Bit Seed. Bei deinem Konstrukt halbiert sich die Periode des Borland-Zufallsgenerators auch noch.


alzaimar - Di 13.09.05 21:12

Blöde Frage, wieso halbiert sich die Periode?
Wenn ich statt 64bit eine 2 bit zufallszahl (0..3) aus zwei 1-bit zufallszahlen (0,1) so konstruiere, bekomme ich doch die doppelte Periode, oder nicht? Na, jedenfalls halbiere ich die ohnehin extrem kurze Periode der 1-bit Zufallszahl nicht :wink:.

Ejal, was verstehst Du unter 64-bit seed, wenn Du keine 64 bit Randomfunktion willst? Es gibt ja standard Randomgeneratoren, wie Du sicherlich weisst. Die kannst Du doch einfach nachprogrammieren... (3-19 Zeilen)


Spaceguide - Di 13.09.05 21:44

Das ist doch ganz klar: Du hast die globale Variable RandSeed, die ist 32-Bit lang. Bei jedem Aufruf von random bekommt diese einen neuen Wert. Da sie nur 32-bit lang ist, wiederholen sich die Werte nach spätestens 2^32 Iterationen, meist früher. Da du in deiner Funktion random zweimal aufrufst, halbiert sich die Anzahl der Iterationen, sprich die Periode halbiert sich.

Schlechte Zufallsgeneratoren sind schnell gemacht, ich hätte aber gerne einen guten, wie z.B. den MersenneTwister (hat leider auch nur 32-bit Seed, arbeitet intern mit wesentlich mehr, um eine Periode von 2^19937 - 1 zu erreichen!). Der Delphi Quellcode des MersenneTwisters ist knapp 2700 Zeilen.


alzaimar - Di 13.09.05 23:08

user profile iconSpaceguide hat folgendes geschrieben:
..Da sie nur 32-bit lang ist, wiederholen sich die Werte nach spätestens 2^32 Iterationen, meist früher. Da du in deiner Funktion random zweimal aufrufst, halbiert sich die Anzahl der Iterationen, sprich die Periode halbiert sich.

Nö. Beispiel: Random mit periode 3.... 0,1,2,0,1,2... wieder mein beispiel mit 2 bits. random liefert also 0,1,0,0,1,0,0,1,0.
Das ergibt 01,00,10,10,01... also 1,0,2,2,1,1,0,2,2,1,1... Das ist doch eine Periode 5 (1,0,2,2,1) oder nicht? Ok, Die 3 kommt erst gar nicht vor, dieser Generator ist Mist, mein Vorschlag ist Schrott, aber halbieren tut er die Periode nur gdw. die Periode eine gerade Zahl ist.

So! Ehre gerettet, Problem nicht gelöst. Muss es denn unbedingt ein Monolith a la Mersenne Twister sein? Der Randomgenerator von Delphi ist doch in Ordnung...


delfiphan - Di 13.09.05 23:12

Delphi-Random-Generator veranschaulicht (Referenz [http://www.delphi-forum.de/viewtopic.php?p=251750#251750]):

Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
var  
 I: DWord;  
 MyRand: DWord;  
begin  
 randseed := 0;  
 for I := 1 to 10 do // 10 Random-Zahlen ausgeben  
  Memo1.Lines.Add(FloatToStr(random));  
 // äquivalent zu:  
 MyRand := 0;  
 for I := 1 to 10 do  
 begin  
  MyRand := MyRand*134775813+1;  
  Memo1.Lines.Add(FloatToStr(MyRand/4294967296.0));  
 end;  
end;


Für random(N) wird iirc Trunc(random*N) gerechnet. Für nen 64 bit Generator kannst du ja statt die Zahlen 134775813 und 4294967296 andere nehmen (aber nicht einfach irgendwelche! Such mal im Internet). Oder du nimmst grad nen Lagged Fibonacci Generator, oder sonst was...

Im Übrigen muss ich Spaceguide Recht geben: Die Periode ist bei deinem Random64 höchstens halb so lang wie das normale random, weil pro random64 zweimal random ausgeführt wird, und die Zahlen deswegen auch doppelt so schnell ausgehen.


Spaceguide - Mi 14.09.05 00:19

Was würde eigentlich passieren, wenn ich den Mersennetwister weiterhin verwende, ihn aber in meinen eigenen kapsele.

Ich habe dann mein int64 (gibt's eigentlich kein word64, also unsigned?) und zerlege das mal in zwei dwords SeedA und SeedB. Jetzt initialisiere ich zwei Mersennetwister A und B mit jeweils SeedA und SeedB. Die Randomfunktion macht dann folgendes: Result:=Frac(A.Random+B.Random). Der Mersennetwister hat wie gesagt eine Periode von 2^19937-1. Irgendwelche Mathematiker hier, die mehr als nur eine Vorlesung Stochastik hatten, wie ich? ;-)


delfiphan - Mi 14.09.05 12:56

Hmm. Ohne Frac wäre es eine Verteilungsfunktion, die wie ein Dreieck aussieht (dazu hab ich irgendwo mal was gepostet...). Mit dem Frac wird dann aber der rechte Teil des Dreiecks auf den linken Teil kopiert und es kommt meiner Meinung nach wieder eine gleichverteilte Zufallsfunktion raus.
Jedoch ohne Garantie. Kannst ja mal einen Test laufen lassen.

Unsigned 64 bit integer = uint64. Irgendwie ist das nirgends dokumentiert ;)


Spaceguide - Mi 14.09.05 13:07

Frac(A+B) funktioniert nicht, da bekomm ich doppelte Sequenzen. Frac(A+B*2) scheint gut zu funktionieren.


delfiphan - Mi 14.09.05 13:10

Doppelte Sequenzen?


Spaceguide - Mi 14.09.05 13:33

Ja, z.b. bei $1 und $100000000, da Seed64 = SeedA + SeedB shl 32;


SMO - Mi 14.09.05 13:37

user profile iconSpaceguide hat folgendes geschrieben:
Ich habe dann mein int64 (gibt's eigentlich kein word64, also unsigned?)

Scheinbar schon: UInt64. Dieser Typ ist aber undokumentiert und deswegen möglicherweise fehlerhaft/riskant. Habe selbst erst vor kurzem davon erfahren, siehe hier [http://www.delphipraxis.net/topic63238_delphis+largeworddivisionsfehler+beheben.html].


delfiphan - Mi 14.09.05 13:45

user profile iconSpaceguide hat folgendes geschrieben:
Ja, z.b. bei $1 und $100000000, da Seed64 = SeedA + SeedB shl 32;

$100000000 ist nicht mehr durch eine 32 bit Zahl darstellbar; nur bis und mit $FFFFFFFF. Oder was meinst du genau?
Edit: Ach so, du meinst du nimmst zwei identische Zufallsgeneratoren mit anderen Seeds? Würd ich aber nicht empfehlen.

@SMO: 5 Posts weiter oben ;)


Spaceguide - Mi 14.09.05 13:54

user profile icondelfiphan hat folgendes geschrieben:

$100000000 ist nicht mehr durch eine 32 bit Zahl darstellbar; nur bis und mit $FFFFFFFF. Oder was meinst du genau?


Man beachte das 64 hinter Seed64. Ich hab doch gesagt, dass ich einen 64-Bit Seed brauche, den in zwei 32-bit Seeds zerlege (höhrere und niedere Bits) und damit den Zufallsgenerator initialisiere. Verwende ich jedoch Frac(A+B) zum zusammensetzen, kommen bei $1 und $100000000 die gleichen Sequenzen raus? Warum? Das ist ganz einfach:

SeedA := Seed64 and $FFFFFFFF;
SeedB := Seed64 shr 32;

$1 => SeedA = 1, SeedB = 0
$100000000 => SeedA = 0, SeedA = 1

Dürfte jetzt trivial genug sein, um es zu verstehen.


alzaimar - Mi 14.09.05 14:32

user profile iconSpaceguide hat folgendes geschrieben:
...Der Delphi Quellcode des MersenneTwisters ist knapp 2700 Zeilen.

Zähl doch nochmal nach... von der homepage des MT

Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
unit MT19937;
{$R-} {range checking 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: matumoto@math.keio.ac.jp
   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-1of longint; // the array for the state vector

procedure sgenrand_MT19937(seed: longint);         // Initialization by seed
procedure lsgenrand_MT19937(const seed_array: tMT19937StateArray); // Initialization by array of seeds
procedure randomize_MT19937;                       // randomization
function  randInt_MT19937(Range: longint):longint; // integer RANDOM with positive range
function  genrand_MT19937: longint;                // random longint (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;


VAR
  mt : tMT19937StateArray;
  mti: integer=MT19937N+1// mti=MT19937N+1 means mt[] is not initialized

{ Initializing the array with a seed }
procedure sgenrand_MT19937(seed: longint);
var
  i: integer;
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 $ffff0000shr 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: integer;
begin
  for i := 0 to MT19937N-1 do mt[i] := seed_array[i];
  mti := MT19937N;
end;

function genrand_MT19937: longint;
const
  mag01 : array [0..1of longint =(0, MT19937MATRIX_A);
var
  y: longint;
  kk: integer;
begin
  if mti >= MT19937N { generate MT19937N longints at one time }
  then begin
     if mti = (MT19937N+1then  // 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+1and MT19937LOWER_MASK);
        mt[kk] := mt[kk+MT19937M] xor (y shr 1xor mag01[y and $00000001];
     end;
     for kk:= MT19937N-MT19937M to MT19937N-2 do begin
       y := (mt[kk] and MT19937UPPER_MASK) or (mt[kk+1and MT19937LOWER_MASK);
       mt[kk] := mt[kk+(MT19937M-MT19937N)] xor (y shr 1xor mag01[y and $00000001];
     end;
     y := (mt[MT19937N-1and MT19937UPPER_MASK) or (mt[0and MT19937LOWER_MASK);
     mt[MT19937N-1] := mt[MT19937M-1xor (y shr 1xor 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 15and TEMPERING_MASK_C;
  y := y xor (y shr 18);
  Result := y;
end;

{ Delphi interface }

procedure Randomize_MT19937;
Var OldRandSeed: longint;
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;

// bug fixed 21.6.2000. 
Function  RandInt_MT19937(Range: longint):longint;
// EAX <- Range
// Result -> EAX
asm
  PUSH  EAX
  CALL  genrand_MT19937
  POP   EDX
  MUL   EDX
  MOV   EAX,EDX
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;
end.


Spaceguide - Mi 14.09.05 14:53

2689 Zeilen, besteht auch hauptsächlich aus Assembler.


alzaimar - Mi 14.09.05 17:54

Wieso komm ich auf 202, bestehend haupsächlich auf Kommentaren und Leerzeilen?


Spaceguide - Mi 14.09.05 18:34

Sind wohl noch ein paar Sachen wie SHA in der Unit und echt jede Menge Kommentare.