1 |
// Upgraded to Delphi 2009: Sebastian Zierer
|
2 |
|
3 |
(* ***** BEGIN LICENSE BLOCK *****
|
4 |
* Version: MPL 1.1
|
5 |
*
|
6 |
* The contents of this file are subject to the Mozilla Public License Version
|
7 |
* 1.1 (the "License"); you may not use this file except in compliance with
|
8 |
* the License. You may obtain a copy of the License at
|
9 |
* http://www.mozilla.org/MPL/
|
10 |
*
|
11 |
* Software distributed under the License is distributed on an "AS IS" basis,
|
12 |
* WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
|
13 |
* for the specific language governing rights and limitations under the
|
14 |
* License.
|
15 |
*
|
16 |
* The Original Code is TurboPower SysTools
|
17 |
*
|
18 |
* The Initial Developer of the Original Code is
|
19 |
* TurboPower Software
|
20 |
*
|
21 |
* Portions created by the Initial Developer are Copyright (C) 1996-2002
|
22 |
* the Initial Developer. All Rights Reserved.
|
23 |
*
|
24 |
* Contributor(s):
|
25 |
*
|
26 |
* ***** END LICENSE BLOCK ***** *)
|
27 |
|
28 |
{*********************************************************}
|
29 |
{* SysTools: StRandom.pas 4.04 *}
|
30 |
{*********************************************************}
|
31 |
{* SysTools: Classes for random number distributions *}
|
32 |
{*********************************************************}
|
33 |
|
34 |
{$I StDefine.inc}
|
35 |
|
36 |
unit StRandom;
|
37 |
|
38 |
interface
|
39 |
|
40 |
uses
|
41 |
Windows, SysUtils, Classes,
|
42 |
StBase;
|
43 |
|
44 |
type
|
45 |
TStRandomBase = class
|
46 |
private
|
47 |
protected
|
48 |
function rbMarsagliaGamma(aShape : double) : double;
|
49 |
function rbMontyPythonNormal : double;
|
50 |
public
|
51 |
{uniform distributions}
|
52 |
function AsFloat : double; virtual; abstract;
|
53 |
function AsInt(aUpperLimit : integer) : integer;
|
54 |
function AsIntInRange(aLowerLimit : integer;
|
55 |
aUpperLimit : integer) : integer;
|
56 |
|
57 |
{continuous non-uniform distributions}
|
58 |
function AsBeta(aShape1, aShape2 : double) : double;
|
59 |
function AsCauchy : double;
|
60 |
function AsChiSquared(aFreedom : integer) : double;
|
61 |
function AsErlang(aMean : double;
|
62 |
aOrder : integer) : double;
|
63 |
function AsExponential(aMean : double) : double;
|
64 |
function AsF(aFreedom1 : integer;
|
65 |
aFreedom2 : integer) : double;
|
66 |
function AsGamma(aShape : double; aScale : double) : double;
|
67 |
function AsLogNormal(aMean : double;
|
68 |
aStdDev : double) : double;
|
69 |
function AsNormal(aMean : double;
|
70 |
aStdDev : double) : double;
|
71 |
function AsT(aFreedom : integer) : double;
|
72 |
function AsWeibull(aShape : double;
|
73 |
aScale : double) : double;
|
74 |
end;
|
75 |
|
76 |
TStRandomSystem = class(TStRandomBase)
|
77 |
private
|
78 |
FSeed : integer;
|
79 |
protected
|
80 |
procedure rsSetSeed(aValue : integer);
|
81 |
public
|
82 |
constructor Create(aSeed : integer);
|
83 |
function AsFloat : double; override;
|
84 |
property Seed : integer read FSeed write rsSetSeed;
|
85 |
end;
|
86 |
|
87 |
TStRandomCombined = class(TStRandomBase)
|
88 |
private
|
89 |
FSeed1 : integer;
|
90 |
FSeed2 : integer;
|
91 |
protected
|
92 |
procedure rcSetSeed1(aValue : integer);
|
93 |
procedure rcSetSeed2(aValue : integer);
|
94 |
public
|
95 |
constructor Create(aSeed1, aSeed2 : integer);
|
96 |
function AsFloat : double; override;
|
97 |
property Seed1 : integer read FSeed1 write rcSetSeed1;
|
98 |
property Seed2 : integer read FSeed2 write rcSetSeed2;
|
99 |
end;
|
100 |
|
101 |
TStRandomMother = class(TStRandomBase)
|
102 |
private
|
103 |
FNminus4 : integer;
|
104 |
FNminus3 : integer;
|
105 |
FNminus2 : integer;
|
106 |
FNminus1 : integer;
|
107 |
FC : integer;
|
108 |
protected
|
109 |
procedure rsSetSeed(aValue : integer);
|
110 |
public
|
111 |
constructor Create(aSeed : integer);
|
112 |
function AsFloat : double; override;
|
113 |
property Seed : integer write rsSetSeed;
|
114 |
end;
|
115 |
|
116 |
implementation
|
117 |
|
118 |
uses
|
119 |
StConst;
|
120 |
|
121 |
var
|
122 |
Root2Pi : double;
|
123 |
InvRoot2Pi : double;
|
124 |
RootLn4 : double;
|
125 |
Ln2 : double;
|
126 |
MPN_s : double;
|
127 |
Ln2MPN_s : double;
|
128 |
MPN_sPlus1 : double;
|
129 |
|
130 |
Mum1 : integer;
|
131 |
Mum2 : integer;
|
132 |
Mum3 : integer;
|
133 |
Mum4 : integer;
|
134 |
|
135 |
{===Helper routines==================================================}
|
136 |
function GetRandomSeed : integer;
|
137 |
var
|
138 |
Hash : integer;
|
139 |
SystemTime: TSystemTime;
|
140 |
G : integer;
|
141 |
begin
|
142 |
{start with the tick count}
|
143 |
Hash := integer(GetTickCount);
|
144 |
|
145 |
{get the current time}
|
146 |
GetLocalTime(SystemTime);
|
147 |
|
148 |
{hash in the milliseconds}
|
149 |
Hash := (Hash shl 4) + SystemTime.wMilliseconds;
|
150 |
G := Hash and longint($F0000000);
|
151 |
if (G <> 0) then
|
152 |
Hash := (Hash xor (G shr 24)) xor G;
|
153 |
|
154 |
{hash in the second}
|
155 |
Hash := (Hash shl 4) + SystemTime.wSecond;
|
156 |
G := Hash and longint($F0000000);
|
157 |
if (G <> 0) then
|
158 |
Hash := (Hash xor (G shr 24)) xor G;
|
159 |
|
160 |
{hash in the minute}
|
161 |
Hash := (Hash shl 4) + SystemTime.wMinute;
|
162 |
G := Hash and longint($F0000000);
|
163 |
if (G <> 0) then
|
164 |
Hash := (Hash xor (G shr 24)) xor G;
|
165 |
|
166 |
{hash in the hour}
|
167 |
Hash := (Hash shl 3) + SystemTime.wHour;
|
168 |
G := Hash and longint($F0000000);
|
169 |
if (G <> 0) then
|
170 |
Hash := (Hash xor (G shr 24)) xor G;
|
171 |
|
172 |
{return the hash}
|
173 |
Result := Hash;
|
174 |
end;
|
175 |
{====================================================================}
|
176 |
|
177 |
|
178 |
{===TStRandomBase====================================================}
|
179 |
function TStRandomBase.AsBeta(aShape1, aShape2 : double) : double;
|
180 |
var
|
181 |
R1, R2 : double;
|
182 |
begin
|
183 |
if not ((aShape1 > 0.0) and (aShape2 > 0.0)) then
|
184 |
raise EStPRNGError.Create(stscPRNGBetaShapeS);
|
185 |
|
186 |
if (aShape2 = 1.0) then begin
|
187 |
repeat
|
188 |
R1 := AsFloat;
|
189 |
until R1 <> 0.0;
|
190 |
Result := exp(ln(R1) / aShape1);
|
191 |
end
|
192 |
else if (aShape1 = 1.0) then begin
|
193 |
repeat
|
194 |
R1 := AsFloat;
|
195 |
until R1 <> 0.0;
|
196 |
Result := 1.0 - exp(ln(R1) / aShape1);
|
197 |
end
|
198 |
else begin
|
199 |
R1 := AsGamma(aShape1, 1.0);
|
200 |
R2 := AsGamma(aShape2, 1.0);
|
201 |
Result := R1 / (R1 + R2);
|
202 |
end;
|
203 |
end;
|
204 |
{--------}
|
205 |
function TStRandomBase.AsCauchy : double;
|
206 |
var
|
207 |
x : double;
|
208 |
y : double;
|
209 |
begin
|
210 |
repeat
|
211 |
repeat
|
212 |
x := AsFloat;
|
213 |
until (x <> 0.0);
|
214 |
y := (AsFloat * 2.0) - 1.0;
|
215 |
until sqr(x) + sqr(y) < 1.0;
|
216 |
Result := y / x;
|
217 |
end;
|
218 |
{--------}
|
219 |
function TStRandomBase.AsChiSquared(aFreedom : integer) : double;
|
220 |
begin
|
221 |
if not (aFreedom > 0) then
|
222 |
raise EStPRNGError.Create(stscPRNGDegFreedomS);
|
223 |
|
224 |
Result := AsGamma(aFreedom * 0.5, 2.0)
|
225 |
end;
|
226 |
{--------}
|
227 |
function TStRandomBase.AsErlang(aMean : double;
|
228 |
aOrder : integer) : double;
|
229 |
var
|
230 |
Product : double;
|
231 |
i : integer;
|
232 |
begin
|
233 |
if not (aMean > 0.0) then
|
234 |
raise EStPRNGError.Create(stscPRNGMeanS);
|
235 |
if not (aOrder > 0) then
|
236 |
raise EStPRNGError.Create(stscPRNGErlangOrderS);
|
237 |
|
238 |
if (aOrder < 10) then begin
|
239 |
Product := 1.0;
|
240 |
for i := 1 to aOrder do
|
241 |
Product := Product * AsFloat;
|
242 |
Result := -aMean * ln(Product) / aOrder;
|
243 |
end
|
244 |
else begin
|
245 |
Result := AsGamma(aOrder, aMean);
|
246 |
end;
|
247 |
end;
|
248 |
{--------}
|
249 |
function TStRandomBase.AsExponential(aMean : double) : double;
|
250 |
var
|
251 |
R : double;
|
252 |
begin
|
253 |
if not (aMean > 0.0) then
|
254 |
raise EStPRNGError.Create(stscPRNGMeanS);
|
255 |
|
256 |
repeat
|
257 |
R := AsFloat;
|
258 |
until (R <> 0.0);
|
259 |
Result := -aMean * ln(R);
|
260 |
end;
|
261 |
{--------}
|
262 |
function TStRandomBase.AsF(aFreedom1 : integer;
|
263 |
aFreedom2 : integer) : double;
|
264 |
begin
|
265 |
Result := (AsChiSquared(aFreedom1) * aFreedom1) /
|
266 |
(AsChiSquared(aFreedom2) * aFreedom2);
|
267 |
end;
|
268 |
{--------}
|
269 |
function TStRandomBase.AsGamma(aShape : double; aScale : double) : double;
|
270 |
var
|
271 |
R : double;
|
272 |
begin
|
273 |
if not (aShape > 0.0) then
|
274 |
raise EStPRNGError.Create(stscPRNGGammaShapeS);
|
275 |
if not (aScale > 0.0) then
|
276 |
raise EStPRNGError.Create(stscPRNGGammaScaleS);
|
277 |
|
278 |
{there are three cases:
|
279 |
..0.0 < shape < 1.0, use Marsaglia's technique of
|
280 |
Gamma(shape) = Gamma(shape+1) * uniform^(1/shape)}
|
281 |
if (aShape < 1.0) then begin
|
282 |
repeat
|
283 |
R := AsFloat;
|
284 |
until (R <> 0.0);
|
285 |
Result := aScale * rbMarsagliaGamma(aShape + 1.0) *
|
286 |
exp(ln(R) / aShape);
|
287 |
end
|
288 |
|
289 |
{..shape = 1.0: this is the same as exponential}
|
290 |
else if (aShape = 1.0) then begin
|
291 |
repeat
|
292 |
R := AsFloat;
|
293 |
until (R <> 0.0);
|
294 |
Result := aScale * -ln(R);
|
295 |
end
|
296 |
|
297 |
{..shape > 1.0: use Marsaglia./Tsang algorithm}
|
298 |
else begin
|
299 |
Result := aScale * rbMarsagliaGamma(aShape);
|
300 |
end;
|
301 |
end;
|
302 |
{--------}
|
303 |
function TStRandomBase.AsInt(aUpperLimit : integer) : integer;
|
304 |
begin
|
305 |
if not (aUpperLimit > 0) then
|
306 |
raise EStPRNGError.Create(stscPRNGLimitS);
|
307 |
|
308 |
Result := Trunc(AsFloat * aUpperLimit);
|
309 |
end;
|
310 |
{--------}
|
311 |
function TStRandomBase.AsIntInRange(aLowerLimit : integer;
|
312 |
aUpperLimit : integer) : integer;
|
313 |
begin
|
314 |
if not (aLowerLimit < aUpperLimit) then
|
315 |
raise EStPRNGError.Create(stscPRNGUpperLimitS);
|
316 |
|
317 |
Result := Trunc(AsFloat * (aUpperLimit - aLowerLimit)) + ALowerLimit;
|
318 |
end;
|
319 |
{--------}
|
320 |
function TStRandomBase.AsLogNormal(aMean : double;
|
321 |
aStdDev : double) : double;
|
322 |
begin
|
323 |
Result := exp(AsNormal(aMean, aStdDev));
|
324 |
end;
|
325 |
{--------}
|
326 |
function TStRandomBase.AsNormal(aMean : double;
|
327 |
aStdDev : double) : double;
|
328 |
begin
|
329 |
if not (aStdDev > 0.0) then
|
330 |
raise EStPRNGError.Create(stscPRNGStdDevS);
|
331 |
|
332 |
Result := (rbMontyPythonNormal * aStdDev) + aMean;
|
333 |
|
334 |
(*** alternative: The Box-Muller transformation
|
335 |
var
|
336 |
R1, R2 : double;
|
337 |
RadiusSqrd : double;
|
338 |
begin
|
339 |
{get two random numbers that define a point in the unit circle}
|
340 |
repeat
|
341 |
R1 := (2.0 * aRandGen.AsFloat) - 1.0;
|
342 |
R2 := (2.0 * aRandGen.AsFloat) - 1.0;
|
343 |
RadiusSqrd := sqr(R1) + sqr(R2);
|
344 |
until (RadiusSqrd < 1.0) and (RadiusSqrd > 0.0);
|
345 |
|
346 |
{apply Box-Muller transformation}
|
347 |
Result := (R1 * sqrt(-2.0 * ln(RadiusSqrd) / RadiusSqrd) * aStdDev)
|
348 |
+ aMean;
|
349 |
***)
|
350 |
end;
|
351 |
{--------}
|
352 |
function TStRandomBase.AsT(aFreedom : integer) : double;
|
353 |
begin
|
354 |
if not (aFreedom > 0) then
|
355 |
raise EStPRNGError.Create(stscPRNGDegFreedomS);
|
356 |
|
357 |
Result := rbMontyPythonNormal / sqrt(AsChiSquared(aFreedom) / aFreedom);
|
358 |
end;
|
359 |
{--------}
|
360 |
function TStRandomBase.AsWeibull(aShape : double;
|
361 |
aScale : double) : double;
|
362 |
var
|
363 |
R : double;
|
364 |
begin
|
365 |
if not (aShape > 0) then
|
366 |
raise EStPRNGError.Create(stscPRNGWeibullShapeS);
|
367 |
if not (aScale > 0) then
|
368 |
raise EStPRNGError.Create(stscPRNGWeibullScaleS);
|
369 |
|
370 |
repeat
|
371 |
R := AsFloat;
|
372 |
until (R <> 0.0);
|
373 |
Result := exp(ln(-ln(R)) / aShape) * aScale;
|
374 |
end;
|
375 |
{--------}
|
376 |
|
377 |
function TStRandomBase.rbMarsagliaGamma(aShape : double) : double;
|
378 |
var
|
379 |
d : double;
|
380 |
c : double;
|
381 |
x : double;
|
382 |
v : double;
|
383 |
u : double;
|
384 |
Done : boolean;
|
385 |
begin
|
386 |
{Notes: implements the Marsaglia/Tsang method of generating random
|
387 |
numbers belonging to the gamma distribution:
|
388 |
|
389 |
Marsaglia & Tsang, "A Simple Method for Generating Gamma
|
390 |
Variables", ACM Transactions on Mathematical Software,
|
391 |
Vol. 26, No. 3, September 2000, Pages 363-372
|
392 |
|
393 |
It is pointless to try and work out what's going on in this
|
394 |
routine without reading this paper :-)
|
395 |
}
|
396 |
|
397 |
d := aShape - (1.0 / 3.0);
|
398 |
c := 1.0 / sqrt(9.0 * d);
|
399 |
Done := false;
|
400 |
{$IFDEF SuppressWarnings}
|
401 |
v := 0.0;
|
402 |
{$ENDIF}
|
403 |
|
404 |
while not Done do begin
|
405 |
repeat
|
406 |
x := rbMontyPythonNormal;
|
407 |
v := 1.0 + (c * x);
|
408 |
until (v > 0.0);
|
409 |
|
410 |
v := v * v * v;
|
411 |
u := AsFloat;
|
412 |
|
413 |
Done := u < (1.0 - 0.0331 * sqr(sqr(x)));
|
414 |
|
415 |
if not Done then
|
416 |
Done := ln(u) < (0.5 * sqr(x)) + d * (1.0 - v + ln(v))
|
417 |
end;
|
418 |
|
419 |
Result := d * v;
|
420 |
end;
|
421 |
{--------}
|
422 |
function TStRandomBase.rbMontyPythonNormal : double;
|
423 |
var
|
424 |
x : double;
|
425 |
y : double;
|
426 |
v : double;
|
427 |
NonZeroRandom : double;
|
428 |
begin
|
429 |
{Notes: implements the Monty Python method of generating random
|
430 |
numbers belonging to the Normal (Gaussian) distribution:
|
431 |
|
432 |
Marsaglia & Tsang, "The Monty Python Method for Generating
|
433 |
Random Variables", ACM Transactions on Mathematical
|
434 |
Software, Vol. 24, No. 3, September 1998, Pages 341-350
|
435 |
|
436 |
It is pointless to try and work out what's going on in this
|
437 |
routine without reading this paper :-)
|
438 |
|
439 |
Some constants:
|
440 |
a = sqrt(ln(4))
|
441 |
b = sqrt(2 * pi)
|
442 |
s = a / (b - a)
|
443 |
}
|
444 |
|
445 |
{step 1: generate a random number x between +/- sqrt(2*Pi) and
|
446 |
return it if its absolute value is less than sqrt(ln(4));
|
447 |
note that this exit will happen about 47% of the time}
|
448 |
x := ((AsFloat * 2.0) - 1.0) * Root2Pi;
|
449 |
if (abs(x) < RootLn4) then begin
|
450 |
Result := x;
|
451 |
Exit;
|
452 |
end;
|
453 |
|
454 |
{step 2a: generate another random number y strictly between 0 and 1}
|
455 |
repeat
|
456 |
y := AsFloat;
|
457 |
until (y <> 0.0);
|
458 |
|
459 |
{step 2b: the first quadratic pretest avoids ln() calculation
|
460 |
calculate v = 2.8658 - |x| * (2.0213 - 0.3605 * |x|)
|
461 |
return x if y < v}
|
462 |
v := 2.8658 - Abs(x) * (2.0213 - 0.3605 * Abs(x));
|
463 |
if (y < v) then begin
|
464 |
Result := x;
|
465 |
Exit;
|
466 |
end;
|
467 |
|
468 |
{step 2c: the second quadratic pretest again avoids ln() calculation
|
469 |
return s * (b - x) if y > v + 0.0506}
|
470 |
if (y > v + 0.0506) then begin
|
471 |
if (x > 0) then
|
472 |
Result := MPN_s * (Root2Pi - x)
|
473 |
else
|
474 |
Result := -MPN_s * (Root2Pi + x);
|
475 |
Exit;
|
476 |
end;
|
477 |
|
478 |
{step 2d: return x if y < f(x) or
|
479 |
ln(y) < ln(2) - (0.5 * x * x) }
|
480 |
if (ln(y) < (Ln2 - (0.5 * x * x))) then begin
|
481 |
Result := x;
|
482 |
Exit;
|
483 |
end;
|
484 |
|
485 |
{step 3: translate x to s * (b - x) and return it if y > g(x) or
|
486 |
ln(1 + s - y) < ln(2 * s) - (0.5 * x * x) }
|
487 |
if (x > 0) then
|
488 |
x := MPN_s * (Root2Pi - x)
|
489 |
else
|
490 |
x := -MPN_s * (Root2Pi + x);
|
491 |
if (ln(MPN_sPlus1 - y) < (Ln2MPN_s - (0.5 * x * x))) then begin
|
492 |
Result := x;
|
493 |
Exit;
|
494 |
end;
|
495 |
|
496 |
{step 4: the iterative process}
|
497 |
repeat
|
498 |
repeat
|
499 |
NonZeroRandom := AsFloat;
|
500 |
until (NonZeroRandom <> 0.0);
|
501 |
x := -ln(NonZeroRandom) * InvRoot2Pi;
|
502 |
repeat
|
503 |
NonZeroRandom := AsFloat;
|
504 |
until (NonZeroRandom <> 0.0);
|
505 |
y := -ln(NonZeroRandom);
|
506 |
until (y + y) > (x * x);
|
507 |
if (NonZeroRandom < 0.5) then
|
508 |
Result := -(Root2Pi + x)
|
509 |
else
|
510 |
Result := Root2Pi + x;
|
511 |
end;
|
512 |
{====================================================================}
|
513 |
|
514 |
|
515 |
{===TStRandomSystem==================================================}
|
516 |
constructor TStRandomSystem.Create(aSeed : integer);
|
517 |
begin
|
518 |
inherited Create;
|
519 |
Seed := aSeed;
|
520 |
end;
|
521 |
{--------}
|
522 |
function TStRandomSystem.AsFloat : double;
|
523 |
var
|
524 |
SaveSeed : integer;
|
525 |
begin
|
526 |
SaveSeed := RandSeed;
|
527 |
RandSeed := FSeed;
|
528 |
Result := System.Random;
|
529 |
FSeed := RandSeed;
|
530 |
RandSeed := SaveSeed;
|
531 |
end;
|
532 |
{--------}
|
533 |
procedure TStRandomSystem.rsSetSeed(aValue : integer);
|
534 |
begin
|
535 |
if (aValue = 0) then
|
536 |
FSeed := GetRandomSeed
|
537 |
else
|
538 |
FSeed := aValue;
|
539 |
end;
|
540 |
{====================================================================}
|
541 |
|
542 |
|
543 |
{===TStRandomCombined================================================}
|
544 |
const
|
545 |
m1 = 2147483563;
|
546 |
m2 = 2147483399;
|
547 |
{--------}
|
548 |
constructor TStRandomCombined.Create(aSeed1, aSeed2 : integer);
|
549 |
begin
|
550 |
inherited Create;
|
551 |
Seed1 := aSeed1;
|
552 |
if (aSeed1 = 0) and (aSeed2 = 0) then
|
553 |
Sleep(10); // a small delay to enable seed to change
|
554 |
Seed2 := aSeed2;
|
555 |
end;
|
556 |
{--------}
|
557 |
function TStRandomCombined.AsFloat : double;
|
558 |
const
|
559 |
a1 = 40014;
|
560 |
q1 = 53668; {equals m1 div a1}
|
561 |
r1 = 12211; {equals m1 mod a1}
|
562 |
|
563 |
a2 = 40692;
|
564 |
q2 = 52774; {equals m2 div a2}
|
565 |
r2 = 3791; {equals m2 mod a2}
|
566 |
|
567 |
OneOverM1 : double = 1.0 / m1;
|
568 |
var
|
569 |
k : longint;
|
570 |
Z : longint;
|
571 |
begin
|
572 |
{advance first PRNG}
|
573 |
k := FSeed1 div q1;
|
574 |
FSeed1 := (a1 * (FSeed1 - (k * q1))) - (k * r1);
|
575 |
if (FSeed1 < 0) then
|
576 |
inc(FSeed1, m1);
|
577 |
|
578 |
{advance second PRNG}
|
579 |
k := FSeed2 div q2;
|
580 |
FSeed2 := (a2 * (FSeed2 - (k * q2))) - (k * r2);
|
581 |
if (FSeed2 < 0) then
|
582 |
inc(FSeed2, m2);
|
583 |
|
584 |
{combine the two seeds}
|
585 |
Z := FSeed1 - FSeed2;
|
586 |
if (Z <= 0) then
|
587 |
Z := Z + m1 - 1;
|
588 |
Result := Z * OneOverM1;
|
589 |
end;
|
590 |
{--------}
|
591 |
procedure TStRandomCombined.rcSetSeed1(aValue : integer);
|
592 |
begin
|
593 |
if (aValue = 0) then
|
594 |
FSeed1 := GetRandomSeed
|
595 |
else
|
596 |
FSeed1 := aValue;
|
597 |
end;
|
598 |
{--------}
|
599 |
procedure TStRandomCombined.rcSetSeed2(aValue : integer);
|
600 |
begin
|
601 |
if (aValue = 0) then
|
602 |
FSeed2 := GetRandomSeed
|
603 |
else
|
604 |
FSeed2 := aValue;
|
605 |
end;
|
606 |
{====================================================================}
|
607 |
|
608 |
|
609 |
{===TStRandomMother==================================================}
|
610 |
constructor TStRandomMother.Create(aSeed : integer);
|
611 |
begin
|
612 |
inherited Create;
|
613 |
Seed := aSeed;
|
614 |
end;
|
615 |
{--------}
|
616 |
function TStRandomMother.AsFloat : double;
|
617 |
const
|
618 |
TwoM31 : double = 1.0 / $7FFFFFFF;
|
619 |
begin
|
620 |
asm
|
621 |
push esi
|
622 |
push edi
|
623 |
push ebx
|
624 |
|
625 |
{get around a compiler bug where it doesn't notice that edx is
|
626 |
being changed in the asm code !!! D5 bug}
|
627 |
push edx
|
628 |
|
629 |
{set ebx to point to self}
|
630 |
mov ebx, eax
|
631 |
|
632 |
{multiply X(n-4) by 21111111}
|
633 |
mov eax, [ebx].TStRandomMother.FNMinus4
|
634 |
mul [Mum1]
|
635 |
mov edi, eax
|
636 |
mov esi, edx
|
637 |
|
638 |
{multiply X(n-3) by 1492 (save it in X(n-4) before though)}
|
639 |
mov eax, [ebx].TStRandomMother.FNMinus3
|
640 |
mov [ebx].TStRandomMother.FNMinus4, eax
|
641 |
mul [Mum2]
|
642 |
add edi, eax
|
643 |
adc esi, edx
|
644 |
|
645 |
{multiply X(n-2) by 1776 (save it in X(n-3) before though)}
|
646 |
mov eax, [ebx].TStRandomMother.FNMinus2
|
647 |
mov [ebx].TStRandomMother.FNMinus3, eax
|
648 |
mul [Mum3]
|
649 |
add edi, eax
|
650 |
adc esi, edx
|
651 |
|
652 |
{multiply X(n-1) by 5115 (save it in X(n-2) before though)}
|
653 |
mov eax, [ebx].TStRandomMother.FNMinus1
|
654 |
mov [ebx].TStRandomMother.FNMinus2, eax
|
655 |
mul [Mum4]
|
656 |
add edi, eax
|
657 |
adc esi, edx
|
658 |
|
659 |
{add in the remainder}
|
660 |
add edi, [ebx].TStRandomMother.FC
|
661 |
adc esi, 0;
|
662 |
|
663 |
{save the lower 32 bits in X(n-1), the upper into the remainder}
|
664 |
mov [ebx].TStRandomMother.FNMinus1, edi
|
665 |
mov [ebx].TStRandomMother.FC, esi
|
666 |
|
667 |
{get around a compiler bug where it doesn't notice that edx was
|
668 |
changed in the asm code !!! D5 bug}
|
669 |
pop edx
|
670 |
|
671 |
pop ebx
|
672 |
pop edi
|
673 |
pop esi
|
674 |
end;
|
675 |
Result := (FNMinus1 shr 1) * TwoM31;
|
676 |
end;
|
677 |
{--------}
|
678 |
{$IFOPT Q+}
|
679 |
{note: TStRandomMother.rsSetSeed expressly overflows integers (it's
|
680 |
equivalent to calculating mod 2^32), so we have to force
|
681 |
overflow checks off}
|
682 |
{$DEFINE SaveQPlus}
|
683 |
{$Q-}
|
684 |
{$ENDIF}
|
685 |
procedure TStRandomMother.rsSetSeed(aValue : integer);
|
686 |
begin
|
687 |
if (aValue = 0) then
|
688 |
aValue := GetRandomSeed;
|
689 |
FNminus4 := aValue;
|
690 |
{note: the following code uses the generator
|
691 |
Xn := (69069 * Xn-1) mod 2^32
|
692 |
from D.E.Knuth, The Art of Computer Programming, Vol. 2
|
693 |
(second edition), Addison-Wesley, 1981, pp.102}
|
694 |
FNminus3 := 69069 * FNminus4;
|
695 |
FNminus2 := 69069 * FNminus3;
|
696 |
FNminus1 := 69069 * FNminus2;
|
697 |
FC := 69069 * FNminus1;
|
698 |
end;
|
699 |
{$IFDEF SaveQPlus}
|
700 |
{$Q+}
|
701 |
{$ENDIF}
|
702 |
{====================================================================}
|
703 |
|
704 |
|
705 |
{====================================================================}
|
706 |
procedure CalcConstants;
|
707 |
begin
|
708 |
{for the normal variates}
|
709 |
Root2Pi := sqrt(2 * Pi);
|
710 |
InvRoot2Pi := 1.0 / Root2Pi;
|
711 |
RootLn4 := sqrt(ln(4.0));
|
712 |
Ln2 := ln(2.0);
|
713 |
MPN_s := RootLn4 / (Root2Pi - RootLn4);
|
714 |
Ln2MPN_s := ln(2.0 * MPN_s);
|
715 |
MPN_sPlus1 := MPN_s + 1.0;
|
716 |
|
717 |
Mum1 := 2111111111;
|
718 |
Mum2 := 1492;
|
719 |
Mum3 := 1776;
|
720 |
Mum4 := 5115;
|
721 |
end;
|
722 |
{====================================================================}
|
723 |
|
724 |
|
725 |
initialization
|
726 |
CalcConstants;
|
727 |
|
728 |
end.
|