package RANDOM_GENERATORS is

--    UDRNRT returns a uniformly distributed random number between
--           0.0 and 1.0

  function UDRNRT return FLOAT ;

--    RANDOM returns a uniformly distributed random number between
--           0.0 and 1.0

  function RANDOM return FLOAT ;

--    GAUSS has parameters of the desired standard deviation and the
--          desired mean. The function returns Gaussian ( Normally )
--          distributed random numbers

  function GAUSS ( STDDEV , MEAN : FLOAT ) return FLOAT ;

--    POISSON has a parameter for the desired mean. The function
--    returns Poisson distributed random numbers

  function POISSON ( MEAN : FLOAT ) return FLOAT ;

  procedure P_TEST ( C : FLOAT ; X : in out FLOAT ; R : in out FLOAT);
    -- just for testing
  function P_DIST( LAMBDA : FLOAT; N : INTEGER; T : FLOAT) return FLOAT ;
    -- for comparison


--    The seed is visable so that repeatable sequences may be obtained
  SEED : INTEGER := 1 ;
  IY : INTEGER := 100001 ;

end RANDOM_GENERATORS ;

with ELEMENTARY_FUNCTIONS; use ELEMENTARY_FUNCTIONS;
package body RANDOM_GENERATORS is

  function UDRNRT return FLOAT is

--           MAGIC = M2 * HALF + M1
--           SEED  = S2 * HALF + S1
--           SEED := SEED * MAGIC  by          S1*M1
--                                          +  S2*M1
--                                          +  S1*M2
--                                        in order to prevent overflow
    MAGIC : constant := 65539 ;
    HALF : constant := 32768 ;
    FULL : constant := HALF * HALF ;
    MAGIC1 : constant := MAGIC mod HALF ;
    MAGIC2 : constant := MAGIC / HALF ;
    SEED1 : INTEGER ;
    SEED2 : INTEGER ;
    TEMP : INTEGER ;
    RANDOM : FLOAT ;

  begin

    SEED1 := SEED mod HALF ;
    SEED2 := SEED / HALF ;
    SEED := SEED1 * MAGIC1 ;
    SEED := SEED + ( ( SEED2 * MAGIC1 ) mod HALF ) * HALF ;
    SEED := SEED + ( ( SEED1 * MAGIC2 ) mod HALF ) * HALF ;
    SEED := SEED mod FULL ;
    RANDOM := FLOAT ( SEED ) / FLOAT ( FULL ) ;
    return RANDOM ;
  end UDRNRT ;


  function RANDOM return FLOAT is
    MAGIC : constant := 2796203 ;
  begin
    IY := IY * 125 ;
    IY := IY mod MAGIC ;
    return FLOAT ( IY ) / FLOAT ( MAGIC ) ;
  end RANDOM ;

  function GAUSS ( STDDEV , MEAN : FLOAT ) return FLOAT is

--     PURPOSE: GENERATE NORMALLY DISTRIBUTED RANDOM NUMBERS
--     INPUTS:
--         STDDEV - THE DESIRED STANDARD DEVIATION
--         MEAN   - THE DESIRED MEAN
--     OUTPUT:
--         GAUSS  - A NORMALLY DISTRIBUTED RANDOM NUMBER (GAUSSIAN DISTRIBUTION)
    A : FLOAT ;

  begin

    A := 0.0 ;
    for I in 1 .. 12 loop
      A := A + UDRNRT ;
    end loop ;
    return ( A - 6.0 ) * STDDEV + MEAN ;
  end GAUSS ;

  function POISSON ( MEAN : FLOAT ) return FLOAT is
    x1, x2, y1, y2, r : float;
  begin
    r := random;
    if r<=0.0 then return 0.0; end if;
    if r>0.9999 then return 12.0; end if;
    x1 := 1.0;
    y2 := r;
    for i in 1..20 loop -- just to prevent run away
      y1 := 1.0-(x1+1.0)*exp(-x1);
      x2 := x1+(y2-y1)/(x1*exp(-x1));
      exit when abs(x2-x1)/(abs(x2)+abs(x1)) < 2.0* float'epsilon;
      x1 := x2;
    end loop;
    return x2;
  end POISSON;

  procedure P_TEST ( C : FLOAT ; X : in out FLOAT ; R : in out FLOAT) is
    Y : FLOAT;
  begin
    R := -1.0;
    X := 1.0;
    while X > C loop
      Y := UDRNRT;
      X := X * Y;
      R := R + 1.0;
    end loop;
  end P_TEST;

  function P_DIST( LAMBDA : FLOAT; N : INTEGER; T : FLOAT) return FLOAT is
    N_FCT : FLOAT := 1.0; -- Q&D in floating point
  begin
    for I in 2..N loop
      N_FCT := N_FCT * FLOAT(I);
    end loop;
    return ((LAMBDA*T)**N) * EXP(-(LAMBDA*T)) / N_FCT;
  end P_DIST;


