/*
**  random.src - pseudo-random number generators
**
** (C) Copyright 1988-1998 by Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
*/

#ifDLLCALL

/*
**  Format                   Purpose                                      Line
** ----------------------------------------------------------------------------
**    x = rndgam(r,c,alpha);     Gamma pseudo-random numbers                31
**    x = rndp(r,c,l);           Poisson pseudo-random numbers             104
**    x = rndnb(r,c,k,p);        Negative binomial pseudo-random numbers   165
**    x = rndbeta(r,c,a,b);      Beta pseudo-random numbers                260
**    x = rndvm(r,c,a,b);        von Mises pseudo-random numbers           355
**
**
**> rndgam
**
**  Purpose:  Computes Gamma pseudo-random numbers.
**
**  Format:  x = rndgam(r,c,alpha);
**
**  Input:  r      scalar, number of rows of resulting matrix.
**
**          c      scalar, number of columns of resulting matrix.
**
**          alpha  rxc matrix, or rx1 vector, or 1Xc vector, or scalar,
**                 shape argument for gamma distribution.
**
**  Output:  x     rxc matrix, gamma distributed random numbers.
**
**
**  Remarks:  The properties of the pseudo-random numbers in x are -
**
**            E(x) = alpha, Var(x) = alpha
**            x > 0,  alpha > 0.
**
**  To generate gamma(alpha,theta) pseudo-random numbers where theta
**  is a scale parameter, multiply the result of rndgam by theta.
**  Thus
**           z = theta * rndgam(1,1,alpha);
**
**  has the properties
**
**            E(z) = alpha * theta, Var(z) = alpha * theta ^ 2
**            z > 0, alpha > 0, theta > 0.
*/


proc rndgam(rows_x,cols_x,a);

   local x,rows_a,cols_a,ret,seed;

   x = zeros(rows_x,cols_x);
   rows_a = rows(a);
   cols_a = cols(a);

   if rows_a /= 1 and rows_a /= rows_x;
       if not trapchk(1);
           errorlog "rndgam: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if cols_a /= 1 and cols_a /= cols_x;
       if not trapchk(1);
           errorlog "rndgam: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if not (a > 0) and not trapchk(1);
       errorlog "rndgam:  parameter out of bounds";
       end;
   endif;

   seed = sysstate(28,0);
   dllcall rndgamdll(seed,x,rows_x,cols_x,a,rows_a,cols_a);
   call sysstate(28,seed);

   retp(x);
endp;


/*
**> rndp
**
**  Purpose:  Computes Poisson pseudo-random numbers
**
**  Format:  x = rndp(r,c,lambda);
**
**  Input:  r        scalar, number of rows of resulting matrix.
**
**          c        scalar, number of columns of resulting matrix.
**
**          lambda   rxc matrix, or rx1 vector, or 1Xc vector, or scalar,
**                   shape argument for Poisson distribution.
**
**  Output:  x       rxc matrix, Poisson distributed random numbers.
**
**  Remarks:  The properties of the pseudo-random numbers in x are -
**
**            E(x) = lambda, Var(x) = lambda
**            x = 0, 1, ....,  lambda > 0.
*/


proc rndp(rows_x,cols_x,l);
   local x,rows_l,cols_l,seed;

   x = zeros(rows_x,cols_x);
   rows_l = rows(l);
   cols_l = cols(l);

   if rows_l /= 1 and rows_l /= rows_x;
       if not trapchk(1);
           errorlog "rndp: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if cols_l /= 1 and cols_l /= cols_x;
       if not trapchk(1);
           errorlog "rndp: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if not (l > 0) and not trapchk(1);
       errorlog "rndgam:  parameter out of bounds";
       end;
   endif;

   seed = sysstate(28,0);
   dllcall rndpdll(seed,x,rows_x,cols_x,l,rows_l,cols_l);
   call sysstate(28,seed);

   retp(x);
endp;


/*
**> rndnb
**
**  Purpose:  Computes negative binomial pseudo-random numbers.
**
**  Format:  x = rndnb(r,c,k,p);
**
**  Input:   r    scalar, number of rows of resulting matrix.
**
**           c    scalar, number of columns of resulting matrix.
**
**           k    rxc matrix, or rx1 vector, or 1Xc vector, or scalar,
**                "event" argument for negative binomial distribution.
**
**           p    rxc matrix, or rx1 vector, or 1Xc vector, or scalar,
**                "probability" argument for negative binomial distribution.
**
**  Output:  x    rxc matrix, negative binomial distributed random numbers.
**
**
**  Remarks:  The properties of the pseudo-random numbers in x are -
**
**                      k * p                  k * p
**             E(x) = --------- , Var(x) = -----------
**                     (1 - p)              (1 - p)^2
**
**             x = 0, 1, ...,   k > 0,  0 < p < 1
*/



proc rndnb(rows_x,cols_x,k,p);
   local x,rows_k,cols_k,rows_p,cols_p,seed;

   x = zeros(rows_x,cols_x);
   rows_k = rows(k);
   cols_k = cols(k);
   rows_p = rows(p);
   cols_p = cols(p);

   if rows_k /= 1 and rows_k /= rows_x;
       if not trapchk(1);
           errorlog "rndnb: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if cols_k /= 1 and cols_k /= cols_x;
       if not trapchk(1);
           errorlog "rndnb: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if not (k > 0) and not trapchk(1);
       errorlog "rndnb:  parameter out of bounds";
       end;
   endif;


   if rows_p /= 1 and rows_p /= rows_x;
       if not trapchk(1);
           errorlog "rndnb: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if cols_p /= 1 and cols_p /= cols_x;
       if not trapchk(1);
           errorlog "rndnb: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if not (p > 0) and not (p < 1) and not trapchk(1);
       errorlog "rndnb:  parameter out of bounds";
       end;
   endif;

   seed = sysstate(28,0);
   dllcall rndnbdll(seed,x,rows_x,cols_x,k,rows_k,cols_k,p,rows_p,cols_p);
   call sysstate(28,seed);

   retp(x);
