unit UNIform;

interface
{

    REAL FUNCTION UNI()
C***BEGIN PROLOGUE  UNI
C***DATE WRITTEN   810915 (YYMMDD)
C***REVISION DATE  871210 (YYMMDD)
C***CATEGORY NO.  L6A21
C***KEYWORDS  RANDOM NUMBERS, UNIFORM RANDOM NUMBERS
C***AUTHOR    KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS
C             MARSAGLIA, GEORGE, SUPERCOMPUTER RES. INST., FLORIDA ST. U.
C
C***PURPOSE  THIS ROUTINE GENERATES REAL (SINGLE PRECISION) UNIFORM
C             RANDOM NUMBERS ON [0,1)
C***DESCRIPTION
C        Computes real (single precision) uniform numbers on [0,1).
C           From the book, "Numerical Methods and Software" by
C                D. Kahaner, C. Moler, S. Nash
C                Prentice Hall, 1988
C
C       USAGE:
C              To initialize the generator
C                   USEED = USTART(ISEED)
C               where: ISEED is any NONZERO integer
C                  will return floating point value of ISEED.
C
C               Subsequently
C                       U = UNI()
C                  will return a real uniform on [0,1)
C
C                One initialization is necessary, but any number of evaluations
C                  of  UNI in any order, are allowed.
C
C           Note: Depending upon the value of K (see below), the output
C                       of UNI may differ from one machine to another.
C
C           Typical usage:
C
C               REAL U,UNI,USTART,USEED
C               INTEGER ISEED
CC                 Set seed
C               ISEED = 305
C               USEED = USTART(ISEED)
C               DO 1 I = 1,1000
C                   U = UNI()
C             1 CONTINUE
CC                 NOTE: If K=24 (the default, see below) the output value of
CC                           U will be 0.1570390462475...
C               WRITE(*,*) U
C               END
C
C          NOTE ON PORTABILITY: Users can choose to run UNI in its default
C               mode (requiring NO user action) which will generate the same
C               sequence of numbers on any computer supporting floating point
C               numbers with at least 24 bit mantissas, or in a mode that
C               will generate numbers with a longer period on computers with
C               larger mantissas.
C          TO EXERCISE THIS OPTION:  B E F O R E  invoking USTART insert
C               the instruction        UBITS = UNIB(K)      K >= 24
C               where K is the number of bits in the mantissa of your floating
C               point word (K=48 for Cray, Cyber 205). UNIB returns the
C               floating point value of K that it actually used.
C                    K input as .LE. 24, then UBITS=24.
C                    K input as .GT. 24, then UBITS=FLOAT(K)
C               If K>24 the sequence of numbers generated by UNI may differ
C               from one computer to another.
C
C
C
C***REFERENCES  MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM
C                 NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U.
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE UNI
      PARAMETER(
     *    CSAVE=362436./16777216.  ,
     *    CD=7654321./16777216.,
     *    CM=16777213./16777216.  )
C                            2**24=16777216
}


type
   TArr17Reals = array[1..17] of real;

  { T_UNIFORM_RND_GENERATOR }

  T_UNIFORM_RND_GENERATOR = class
  // inicializa el generador con la semilla
    constructor Create( ISEED: integer );

  // reinicializa con la semilla (retorna el mismo valor de ISEED)
    function Call_USTART(ISEED: integer): real;

    function Call_UNI: real;

    function Call_UNIB(KK: integer): real;

  private
    U_: TArr17Reals;
    I, J, K: integer;
    C: real;
  end;


implementation





const
  CSAVE = 362436.0 / 16777216.0;
  CD = 7654321.0 / 16777216.0;
  CM = 16777213.0 / 16777216.0;

  U_inicial: TArr17Reals = (0.8668672834288, 0.3697986366357, 0.8008968294805,
    0.4173889774680, 0.8254561579836, 0.9640965269077,
    0.4508667414265, 0.6451309529668, 0.1645456024730,
    0.2787901807898, 0.06761531340295, 0.9663226330820,
    0.01963343943798, 0.02947398211399, 0.1636231515294,
    0.3976343250467, 0.2631008574685);


{ T_UNIFORM_RND_GENERATOR }

constructor T_UNIFORM_RND_GENERATOR.Create;
var
  iPos: integer;

begin
  inherited Create;
  for ipos:= 1 to 17 do  U_[ipos]:= U_inicial[ipos];
  I:= 17;
  J:= 5;
  K:= 24;
  C:= CSAVE;
  Call_USTART( ISEED );
end;


function T_UNIFORM_RND_GENERATOR.Call_USTART(ISEED: integer): real;
{
      ENTRY USTART(ISEED)
C
C          Set up ...
C          Convert ISEED to four smallish positive integers.
C
}


var
  I1, J1, K1, L1: integer;
  II, JJ: integer;
  S, T: real;
  M1: integer;

begin
  I1 := ( ABS(ISEED) mod 177) + 1;
  J1 := ( ABS(ISEED) mod 167) + 1;
  K1 := ( ABS(ISEED) mod 157) + 1;
  L1 := ( ABS(ISEED) mod 147) + 1;


{C
C              Generate random bit pattern in array based on given seed.
C}

  for II := 1 to 17 do
  begin
    S := 0.0;
    T := 0.5;
{C             Do for each of the bits of mantissa of word
C             Loop  over K bits, where K is defaulted to 24 but can
C               be changed by user call to UNIB(K)}
    for JJ := 1 to K do
    begin
      M1 := ((( I1 * J1) mod  179) * K1) mod 179;
      I1 := J1;
      J1 := K1;
      K1 := M1;
      L1 := (53 * L1 + 1) mod 169;
      if ((L1 * M1) mod 64) >= 32 then
        S := S + T;
      T := 0.5 * T;
    end;
    U_[II] := S;
  end;
  Call_USTART := ISEED;
end;

function T_UNIFORM_RND_GENERATOR.Call_UNI: real;
var
  UNI: real;
  S, T, USTART, UNIB: real;
  II, JJ, KK, I1, J1, K1, L1, M1, ISEED: integer;
{
C
      SAVE U,I,J,K,C
C      Load data array in case user forgets to initialize.
C      This array is the result of calling UNI 100000 times
C         with ISEED=305 and K=64.
}
begin
{C
C   Basic generator is Fibonacci
C
}
  UNI := U_[I] - U_[J];
  if (UNI < 0.0) then
    UNI := UNI + 1.0;
  U_[I] := UNI;
  I := I - 1;
  if (I = 0) then
    I := 17;
  J := J - 1;
  if (J = 0) then
    J := 17;
{C
C   Second generator is congruential
C}
  C := C - CD;
  if (C < 0.0) then
    C := C + CM;
{C
C   Combination generator
C}
  UNI := UNI - C;
  if (UNI < 0.0) then
    UNI := UNI + 1.0;
  Call_UNI := UNI;
end;


function T_UNIFORM_RND_GENERATOR.Call_UNIB(KK: integer): real;
begin
  if KK <= 24 then
    K := 24
  else
    K := KK;
  Call_UNIB := K;
end;



end.
