Programmer 7500/MAX_PROGRAMMERS.iso/PASCAL/RANDOMTP.ZIP/RAND.PAS

Pascal/Delphi Source File | 1993-01-04 | 7.0 KB | 192 lines

{$A+,B-,D-,E-,F-,I-,L-,N-,O-,R-,S-,V+}
Unit Rand;
Interface
{ RAND.PAS - This unit implements the "minimal standard" random number
  generator as described in the October '88 "Communications of the acm"
  article by Park and Miller.  The algorithm is coded in Inline assembler
  using a variation of the technique described in the January, 1990
  CACM article by Carta.  The floating point result returned by function
  Random differs from the Park and Miller function only in the interval
  bounds: [0.0 <= Random < 1.0) for this unit as opposed to
  (0.0 < Random < 1.0) for Park/Miller.  The IRand Inline macro (see below)
  has been verified to yield the same seed value as the Park/Miller function
  after 3 million iterations, including several "overflow" conditions
  (see Carta).

  To eliminate the "parallel [hyper]planes" pattern in supercritical
  applications, a variation is provided that "Shuffles" the output of the
  minimal standard generator as described in Knuth, "The Art of Computer
  Programming", Vol 2 (Seminumerical Algorithms), page 32, Algorithm "B".

  Also provided is a routine to assist in selection of an unbiased sample
  from a population of known size.

  No attempt was made to do a statistical analysis of the output of these
  generators.  However, their performance has been analyzed extensively in
  the works cited above and their references.

  This work is contributed to the public domain by the author:

  Edwin T. Floyd
  #9 Adams Park Court
  Columbus, GA 31909
  [phone] (home)
  [phone] (work)

  First release, 1-15-90.
}
Const
  RandSeed : LongInt = 1;{ Random number seed: 0 < RandSeed < RandModulo }
  Shuffled : Boolean = False; { True if shuffle table is full }
  ShuffleBits = 5;       { 5 => 32-way shuffle, 6 => 64-way, etc }
  RandMult = 16807;      { =7**5; RandMult must be expressable in 16 bits.
                           48271 may give better "randomness" (see ACM ref.),
                           but shuffling is much more effective (see Knuth). }
  RandModulo = $7FFFFFFF;{ =2**31-1 = 2147483647 }

Function IRand : LongInt;
{ Return next "minimal standard", 31 bit pseudo-random integer.  This function
  actually computes (RandSeed * RandMult) Mod (2**31-1) where RandMult is
  a 16 bit quantity and RandSeed is 32 bits (See Carta, CACM 1/90).  Apart
  from the load/store of RandSeed, this may be the fastest possible
  implementation of this function with standard 8086 instructions on 16 bit
  registers. (The glove is down!) :-}
Inline(
  $A1/>RandSeed+2/       {         mov     ax,[>RandSeed+2]}
  $BF/>RandMult/         {         mov     di,>RandMult}
  $F7/$E7/               {         mul     di}
  $89/$C3/               {         mov     bx,ax}
  $89/$D1/               {         mov     cx,dx}
  $A1/>RandSeed/         {         mov     ax,[>RandSeed]}
  $F7/$E7/               {         mul     di}
  $01/$DA/               {         add     dx,bx}
  $83/$D1/$00/           {         adc     cx,0 ; cx:dx:ax = Seed * Mult }
  $D0/$E6/               {         shl     dh,1 ; split p & q at 31 bits }
  $D1/$D1/               {         rcl     cx,1}
  $D0/$EE/               {         shr     dh,1 ; cx = p, dx:ax = q }
  $01/$C8/               {         add     ax,cx}
  $83/$D2/$00/           {         adc     dx,0 ; dx:ax = p + q }
  $71/$09/               {         jno     done}
  $05/$01/$00/           {         add     ax,1 ; overflow, inc(p + q) }
  $83/$D2/$00/           {         adc     dx,0}
  $80/$E6/$7F/           {         and     dh,$7F ; limit to 31 bits }
                         {done:}
  $A3/>RandSeed/         {         mov     [>RandSeed],ax}
  $89/$16/>RandSeed+2);  {         mov     [>RandSeed+2],dx}

