systools: Add units StRandom and StStat (and demo for StRandom).

git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@6141 8e941d3f-bd1b-0410-a28a-d453659cc2b4
This commit is contained in:
wp_xxyyzz
2018-01-17 08:04:35 +00:00
parent 93e37e8e76
commit 560fd631fa
9 changed files with 3885 additions and 4 deletions

View File

@@ -0,0 +1,84 @@
<?xml version="1.0" encoding="UTF-8"?>
<CONFIG>
<ProjectOptions>
<Version Value="11"/>
<PathDelim Value="\"/>
<General>
<Flags>
<UseDefaultCompilerOptions Value="True"/>
</Flags>
<SessionStorage Value="InProjectDir"/>
<MainUnit Value="0"/>
<Title Value="exrandom"/>
<Scaled Value="True"/>
<ResourceType Value="res"/>
<UseXPManifest Value="True"/>
<XPManifest>
<DpiAware Value="True"/>
</XPManifest>
<Icon Value="0"/>
</General>
<BuildModes Count="1">
<Item1 Name="Default" Default="True"/>
</BuildModes>
<PublishOptions>
<Version Value="2"/>
</PublishOptions>
<RunParams>
<FormatVersion Value="2"/>
<Modes Count="0"/>
</RunParams>
<RequiredPackages Count="2">
<Item1>
<PackageName Value="laz_systools"/>
</Item1>
<Item2>
<PackageName Value="LCL"/>
</Item2>
</RequiredPackages>
<Units Count="2">
<Unit0>
<Filename Value="exrandom.lpr"/>
<IsPartOfProject Value="True"/>
</Unit0>
<Unit1>
<Filename Value="exrndu.pas"/>
<IsPartOfProject Value="True"/>
<ComponentName Value="Form1"/>
<HasResources Value="True"/>
<ResourceBaseClass Value="Form"/>
<UnitName Value="ExRndU"/>
</Unit1>
</Units>
</ProjectOptions>
<CompilerOptions>
<Version Value="11"/>
<PathDelim Value="\"/>
<Target>
<Filename Value="exrandom"/>
</Target>
<SearchPaths>
<UnitOutputDirectory Value="lib\$(TargetCPU)-$(TargetOS)"/>
</SearchPaths>
<Linking>
<Options>
<Win32>
<GraphicApplication Value="True"/>
</Win32>
</Options>
</Linking>
</CompilerOptions>
<Debugging>
<Exceptions Count="3">
<Item1>
<Name Value="EAbort"/>
</Item1>
<Item2>
<Name Value="ECodetoolError"/>
</Item2>
<Item3>
<Name Value="EFOpenError"/>
</Item3>
</Exceptions>
</Debugging>
</CONFIG>

View File

@@ -0,0 +1,43 @@
(* ***** 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 ***** *)
{$IFDEF FPC}
{$mode DELPHI}
{$ENDIF}
program ExRandom;
uses
Interfaces,
Forms,
ExRndU in 'ExRndU.pas' {Form1};
{$R *.res}
begin
Application.Initialize;
Application.CreateForm(TForm1, Form1);
Application.Run;
end.

View File

