tpsystools_4.04/source/StRandom.pas

// Upgraded to Delphi 2009: Sebastian Zierer

(* ***** BEGIN LICENSE BLOCK *****
 * Version: MPL 1.1
 *
 * The contents of this file are subject to the Mozilla Public License Version
 * 1.1 (the "License"); you may not use this file except in compliance with
 * the License. You may obtain a copy of the License at
 * http://www.mozilla.org/MPL/
 *
 * Software distributed under the License is distributed on an "AS IS" basis,
 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
 * for the specific language governing rights and limitations under the
 * License.
 *
 * The Original Code is TurboPower SysTools
 *
 * The Initial Developer of the Original Code is
 * TurboPower Software
 *
 * Portions created by the Initial Developer are Copyright (C) 1996-2002
 * the Initial Developer. All Rights Reserved.
 *
 * Contributor(s):
 *
 * ***** END LICENSE BLOCK ***** *)

{*********************************************************}
{* SysTools: StRandom.pas 4.04                           *}
{*********************************************************}
{* SysTools: Classes for random number distributions     *}
{*********************************************************}

{$I StDefine.inc}

unit StRandom;

interface

uses
  Windows, SysUtils, Classes,
  StBase;

type
  TStRandomBase = class
    private
    protected
      function rbMarsagliaGamma(aShape   : double) : double;
      function rbMontyPythonNormal : double;
    public
      {uniform distributions}
      function AsFloat : double; virtual; abstract;
      function AsInt(aUpperLimit : integer) : integer;
      function AsIntInRange(aLowerLimit : integer;
                            aUpperLimit : integer) : integer;

      {continuous non-uniform distributions}
      function AsBeta(aShape1, aShape2 : double) : double;
      function AsCauchy : double;
      function AsChiSquared(aFreedom : integer) : double;
      function AsErlang(aMean  : double;
                        aOrder : integer) : double;
      function AsExponential(aMean : double) : double;
      function AsF(aFreedom1 : integer;
                   aFreedom2 : integer) : double;
      function AsGamma(aShape : double; aScale : double) : double;
      function AsLogNormal(aMean   : double;
                           aStdDev : double) : double;
      function AsNormal(aMean   : double;
                        aStdDev : double) : double;
      function AsT(aFreedom : integer) : double;
      function AsWeibull(aShape : double;
                         aScale : double) : double;
  end;

  TStRandomSystem = class(TStRandomBase)
    private
      FSeed : integer;
    protected
      procedure rsSetSeed(aValue : integer);
    public
      constructor Create(aSeed : integer);
      function AsFloat : double; override;
      property Seed : integer read FSeed write rsSetSeed;
  end;

  TStRandomCombined = class(TStRandomBase)
    private
      FSeed1 : integer;
      FSeed2 : integer;
    protected
      procedure rcSetSeed1(aValue : integer);
      procedure rcSetSeed2(aValue : integer);
    public
      constructor Create(aSeed1, aSeed2 : integer);
      function AsFloat : double; override;
      property Seed1 : integer read FSeed1 write rcSetSeed1;
      property Seed2 : integer read FSeed2 write rcSetSeed2;
  end;

  TStRandomMother = class(TStRandomBase)
    private
      FNminus4 : integer;
      FNminus3 : integer;
      FNminus2 : integer;
      FNminus1 : integer;
      FC       : integer;
    protected
      procedure rsSetSeed(aValue : integer);
    public
      constructor Create(aSeed : integer);
      function AsFloat : double; override;
      property Seed : integer write rsSetSeed;
  end;

implementation

uses
  StConst;

var
  Root2Pi    : double;
  InvRoot2Pi : double;
  RootLn4    : double;
  Ln2        : double;
  MPN_s      : double;
  Ln2MPN_s   : double;
  MPN_sPlus1 : double;

  Mum1       : integer;
  Mum2       : integer;
  Mum3       : integer;
  Mum4       : integer;

{===Helper routines==================================================}
function GetRandomSeed : integer;
var
  Hash : integer;
  SystemTime: TSystemTime;
  G : integer;