end RANDOM_GENERATORS ;


with RANDOM_GENERATORS; use RANDOM_GENERATORS;
with TEXT_IO; use TEXT_IO;
with INTEGER_TEXT_IO; use INTEGER_TEXT_IO;
with FLOAT_TEXT_IO; use FLOAT_TEXT_IO;
with REAL_ARRAYS; use REAL_ARRAYS;
with INTEGER_ARRAYS; use INTEGER_ARRAYS;
with INTEGER_ARRAYS_IO; use INTEGER_ARRAYS_IO;
with ELEMENTARY_FUNCTIONS; use ELEMENTARY_FUNCTIONS;
procedure TEST_POISSON is

  EXPD       : INTEGER_VECTOR(0..30) := (others => 0);
  EXPD_SCALE : FLOAT := 4.0;
  EXPD_LIM   : INTEGER := 1000*EXPD'LENGTH;

  Y : FLOAT;
  J : INTEGER;
  LAMBDA : FLOAT; -- DESIRED MEAN
  DATA_SAMPLE   : REAL_VECTOR(1..63_000);
  MEAN          : FLOAT; -- MEASURED MEAN
  STD_DEVIATION : FLOAT; -- MEASURED standard deviation
  R : FLOAT;
  C : FLOAT;

  procedure stddev(data    : real_vector;
                   mean    : out float;
                   std_dev : out float) is
    t_mean : float := 0.0;
    sum : float := 0.0;
    N : float := float(data'length);
  begin
    for i in data'range loop -- accumulate sum to get mean
      t_mean := t_mean + data(i);
    end loop;
    t_mean := t_mean / N; -- have mean
    mean := t_mean;
    for i in data'range loop -- accumulate deviations
      sum := sum + (data(i)-t_mean)*(data(i)-t_mean);
    end loop;
    std_dev := SQRT( sum / (N-1.0) );
  end stddev;

begin

  LAMBDA := 3.0;
  EXPD := (others => 0);
  PUT_LINE("POISSON DIST in" & INTEGER'IMAGE(EXPD'LENGTH) & " bins, " &
           INTEGER'IMAGE(EXPD_LIM) & " samples, mean at " &
           INTEGER'IMAGE(INTEGER(LAMBDA)) & "  scale=" &
           INTEGER'IMAGE(INTEGER(EXPD_SCALE)));
  PUT_LINE("              last bin has all larger values");
  for I in 1..EXPD_LIM loop
    Y := POISSON(3.0);
    DATA_SAMPLE(I) := Y;
    J := INTEGER(Y*EXPD_SCALE-0.5); -- simulate truncation
    if J < EXPD'FIRST then J := EXPD'FIRST; end if;
    if J > EXPD'LAST  then J := EXPD'LAST;  end if;
    EXPD(J) := EXPD(J) + 1;
  end loop;
  PUT(EXPD,5);
  NEW_LINE;
  EXPD := (others => 0);
  for I in 1..EXPD_LIM loop
    Y := DATA_SAMPLE(I);
    J := INTEGER(Y*EXPD_SCALE); -- normal averaging
    if J < EXPD'FIRST then J := EXPD'FIRST; end if;
    if J > EXPD'LAST  then J := EXPD'LAST;  end if;
    EXPD(J) := EXPD(J) + 1;
  end loop;
  PUT_LINE("Average around into ranges");
  PUT(EXPD,5);
  STDDEV(DATA_SAMPLE(1..EXPD_LIM), MEAN, STD_DEVIATION);
  PUT(MEAN); PUT(" = mean,   ");
  PUT(STD_DEVIATION); PUT_LINE(" = standard deviation");
  NEW_LINE;

  LAMBDA := 10.0;
  EXPD := (others => 0);
  EXPD_SCALE := 1.0;
  PUT_LINE("POISSON DIST in" & INTEGER'IMAGE(EXPD'LENGTH) & " bins, " &
           INTEGER'IMAGE(EXPD_LIM) & " samples, mean at " &
           INTEGER'IMAGE(INTEGER(LAMBDA)) & "  scale=" &
           INTEGER'IMAGE(INTEGER(EXPD_SCALE)));
  PUT_LINE("              last bin has all larger values");
  for I in 1..EXPD_LIM loop
    Y := POISSON(10.0);
    DATA_SAMPLE(I) := Y;
    J := INTEGER(Y*EXPD_SCALE-0.5); -- simulate truncation
    if J < EXPD'FIRST then J := EXPD'FIRST; end if;
    if J > EXPD'LAST  then J := EXPD'LAST;  end if;
    EXPD(J) := EXPD(J) + 1;
  end loop;
  PUT_LINE("Truncated into ranges");
  PUT(EXPD,5);
  NEW_LINE;
  EXPD := (others => 0);
  for I in 1..EXPD_LIM loop
    Y := DATA_SAMPLE(I);
    J := INTEGER(Y*EXPD_SCALE); -- normal averaging
    if J < EXPD'FIRST then J := EXPD'FIRST; end if;
    if J > EXPD'LAST  then J := EXPD'LAST;  end if;
    EXPD(J) := EXPD(J) + 1;
  end loop;
  PUT_LINE("Average around into ranges");
  for I in EXPD'RANGE loop
    INTEGER_TEXT_IO.PUT(I); PUT("   ");
    INTEGER_TEXT_IO.PUT(EXPD(I)); PUT("   ");
    PUT(1000.0*P_DIST(1.0,1,FLOAT(I)/EXPD_SCALE)); NEW_LINE;
  end loop;
  STDDEV(DATA_SAMPLE(1..EXPD_LIM), MEAN, STD_DEVIATION);
  PUT(MEAN); PUT(" = mean,   ");
  PUT(STD_DEVIATION); PUT_LINE(" = standard deviation");
  NEW_LINE;


  LAMBDA := 3.0;
  C := EXP(-LAMBDA);
  EXPD := (others => 0);
  EXPD_SCALE := 15.0*C;
  PUT_LINE("POISSON DIST in" & INTEGER'IMAGE(EXPD'LENGTH) & " bins, " &
           INTEGER'IMAGE(EXPD_LIM) & " samples, mean at " &
           INTEGER'IMAGE(INTEGER(LAMBDA)) & "  scale=" &
           INTEGER'IMAGE(INTEGER(EXPD_SCALE)));
  for I in 1..EXPD_LIM loop
    P_TEST(C, Y, R);
    DATA_SAMPLE(I) := Y;
    J := INTEGER(500.0*(Y-0.0248)+15.0); 
    if J < EXPD'FIRST then J := EXPD'FIRST; end if;
    if J > EXPD'LAST  then J := EXPD'LAST;  end if;
    EXPD(J) := EXPD(J) + 1;
  end loop;
  PUT(EXPD,5);
  STDDEV(DATA_SAMPLE(1..EXPD_LIM), MEAN, STD_DEVIATION);
  PUT(MEAN); PUT(" = mean,   ");
  PUT(STD_DEVIATION); PUT_LINE(" = standard deviation");
  NEW_LINE;

end TEST_POISSON;