@@ -0,0 +1,146 @@
object Form1: TForm1
Left = 192
Height = 498
Top = 114
Width = 462
Caption = 'Random Distributions '
ClientHeight = 498
ClientWidth = 462
Color = clBtnFace
Font.Color = clWindowText
OnCreate = FormCreate
OnDestroy = FormDestroy
LCLVersion = '1.9.0.0'
object imgGraph: TImage
Left = 32
Height = 250
Top = 200
Width = 400
end
object lblPrompt: TLabel
Left = 32
Height = 15
Top = 24
Width = 159
Caption = 'Select the distribution to view:'
ParentColor = False
end
object lblGraphTitle: TLabel
Left = 32
Height = 1
Top = 184
Width = 1
ParentColor = False
end
object lblParms: TLabel
Left = 32
Height = 15
Top = 56
Width = 62
Caption = 'Parameters:'
ParentColor = False
end
object lblParm1: TLabel
Left = 48
Height = 1
Top = 80
Width = 1
ParentColor = False
end
object lblParm2: TLabel
Left = 48
Height = 1
Top = 104
Width = 1
ParentColor = False
end
object lblLeft: TLabel
Left = 32
Height = 1
Top = 456
Width = 1
ParentColor = False
end
object lblRight: TLabel
Left = 429
Height = 1
Top = 456
Width = 1
Alignment = taRightJustify
ParentColor = False
end
object lblMaxY: TLabel
Left = 8
Height = 1
Top = 184
Width = 1
ParentColor = False
end
object cboDist: TComboBox
Left = 200
Height = 23
Top = 21
Width = 145
ItemHeight = 15
OnChange = cboDistChange
Sorted = True
Style = csDropDownList
TabOrder = 0
end
object btnGenerate: TButton
Left = 32
Height = 25
Top = 144
Width = 145
Caption = 'Generate graph'
Default = True
OnClick = btnGenerateClick
TabOrder = 3
end
object prgGenProgress: TProgressBar
Left = 32
Height = 16
Top = 472
Width = 400
TabOrder = 6
end
object edtParm1: TEdit
Left = 200
Height = 23
Top = 74
Width = 121
MaxLength = 10
TabOrder = 1
end
object edtParm2: TEdit
Left = 200
Height = 23
Top = 98
Width = 121
MaxLength = 10
TabOrder = 2
end
object updRight: TUpDown
Left = 440
Height = 24
Top = 448
Width = 16
Min = 1
OnClick = updRightClick
Position = 1
TabOrder = 5
Wrap = False
end
object updLeft: TUpDown
Left = 8
Height = 24
Top = 448
Width = 16
Max = 0
Min = -100
OnClick = updLeftClick
Position = 0
TabOrder = 4
Wrap = False
end
end

View File