begin
  {start with the tick count}
  Hash := integer(GetTickCount);

  {get the current time}
  GetLocalTime(SystemTime);

  {hash in the milliseconds}
  Hash := (Hash shl 4) + SystemTime.wMilliseconds;
  G := Hash and longint($F0000000);
  if (G <> 0) then
    Hash := (Hash xor (G shr 24)) xor G;

  {hash in the second}
  Hash := (Hash shl 4) + SystemTime.wSecond;
  G := Hash and longint($F0000000);
  if (G <> 0) then
    Hash := (Hash xor (G shr 24)) xor G;

  {hash in the minute}
  Hash := (Hash shl 4) + SystemTime.wMinute;
  G := Hash and longint($F0000000);
  if (G <> 0) then
    Hash := (Hash xor (G shr 24)) xor G;

  {hash in the hour}
  Hash := (Hash shl 3) + SystemTime.wHour;
  G := Hash and longint($F0000000);
  if (G <> 0) then
    Hash := (Hash xor (G shr 24)) xor G;

  {return the hash}
  Result := Hash;
end;
{====================================================================}


{===TStRandomBase====================================================}
function TStRandomBase.AsBeta(aShape1, aShape2 : double) : double;
var
  R1, R2 : double;
begin
  if not ((aShape1 > 0.0) and (aShape2 > 0.0)) then
    raise EStPRNGError.Create(stscPRNGBetaShapeS);

  if (aShape2 = 1.0) then begin
    repeat
      R1 := AsFloat;
    until R1 <> 0.0;
    Result := exp(ln(R1) / aShape1);
  end
  else if (aShape1 = 1.0) then begin
    repeat
      R1 := AsFloat;
    until R1 <> 0.0;
    Result := 1.0 - exp(ln(R1) / aShape1);
  end
  else begin
    R1 := AsGamma(aShape1, 1.0);
    R2 := AsGamma(aShape2, 1.0);
    Result := R1 / (R1 + R2);
  end;
end;
{--------}
function TStRandomBase.AsCauchy : double;
var
  x : double;
  y : double;
begin
  repeat
    repeat
      x := AsFloat;
    until (x <> 0.0);
    y := (AsFloat * 2.0) - 1.0;
  until sqr(x) + sqr(y) < 1.0;
  Result := y / x;
end;
{--------}
function TStRandomBase.AsChiSquared(aFreedom : integer) : double;
begin
  if not (aFreedom > 0) then
    raise EStPRNGError.Create(stscPRNGDegFreedomS);

  Result := AsGamma(aFreedom * 0.5, 2.0)
end;
{--------}
function TStRandomBase.AsErlang(aMean  : double;
                                aOrder : integer) : double;
var
  Product : double;
  i       : integer;
begin
  if not (aMean > 0.0) then
    raise EStPRNGError.Create(stscPRNGMeanS);
  if not (aOrder > 0) then
    raise EStPRNGError.Create(stscPRNGErlangOrderS);

  if (aOrder < 10) then begin
    Product := 1.0;
    for i := 1 to aOrder do
      Product := Product * AsFloat;
    Result := -aMean * ln(Product) / aOrder;
  end
  else begin
    Result := AsGamma(aOrder, aMean);
  end;
end;
{--------}
function TStRandomBase.AsExponential(aMean : double) : double;
var
  R : double;
begin
  if not (aMean > 0.0) then
    raise EStPRNGError.Create(stscPRNGMeanS);

  repeat
    R := AsFloat;
  until (R <> 0.0);
  Result := -aMean * ln(R);
end;
{--------}
function TStRandomBase.AsF(aFreedom1 : integer;
                           aFreedom2 : integer) : double;
begin
  Result := (AsChiSquared(aFreedom1) * aFreedom1) /
            (AsChiSquared(aFreedom2) * aFreedom2);
end;
{--------}
function TStRandomBase.AsGamma(aShape : double; aScale : double) : double;
var
  R : double;
begin
  if not (aShape > 0.0) then
    raise EStPRNGError.Create(stscPRNGGammaShapeS);
  if not (aScale > 0.0) then
    raise EStPRNGError.Create(stscPRNGGammaScaleS);

  {there are three cases:
   ..0.0 < shape < 1.0, use Marsaglia's technique of
       Gamma(shape) = Gamma(shape+1) * uniform^(1/shape)}
  if (aShape < 1.0) then begin
    repeat
      R := AsFloat;
    until (R <> 0.0);
    Result := aScale * rbMarsagliaGamma(aShape + 1.0) *
              exp(ln(R) / aShape);
  end

  {..shape = 1.0: this is the same as exponential}
  else if (aShape = 1.0) then begin
    repeat
      R := AsFloat;
    until (R <> 0.0);
    Result := aScale * -ln(R);
  end

  {..shape > 1.0: use Marsaglia./Tsang algorithm}
  else begin
    Result := aScale * rbMarsagliaGamma(aShape);
  end;