Procedure Randomize;
{ Initialize RandSeed using BIOS clock tick counter }

Function Random : Real;
{ Return pseudo-random, floating point number uniformly distributed on
  the interval [0.0 <= Random < 1.0) }

Function SIRand : LongInt;
{ Return shuffled, 31 bit, pseudo-random integer }

Function SRandom : Real;
{ Return shuffled, pseudo-random, floating point number uniformly distributed
  on the interval [0.0 <= SRandom < 1.0) }

Function SelectRecord(Sample, Population : LongInt;
  Var Selected, Index : LongInt) : Boolean;
{ Use this function to select an unbiased, random sample from a population of
  known size, without replacement. Set Index := 0 before the first reference to
  SelectRecord and test each record in the population for membership in
  sample via reference to SelectRecord with constant Sample <= Population
  parameters.  SelectRecord is guaranteed to return True for exactly "Sample"
  calls, randomly distributed over the population.  This routine will re-select
  the same sample given the same initial RandSeed and shuffle table state.
  For a discussion of the algorithm, see Knuth, "Art of Computer Programming",
  Vol 2 (Seminumerical Algorithms), page 137. Contrived example:

  Program DealOnePokerHand;
  Uses Rand;
  Var
    i, j, k : LongInt;
    CardFile : Text;
    Card : String;

  Begin
    Assign(CardFile, 'CARDNAME.TXT');
    Reset(CardFile);
    k := 0; (Start with zero index)
    Randomize;
    For i := 1 to 52 Do Begin
      ReadLn(CardFile, Card);
      If SelectRecord(5, 52, j, k) Then WriteLn(Card);
      (Select unbiased sample of 5 from a population of 52)
    End;
  End.
}

Implementation
Const
  RandomizeCallCount : Word = 1;
  ShuffleShift = 16 - ShuffleBits;
  ShufTableEnd = $FFFF Shr ShuffleShift;

Var
  ShufTable : Array[0..ShufTableEnd] Of LongInt;
  NextOut : Word;

Procedure Randomize;
Var
  BiosClock : LongInt Absolute $40:$6C;
  ClockLow : Byte Absolute $40:$6C;
  i : Byte;
Begin
  RandSeed := BiosClock And RandModulo;
  If RandSeed < 1 Then RandSeed := 1;
  If RandSeed >= RandModulo Then RandSeed := Pred(RandModulo);
  Shuffled := False;
  Inc(RandomizeCallCount);
  For i := 1 To (RandomizeCallCount + ClockLow) And $FF Do
    If Boolean(IRand) Then;
End;

Function Random : Real;
Begin
  Random := Pred(IRand) / Pred(RandModulo);
End;

Function SIRand : LongInt;
Var
  y : LongInt;
Begin
  If Not Shuffled Then Begin
    For NextOut := 0 To ShufTableEnd Do ShufTable[NextOut] := IRand;
    y := IRand;
    NextOut := Word(y) Shr ShuffleShift;
    Shuffled := True;
  End;
  y := ShufTable[NextOut];
  ShufTable[NextOut] := IRand;
  NextOut := Word(y) Shr ShuffleShift;
  SIRand := y;
End;

Function SRandom : Real;
Begin
  SRandom := Pred(SIRand) / Pred(RandModulo);
End;

Function SelectRecord(Sample, Population : LongInt;
  Var Selected, Index : LongInt) : Boolean;
Begin
  If Index = 0 Then Begin
    Selected := 0;
    If Sample > Population Then Begin
      WriteLn('Sample must be <= Population in RAND.SelectRecord');
      Halt(1);
    End;
  End;
  If (SRandom * (Population - Index) < (Sample - Selected))
  And (Selected < Sample) Then Begin
    SelectRecord := True;
    Inc(Selected);
  End Else SelectRecord := False;
  Inc(Index);
End;

End.