strandom.pas

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

  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.