end;
{--------}
function TStRandomBase.AsInt(aUpperLimit : integer) : integer;
begin
  if not (aUpperLimit > 0) then
    raise EStPRNGError.Create(stscPRNGLimitS);

  Result := Trunc(AsFloat * aUpperLimit);
end;
{--------}
function TStRandomBase.AsIntInRange(aLowerLimit : integer;
                                    aUpperLimit : integer) : integer;
begin
  if not (aLowerLimit < aUpperLimit) then
    raise EStPRNGError.Create(stscPRNGUpperLimitS);

  Result := Trunc(AsFloat * (aUpperLimit - aLowerLimit)) + ALowerLimit;
end;
{--------}
function TStRandomBase.AsLogNormal(aMean   : double;
                                   aStdDev : double) : double;
begin
  Result := exp(AsNormal(aMean, aStdDev));
end;
{--------}
function TStRandomBase.AsNormal(aMean   : double;
                                aStdDev : double) : double;
begin
  if not (aStdDev > 0.0) then
    raise EStPRNGError.Create(stscPRNGStdDevS);

  Result := (rbMontyPythonNormal * aStdDev) + aMean;

(*** alternative: The Box-Muller transformation
var
  R1, R2     : double;
  RadiusSqrd : double;
begin
  {get two random numbers that define a point in the unit circle}
  repeat
    R1 := (2.0 * aRandGen.AsFloat) - 1.0;
    R2 := (2.0 * aRandGen.AsFloat) - 1.0;
    RadiusSqrd := sqr(R1) + sqr(R2);
  until (RadiusSqrd < 1.0) and (RadiusSqrd > 0.0);

  {apply Box-Muller transformation}
  Result := (R1 * sqrt(-2.0 * ln(RadiusSqrd) / RadiusSqrd) * aStdDev)
            + aMean;
 ***)
end;
{--------}
function TStRandomBase.AsT(aFreedom : integer) : double;
begin
  if not (aFreedom > 0) then
    raise EStPRNGError.Create(stscPRNGDegFreedomS);

  Result := rbMontyPythonNormal / sqrt(AsChiSquared(aFreedom) / aFreedom);
end;
{--------}
function TStRandomBase.AsWeibull(aShape : double;
                   aScale : double) : double;
var
  R : double;
begin
  if not (aShape > 0) then
    raise EStPRNGError.Create(stscPRNGWeibullShapeS);
  if not (aScale > 0) then
    raise EStPRNGError.Create(stscPRNGWeibullScaleS);

  repeat
    R := AsFloat;
  until (R <> 0.0);
  Result := exp(ln(-ln(R)) / aShape) * aScale;
end;
{--------}

function TStRandomBase.rbMarsagliaGamma(aShape   : double) : double;
var
  d : double;
  c : double;
  x : double;
  v : double;
  u : double;
  Done : boolean;
begin
  {Notes: implements the Marsaglia/Tsang method of generating random
          numbers belonging to the gamma distribution:

          Marsaglia & Tsang, "A Simple Method for Generating Gamma
          Variables", ACM Transactions on Mathematical Software,
          Vol. 26, No. 3, September 2000, Pages 363-372

          It is pointless to try and work out what's going on in this
          routine without reading this paper :-)
          }

  d := aShape - (1.0 / 3.0);
  c := 1.0 / sqrt(9.0 * d);
  Done := false;
  {$IFDEF SuppressWarnings}
  v := 0.0;
  {$ENDIF}

  while not Done do begin
    repeat
      x := rbMontyPythonNormal;
      v := 1.0 + (c * x);
    until (v > 0.0);

    v := v * v * v;
    u := AsFloat;

    Done := u < (1.0 - 0.0331 * sqr(sqr(x)));

    if not Done then
      Done := ln(u) < (0.5 * sqr(x)) + d * (1.0 - v + ln(v))
  end;

  Result := d * v;
end;
{--------}
function TStRandomBase.rbMontyPythonNormal : double;
var
  x : double;
  y : double;
  v : double;
  NonZeroRandom : double;