endp;


/*
**> rndbeta
**
**  Purpose:  Computes beta pseudo-random numbers.
**
**  Format:  x = rndbeta(r,c,a,b);
**
**  Input:   r    scalar, number of rows of resulting matrix.
**
**           c    scalar, number of columns of resulting matrix.
**
**           a    rxc matrix, or rx1 vector, or 1Xc vector, or scalar,
**                first shape argument for beta distribution.
**
**           b    rxc matrix, or rx1 vector, or 1Xc vector, or scalar,
**                second shape argument for beta distribution.
**
**  Output:  x    rxc matrix, beta distributed random numbers.
**
**
**  Remarks:  The properties of the pseudo-random numbers in x are -
**
**                       a                      (a * b)
**             E(x) = ------- ,   Var(x) = -------------------------
**                     a + b                (a + b + 1) * (a + b)^2
**
**             0 < x < 1,  a > 0,  b > 0
*/



proc rndbeta(rows_x,cols_x,a,b);
   local x,rows_a,cols_a,rows_b,cols_b,seed;

   x = zeros(rows_x,cols_x);
   rows_a = rows(a);
   cols_a = cols(a);
   rows_b = rows(b);
   cols_b = cols(b);

   if rows_a /= 1 and rows_a /= rows_x;
       if not trapchk(1);
           errorlog "rndbeta: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if cols_a /= 1 and cols_a /= cols_x;
       if not trapchk(1);
           errorlog "rndbeta: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if not (a > 0) and not trapchk(1);
       errorlog "rndbeta:  parameter out of bounds";
       end;
   endif;


   if rows_b /= 1 and rows_b /= rows_x;
       if not trapchk(1);
           errorlog "rndbeta: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if cols_b /= 1 and cols_b /= cols_x;
       if not trapchk(1);
           errorlog "rndbeta: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if not (b > 0) and not (b < 1) and not trapchk(1);
       errorlog "rndbeta:  parameter out of bounds";
       end;
   endif;

   seed = sysstate(28,0);
   dllcall rndbetadll(seed,x,rows_x,cols_x,a,rows_a,cols_a,b,rows_b,cols_b);
   call sysstate(28,seed);

   retp(x);
endp;


/*
**> rndvm
**
**  Purpose:  Computes von Mises pseudo-random numbers.
**
**  Format:  x = rndvm(r,c,m,k);
**
**  Input:   r    scalar, number of rows of resulting matrix.
**
**           c    scalar, number of columns of resulting matrix.
**
**           a    rxc matrix, or rx1 vector, or 1Xc vector, or scalar,
**                means for vm distribution.
**
**           b    rxc matrix, or rx1 vector, or 1Xc vector, or scalar,
**                shape argument for vm distribution.
**
**  Output:  x    rxc matrix, von Mises distributed random numbers.
**
*/



proc rndvm(rows_x,cols_x,a,b);
   local x,rows_a,cols_a,rows_b,cols_b,seed;

   x = zeros(rows_x,cols_x);
   rows_a = rows(a);
   cols_a = cols(a);
   rows_b = rows(b);
   cols_b = cols(b);

   if rows_a /= 1 and rows_a /= rows_x;
       if not trapchk(1);
           errorlog "rndvm: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if cols_a /= 1 and cols_a /= cols_x;
       if not trapchk(1);
           errorlog "rndvm: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if not (a > 0) and not trapchk(1);
       errorlog "rndvm:  parameter out of bounds";
       end;
   endif;


   if rows_b /= 1 and rows_b /= rows_x;
       if not trapchk(1);
           errorlog "rndvm: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if cols_b /= 1 and cols_b /= cols_x;
       if not trapchk(1);
           errorlog "rndvm: parameter matrix not conformable";
           end;
       else;
           retp(error(0));
       endif;
   endif;

   if not (b > 0) and not (b < 1) and not trapchk(1);
       errorlog "rndvm:  parameter out of bounds";
       end;
   endif;

   seed = sysstate(28,0);
   dllcall rndvmdll(seed,x,rows_x,cols_x,a,rows_a,cols_a,b,rows_b,cols_b);
   call sysstate(28,seed);

   retp(x);
endp;

#else

/*
**  rndp
**
**  Format      y = rndp(r,c,l);
**
**  Input       r   scalar, number of rows.
**              c   scalar, number of columns.
**
**              l   expected value(s) of Poisson variates.
**                  if 1x1 scalar, used for all variates.
**                  if rx1 vector, row of l used for each row of result.
**                  if cx1 vector, col of l used for each col of result.
**                  if rxc matrix, element of l used for each element of result.
**
**   Output     y   RxC matrix of random poisson numbers with parameters l.
*/


proc rndp(r,c,l);
    local i,j,z,l0;
    z = zeros(r,c);
    i = 1;
    do until i > r;
        j = 1;
        do until j > c;
            if cols(l) == c and rows(l) == r;
                l0 = l[i,j];
            elseif rows(l) == r;
                l0 = l[i,.];
            elseif cols(l) == c;
                l0 = l[.,j];
            else;
                l0 = l[1,1];
            endif;
            z[i,j] = _rndp(l0);
            j = j + 1;
       endo;
       i = i + 1;
    endo;
    retp(z);
endp;

proc _rndp(l);
    local x,r,r0;
    x = 0;
    r0 = rndu(1,1);
    r = exp(-l);
    do while r < R0;
       r0 = r0 - r;
       x = x + 1;
       r = r*l/x;
    endo;
    retp(x);
endp;

#endif