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