begin
  {Notes: implements the Monty Python method of generating random
          numbers belonging to the Normal (Gaussian) distribution:

          Marsaglia & Tsang, "The Monty Python Method for Generating
          Random Variables", ACM Transactions on Mathematical
          Software, Vol. 24, No. 3, September 1998, Pages 341-350

          It is pointless to try and work out what's going on in this
          routine without reading this paper :-)

          Some constants:
          a = sqrt(ln(4))
          b = sqrt(2 * pi)
          s = a / (b - a)
          }

  {step 1: generate a random number x between +/- sqrt(2*Pi) and
           return it if its absolute value is less than sqrt(ln(4));
           note that this exit will happen about 47% of the time}
  x := ((AsFloat * 2.0) - 1.0) * Root2Pi;
  if (abs(x) < RootLn4) then begin
    Result := x;
    Exit;
  end;

  {step 2a: generate another random number y strictly between 0 and 1}
  repeat
    y := AsFloat;
  until (y <> 0.0);

  {step 2b: the first quadratic pretest avoids ln() calculation
            calculate v = 2.8658 - |x| * (2.0213 - 0.3605 * |x|)
            return x if y < v}
  v := 2.8658 - Abs(x) * (2.0213 - 0.3605 * Abs(x));
  if (y < v) then begin
    Result := x;
    Exit;
  end;

  {step 2c: the second quadratic pretest again avoids ln() calculation
            return s * (b - x) if y > v + 0.0506}
  if (y > v + 0.0506) then begin
    if (x > 0) then
      Result := MPN_s * (Root2Pi - x)
    else
      Result := -MPN_s * (Root2Pi + x);
    Exit;
  end;

  {step 2d: return x if y < f(x) or
              ln(y) < ln(2) - (0.5 * x * x) }
  if (ln(y) < (Ln2 - (0.5 * x * x))) then begin
    Result := x;
    Exit;
  end;

  {step 3: translate x to s * (b - x) and return it if y > g(x) or
             ln(1 + s - y) < ln(2 * s) - (0.5 * x * x) }
  if (x > 0) then
    x := MPN_s * (Root2Pi - x)
  else
    x := -MPN_s * (Root2Pi + x);
  if (ln(MPN_sPlus1 - y) < (Ln2MPN_s - (0.5 * x * x))) then begin
    Result := x;
    Exit;
  end;

  {step 4: the iterative process}
  repeat
    repeat
      NonZeroRandom := AsFloat;
    until (NonZeroRandom <> 0.0);
    x := -ln(NonZeroRandom) * InvRoot2Pi;
    repeat
      NonZeroRandom := AsFloat;
    until (NonZeroRandom <> 0.0);
    y := -ln(NonZeroRandom);
  until (y + y) > (x * x);
  if (NonZeroRandom < 0.5) then
    Result := -(Root2Pi + x)
  else
    Result := Root2Pi + x;
end;
{====================================================================}


{===TStRandomSystem==================================================}
constructor TStRandomSystem.Create(aSeed : integer);
begin
  inherited Create;
  Seed := aSeed;
end;
{--------}
function TStRandomSystem.AsFloat : double;
var
  SaveSeed : integer;
begin
  SaveSeed := RandSeed;
  RandSeed := FSeed;
  Result := System.Random;
  FSeed := RandSeed;
  RandSeed := SaveSeed;
end;
{--------}
procedure TStRandomSystem.rsSetSeed(aValue : integer);
begin
  if (aValue = 0) then
    FSeed := GetRandomSeed
  else
    FSeed := aValue;
end;
{====================================================================}


{===TStRandomCombined================================================}
const
  m1 = 2147483563;
  m2 = 2147483399;
{--------}
constructor TStRandomCombined.Create(aSeed1, aSeed2 : integer);
begin
  inherited Create;
  Seed1 := aSeed1;
  if (aSeed1 = 0) and (aSeed2 = 0) then
    Sleep(10); // a small delay to enable seed to change
  Seed2 := aSeed2;
end;
{--------}
function TStRandomCombined.AsFloat : double;
const
  a1 = 40014;
  q1 = 53668;  {equals m1 div a1}
  r1 = 12211;  {equals m1 mod a1}

  a2 = 40692;
  q2 = 52774;  {equals m2 div a2}
  r2 = 3791;   {equals m2 mod a2}

  OneOverM1 : double = 1.0 / m1;
var
  k : longint;
  Z : longint;
begin
  {advance first PRNG}
  k := FSeed1 div q1;
  FSeed1 := (a1 * (FSeed1 - (k * q1))) - (k * r1);
  if (FSeed1 < 0) then
    inc(FSeed1, m1);

  {advance second PRNG}
  k := FSeed2 div q2;
  FSeed2 := (a2 * (FSeed2 - (k * q2))) - (k * r2);
  if (FSeed2 < 0) then
    inc(FSeed2, m2);

  {combine the two seeds}
  Z := FSeed1 - FSeed2;
  if (Z <= 0) then
    Z := Z + m1 - 1;
  Result := Z * OneOverM1;