@@ -0,0 +1,523 @@
(* ***** 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 ***** *)
{$IFDEF FPC}
{$mode DELPHI}
{$ENDIF}
unit ExRndU;
interface
uses
{$IFNDEF FPC}
Windows, Messages,
{$ENDIF}
SysUtils, Classes, Graphics, Controls, Forms,
Dialogs, ComCtrls, StdCtrls, ExtCtrls,
StRandom;
type
TGetRandom = function : double of object;
type
TForm1 = class(TForm)
imgGraph: TImage;
cboDist: TComboBox;
lblPrompt: TLabel;
btnGenerate: TButton;
prgGenProgress: TProgressBar;
lblGraphTitle: TLabel;
lblParms: TLabel;
lblParm1: TLabel;
lblParm2: TLabel;
edtParm1: TEdit;
edtParm2: TEdit;
lblLeft: TLabel;
lblRight: TLabel;
updRight: TUpDown;
updLeft: TUpDown;
lblMaxY: TLabel;
procedure btnGenerateClick(Sender: TObject);
procedure cboDistChange(Sender: TObject);
procedure FormCreate(Sender: TObject);
procedure updRightClick(Sender: TObject; Button: TUDBtnType);
procedure updLeftClick(Sender: TObject; Button: TUDBtnType);
procedure FormDestroy(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
GraphLeft : double;
GraphRight : double;
Value1 : double;
Value2 : double;
PRNG : TStRandomBase;
GetRandom : TGetRandom;
procedure GenerateGraph(aDistInx : integer);
procedure PrepForBeta;
procedure PrepForCauchy;
procedure PrepForChiSquared;
procedure PrepForErlang;
procedure PrepForExponential;
procedure PrepForF;
procedure PrepForGamma;
procedure PrepForLogNormal;
procedure PrepForNormal;
procedure PrepForT;
procedure PrepForUniform;
procedure PrepForWeibull;
function GetBeta : double;
function GetCauchy : double;
function GetChiSquared : double;
function GetErlang : double;
function GetExponential : double;
function GetF : double;
function GetGamma : double;
function GetLogNormal : double;
function GetNormal : double;
function GetT : double;
function GetUniform : double;
function GetWeibull : double;
end;
var
Form1: TForm1;
implementation
{$R *.lfm}
const
DistNames : array [0..11] of string = (
'Beta', 'Cauchy', 'ChiSquared', 'Erlang', 'Exponential',
'F', 'Gamma', 'LogNormal', 'Normal', 'Student''s t',
'Uniform', 'Weibull');
const
RandomCount = 1000000;
procedure TForm1.GenerateGraph(aDistInx : integer);
var
Buckets : array[0..400] of integer;
i : integer;
R : double;
Inx : integer;
MaxHt : integer;
MaxLineFactor : double;
GraphWidth : double;
OldPercent : integer;
NewPercent : integer;
begin
{zero out the buckets}
FillChar(Buckets, sizeof(Buckets), 0);
{calculate random numbers according to distribution, convert to a
bucket index, and increment that bucket count}
OldPercent := -1;
GraphWidth := imgGraph.Width;
for i := 1 to RandomCount do begin
NewPercent := (i * 100) div RandomCount;
if (NewPercent <> OldPercent) then begin
prgGenProgress.Position := NewPercent;
OldPercent := NewPercent;
end;
R := GetRandom;
if (GraphLeft <= R) and (R <= GraphRight) then begin
Inx := trunc((R - GraphLeft) * GraphWidth / (GraphRight - GraphLeft));
if (0 <= Inx) and (Inx <= 400) then
inc(Buckets[Inx]);
end;
end;
{calculate the largest bucket}
MaxHt := 1;
for i := 0 to 400 do
if (MaxHt < Buckets[i]) then
MaxHt := Buckets[i];
{draw the graph}
imgGraph.Canvas.Lock;
try
imgGraph.Canvas.FillRect(Rect(0, 0, imgGraph.Width, imgGraph.Height));
MaxLineFactor := imgGraph.Height / MaxHt;
imgGraph.Canvas.Pen.Color := clRed;
for i := 0 to 400 do begin
imgGraph.Canvas.PenPos := Point(i, imgGraph.Height);
imgGraph.Canvas.LineTo(i, imgGraph.Height - trunc(Buckets[i] * MaxLineFactor));
end;
finally
imgGraph.Canvas.Unlock;
end;
lblMaxY.Caption := Format('Max: %8.6f', [MaxHt / RandomCount]);
end;
procedure TForm1.btnGenerateClick(Sender: TObject);
begin
if (edtParm1.Text = '') then
Value1 := 0.0
else
Value1 := StrToFloat(edtParm1.Text);
if (edtParm2.Text = '') then
Value2 := 0.0
else
Value2 := StrToFloat(edtParm2.Text);
GenerateGraph(cboDist.ItemIndex);
end;
procedure TForm1.cboDistChange(Sender: TObject);
begin
case cboDist.ItemIndex of
0 : PrepForBeta;
1 : PrepForCauchy;
2 : PrepForChiSquared;
3 : PrepForErlang;
4 : PrepForExponential;
5 : PrepForF;
6 : PrepForGamma;
7 : PrepForLogNormal;
8 : PrepForNormal;
9 : PrepForT;
10: PrepForUniform;
11: PrepForWeibull
end;
updRightClick(Self, btNext);
updLeftClick(Self, btNext);
edtParm1.Text := FloatToStr(Value1);
edtParm2.Text := FloatToStr(Value2);
end;
procedure TForm1.PrepForBeta;
begin
lblParm1.Caption := 'Shape 1:';
lblParm1.Visible := true;
lblParm2.Caption := 'Shape 2:';
lblParm2.Visible := true;
edtParm1.Visible := true;
edtParm1.Enabled := true;
edtParm2.Visible := true;
edtParm2.Enabled := true;
updLeft.Position := 0;
updRight.Position := 1;
Value1 := 2.0;
Value2 := 4.0;
GetRandom := GetBeta;
end;
procedure TForm1.PrepForCauchy;
begin
lblParm1.Caption := '(none)';
lblParm1.Visible := true;
lblParm2.Visible := false;
edtParm1.Visible := false;
edtParm1.Enabled := false;
edtParm2.Visible := false;
edtParm2.Enabled := false;
updLeft.Position := -5;
updRight.Position := 5;
Value1 := 0.0;
Value2 := 0.0;
GetRandom := GetCauchy;
end;
procedure TForm1.PrepForChiSquared;
begin
lblParm1.Caption := 'Degrees of freedom:';
lblParm1.Visible := true;
lblParm2.Visible := false;
edtParm1.Visible := true;
edtParm1.Enabled := true;
edtParm2.Visible := false;
edtParm2.Enabled := false;
updLeft.Position := 0;
updRight.Position := 20;
Value1 := 5.0;
Value2 := 0.0;
GetRandom := GetChiSquared;
end;
procedure TForm1.PrepForErlang;
begin
lblParm1.Caption := 'Mean:';
lblParm1.Visible := true;
lblParm2.Caption := 'Order:';
lblParm2.Visible := true;
edtParm1.Visible := true;
edtParm1.Enabled := true;
edtParm2.Visible := true;
edtParm2.Enabled := true;
updLeft.Position := 0;
updRight.Position := 5;
Value1 := 1.0;
Value2 := 4.0;
GetRandom := GetErlang;
end;
procedure TForm1.PrepForExponential;
begin
lblParm1.Caption := 'Mean:';
lblParm1.Visible := true;
lblParm2.Visible := false;
edtParm1.Visible := true;
edtParm1.Enabled := true;
edtParm2.Visible := false;
edtParm2.Enabled := false;
updLeft.Position := 0;
updRight.Position := 10;
Value1 := 1.0;
Value2 := 0.0;
GetRandom := GetExponential;
end;
procedure TForm1.PrepForF;
begin
lblParm1.Caption := 'Degrees of freedom 1:';
lblParm1.Visible := true;
lblParm2.Caption := 'Degrees of freedom 2:';
lblParm2.Visible := true;
edtParm1.Visible := true;
edtParm1.Enabled := true;
edtParm2.Visible := true;
edtParm2.Enabled := true;
updLeft.Position := 0;
updRight.Position := 20;
Value1 := 10.0;
Value2 := 5.0;
GetRandom := GetF;
end;
procedure TForm1.PrepForGamma;
begin
lblParm1.Caption := 'Shape:';
lblParm1.Visible := true;
lblParm2.Caption := 'Scale:';
lblParm2.Visible := true;
edtParm1.Visible := true;
edtParm1.Enabled := true;
edtParm2.Visible := true;
edtParm2.Enabled := true;
updLeft.Position := 0;
updRight.Position := 10;
Value1 := 2.0;
Value2 := 1.0;
GetRandom := GetGamma;
end;
procedure TForm1.PrepForLogNormal;
begin
lblParm1.Caption := 'Mean:';
lblParm1.Visible := true;
lblParm2.Caption := 'Standard deviation:';
lblParm2.Visible := true;
edtParm1.Visible := true;
edtParm1.Enabled := true;
edtParm2.Visible := true;
edtParm2.Enabled := true;
updLeft.Position := 0;
updRight.Position := 10;
Value1 := 0.0;
Value2 := 1.0;
GetRandom := GetLogNormal;
end;
procedure TForm1.PrepForNormal;
begin
lblParm1.Caption := 'Mean:';
lblParm1.Visible := true;
lblParm2.Caption := 'Standard deviation:';
lblParm2.Visible := true;
edtParm1.Visible := true;
edtParm1.Enabled := true;
edtParm2.Visible := true;
edtParm2.Enabled := true;
updLeft.Position := -5;
updRight.Position := 5;
Value1 := 0.0;
Value2 := 1.0;
GetRandom := GetNormal;
end;
procedure TForm1.PrepForT;
begin
lblParm1.Caption := 'Degrees of freedom:';
lblParm1.Visible := true;
lblParm2.Visible := false;
edtParm1.Visible := true;
edtParm1.Enabled := true;
edtParm2.Visible := false;
edtParm2.Enabled := false;
updLeft.Position := -10;
updRight.Position := 10;
Value1 := 10.0;
Value2 := 0.0;
GetRandom := GetT;
end;
procedure TForm1.PrepForUniform;
begin
lblParm1.Caption := '(none)';
lblParm1.Visible := true;
lblParm2.Visible := false;
edtParm1.Visible := false;
edtParm1.Enabled := false;
edtParm2.Visible := false;
edtParm2.Enabled := false;
updLeft.Position := 0;
updRight.Position := 1;
Value1 := 0.0;
Value2 := 0.0;
GetRandom := GetUniform;
end;
procedure TForm1.PrepForWeibull;
begin
lblParm1.Caption := 'Shape:';
lblParm1.Visible := true;
lblParm2.Caption := 'Scale:';
lblParm2.Visible := true;
edtParm1.Visible := true;
edtParm1.Enabled := true;
edtParm2.Visible := true;
edtParm2.Enabled := true;
updLeft.Position := 0;
updRight.Position := 10;
Value1 := 2.0;
Value2 := 3.0;
GetRandom := GetWeibull;
end;
function TForm1.GetBeta : double;
begin
Result := PRNG.AsBeta(Value1, Value2)
end;
function TForm1.GetCauchy : double;
begin
Result := PRNG.AsCauchy
end;
function TForm1.GetChiSquared : double;
begin
if (Value1 > 65535.0) then
raise Exception.Create(
'TForm1.GetChiSquared: the degrees of freedom value 1 is too large for this example program');
Result := PRNG.AsChiSquared(trunc(Value1))
end;
function TForm1.GetErlang : double;
begin
Result := PRNG.AsErlang(Value1, trunc(Value2))
end;
function TForm1.GetExponential : double;
begin
Result := PRNG.AsExponential(Value1)
end;
function TForm1.GetF : double;
begin
if (Value1 > 65535.0) then
raise Exception.Create(
'TForm1.GetF: the degrees of freedom value 1 is too large for this example program');
if (Value2 > 65535.0) then
raise Exception.Create(
'TForm1.GetF: the degrees of freedom value 2 is too large for this example program');
Result := PRNG.AsF(trunc(Value1), trunc(Value2))
end;
function TForm1.GetGamma : double;
begin
Result := PRNG.AsGamma(Value1, Value2)
end;
function TForm1.GetLogNormal : double;
begin
Result := PRNG.AsLogNormal(Value1, Value2)
end;
function TForm1.GetNormal : double;
begin
Result := PRNG.AsNormal(Value1, Value2)
end;
function TForm1.GetT : double;
begin
if (Value1 > 65535.0) then
raise Exception.Create(
'TForm1.GetT: the degrees of freedom value is too large for this example program');
Result := PRNG.AsT(trunc(Value1))
end;
function TForm1.GetUniform : double;
begin
Result := PRNG.AsFloat
end;
function TForm1.GetWeibull : double;
begin
Result := PRNG.AsWeibull(Value1, Value2)
end;
procedure TForm1.FormCreate(Sender: TObject);
var
i : integer;
UniformInx : integer;
begin
cboDist.Items.Clear;
UniformInx := -1;
for i := 0 to high(DistNames) do begin
cboDist.Items.Add(DistNames[i]);
if (Copy(DistNames[i], 1, 7) = 'Uniform') then
UniformInx := i;
end;
cboDist.ItemIndex := UniformInx;
cboDistChange(Self);
PRNG := TStRandomSystem.Create(0);
end;
procedure TForm1.updRightClick(Sender: TObject; Button: TUDBtnType);
begin
lblRight.Caption := IntToStr(updRight.Position);
GraphRight := updRight.Position;
end;
procedure TForm1.updLeftClick(Sender: TObject; Button: TUDBtnType);
begin
lblLeft.Caption := IntToStr(updLeft.Position);
GraphLeft := updLeft.Position;
end;
procedure TForm1.FormDestroy(Sender: TObject);
begin
PRNG.Free;
end;
end.

View File

@@ -16,7 +16,7 @@
<Description Value="Lazarus port of TurboPower SysTools (Delphi version at https://sourceforge.net/projects/tpsystools/) - runtime package"/> <Description Value="Lazarus port of TurboPower SysTools (Delphi version at https://sourceforge.net/projects/tpsystools/) - runtime package"/>
<License Value="MPL 1.1"/> <License Value="MPL 1.1"/>
<Version Major="4" Release="4"/> <Version Major="4" Release="4"/>
<Files Count="21"> <Files Count="23">
<Item1> <Item1>
<Filename Value="source\run\stbarc.pas"/> <Filename Value="source\run\stbarc.pas"/>
<UnitName Value="StBarC"/> <UnitName Value="StBarC"/>
@@ -101,6 +101,14 @@
<Filename Value="source\run\stmoney.pas"/> <Filename Value="source\run\stmoney.pas"/>
<UnitName Value="StMoney"/> <UnitName Value="StMoney"/>
</Item21> </Item21>
<Item22>
<Filename Value="source\run\strandom.pas"/>
<UnitName Value="StRandom"/>
</Item22>
<Item23>
<Filename Value="source\run\ststat.pas"/>
<UnitName Value="StStat"/>
</Item23>
</Files> </Files>
<RequiredPkgs Count="2"> <RequiredPkgs Count="2">
<Item1> <Item1>

View File

@@ -10,7 +10,7 @@ interface
uses uses
StBarC, StBase, StConst, StBarPN, StStrL, St2DBarC, StDate, StUtils, StCRC, StBarC, StBase, StConst, StBarPN, StStrL, St2DBarC, StDate, StUtils, StCRC,
StHASH, StToHTML, StStrms, StDict, StIniStm, StDecMth, StExpr, StMath, StHASH, StToHTML, StStrms, StDict, StIniStm, StDecMth, StExpr, StMath,
StFIN, StDateSt, StMoney; StFIN, StDateSt, StMoney, StRandom, StStat;
implementation implementation

View File

@@ -126,10 +126,12 @@ uses
StRegIni, StRegIni,
StSaturn, StSaturn,
StSort, StSort,
*)
StStat, StStat,
StStrL, StStrL,
StStrms, StStrms,
StStrS, StStrS,
(*
StStrW, StStrW,
StStrZ, StStrZ,
StText, StText,
@@ -141,15 +143,17 @@ uses
StVArr, StVArr,
StVenus, StVenus,
{ new units in ver 4: } { new units in ver 4: }
*)
StIniStm, StIniStm,
(*
StMerge, StMerge,
StSystem, StSystem,
StTxtDat, StTxtDat,
StDecMth, StDecMth,
*) *)
StMoney StMoney,
StRandom
(* (*
StRandom,
StNTLog, StNTLog,
{ !!! StExpEng unit designed to handle problem with initialization } { !!! StExpEng unit designed to handle problem with initialization }
{ section in C++Builder; should NOT be included in Registration unit } { section in C++Builder; should NOT be included in Registration unit }

View File

@@ -0,0 +1,735 @@
// 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 *}
{*********************************************************}
{$IFDEF FPC}
{$mode DELPHI}
{$ENDIF}
//{$I StDefine.inc}
unit StRandom;
interface
uses
{$IFNDEF FPC}
Windows,
{$ENDIF}
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.

File diff suppressed because it is too large Load Diff