strandom.pas

条码控件： 一维条码控件 二维条码控件 PDF417Barcode MaxiCodeBarcode

(* ***** 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.03                           *}
{*********************************************************}
{* 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);