end;
{--------}
procedure TStRandomCombined.rcSetSeed1(aValue : integer);
begin
  if (aValue = 0) then
    FSeed1 := GetRandomSeed
  else
    FSeed1 := aValue;
end;
{--------}
procedure TStRandomCombined.rcSetSeed2(aValue : integer);
begin
  if (aValue = 0) then
    FSeed2 := GetRandomSeed
  else
    FSeed2 := aValue;
end;
{====================================================================}


{===TStRandomMother==================================================}
constructor TStRandomMother.Create(aSeed : integer);
begin
  inherited Create;
  Seed := aSeed;
end;
{--------}
function TStRandomMother.AsFloat : double;
const
  TwoM31 : double = 1.0 / $7FFFFFFF;
begin
  asm
    push esi
    push edi
    push ebx

    {get around a compiler bug where it doesn't notice that edx is
     being changed in the asm code !!! D5 bug}
    push edx

    {set ebx to point to self}
    mov ebx, eax

    {multiply X(n-4) by 21111111}
    mov eax, [ebx].TStRandomMother.FNMinus4
    mul [Mum1]
    mov edi, eax
    mov esi, edx

    {multiply X(n-3) by 1492 (save it in X(n-4) before though)}
    mov eax, [ebx].TStRandomMother.FNMinus3
    mov [ebx].TStRandomMother.FNMinus4, eax
    mul [Mum2]
    add edi, eax
    adc esi, edx

    {multiply X(n-2) by 1776 (save it in X(n-3) before though)}
    mov eax, [ebx].TStRandomMother.FNMinus2
    mov [ebx].TStRandomMother.FNMinus3, eax
    mul [Mum3]
    add edi, eax
    adc esi, edx

    {multiply X(n-1) by 5115 (save it in X(n-2) before though)}
    mov eax, [ebx].TStRandomMother.FNMinus1
    mov [ebx].TStRandomMother.FNMinus2, eax
    mul [Mum4]
    add edi, eax
    adc esi, edx

    {add in the remainder}
    add edi, [ebx].TStRandomMother.FC
    adc esi, 0;

    {save the lower 32 bits in X(n-1), the upper into the remainder}
    mov [ebx].TStRandomMother.FNMinus1, edi
    mov [ebx].TStRandomMother.FC, esi

    {get around a compiler bug where it doesn't notice that edx was
     changed in the asm code !!! D5 bug}
    pop edx

    pop ebx
    pop edi
    pop esi
  end;
  Result := (FNMinus1 shr 1) * TwoM31;
end;
{--------}
{$IFOPT Q+}
{note: TStRandomMother.rsSetSeed expressly overflows integers (it's
       equivalent to calculating mod 2^32), so we have to force
       overflow checks off}
{$DEFINE SaveQPlus}
{$Q-}
{$ENDIF}
procedure TStRandomMother.rsSetSeed(aValue : integer);
begin
  if (aValue = 0) then
    aValue := GetRandomSeed;
  FNminus4 := aValue;
  {note: the following code uses the generator
            Xn := (69069 * Xn-1) mod 2^32
         from D.E.Knuth, The Art of Computer Programming, Vol. 2
         (second edition), Addison-Wesley, 1981, pp.102}
  FNminus3 := 69069 * FNminus4;
  FNminus2 := 69069 * FNminus3;
  FNminus1 := 69069 * FNminus2;
  FC := 69069 * FNminus1;
end;
{$IFDEF SaveQPlus}
{$Q+}
{$ENDIF}
{====================================================================}


{====================================================================}
procedure CalcConstants;
begin
  {for the normal variates}
  Root2Pi := sqrt(2 * Pi);
  InvRoot2Pi := 1.0 / Root2Pi;
  RootLn4 := sqrt(ln(4.0));
  Ln2 := ln(2.0);
  MPN_s := RootLn4 / (Root2Pi - RootLn4);
  Ln2MPN_s := ln(2.0 * MPN_s);
  MPN_sPlus1 := MPN_s + 1.0;

  Mum1 := 2111111111;
  Mum2 := 1492;
  Mum3 := 1776;
  Mum4 := 5115;
end;
{====================================================================}


initialization
  CalcConstants;

end.
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.