1 |
torben |
2671 |
// 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.
|