/**************************************************************************
 *                                     file pruec                         *
 **************************************************************************/
#include<math.h>

extern double drand(unsigned long *seed);

/******************************************************************
 *                                                                *
 *       Poisson Distribution - Ratio of Uniforms/Inversion       *
 *                                                                *
 ******************************************************************
 *                                                                *
 * For parameter  my < 5.5 the Inversion method is applied:       *
 * The random numbers are generated via sequential search,        *
 * starting at the lowest index k=0. The cumulative probabilities *
 * are avoided by using the technique of chop-down.               *
 * For my >= 5.5  Ratio of Uniforms is employed:                  *
 * A table mountain hat function h(x) with optimal scale          *
 * parameter s for fixed location parameter  a = mu+1/2  is used. *
 *                                                                *
 ******************************************************************
 *                                                                *
 * FUNCTION:    - pruec samples a random number from the Poisson  *
 *                distribution with parameter my > 0.             *
 * REFERENCE:   - E. Stadlober (1989): Sampling from Poisson,     *
 *                binomial and hypergeometric distributions:      *
 *                ratio of uniforms as a simple and fast          *
 *                alternative, Bericht 303, Math. Stat. Sektion,  *
 *                Forschungsgesellschaft Joanneum, Graz.          *
 * SUBPROGRAMS: - drand(seed) ... (0,1)-Uniform generator with    *
 *                unsigned long integer *seed.                    *
 *                                                                *
 * Implemented by R.Kremer 1990, revised by P.Busswald, July 1992 *
 ******************************************************************/

long pruec(unsigned long *seed, double my)

{
  #define       C0      9.18938533204672742e-01       /* constants for */
  #define       C1      8.33333333333333333e-02       /* Stirling's    */
  #define       C3     -2.77777777777777778e-03       /* approximation */
  #define       C5      7.93650793650793651e-04       /* of ln(k!)     */
  #define       C7     -5.95238095238095238e-04

  static long m,b;

  static double logfak[30L] = {
          0.00000000000000000,   0.00000000000000000,   0.69314718055994531,
          1.79175946922805500,   3.17805383034794562,   4.78749174278204599,
          6.57925121201010100,   8.52516136106541430,  10.60460290274525023,
         12.80182748008146961,  15.10441257307551530,  17.50230784587388584,
         19.98721449566188615,  22.55216385312342289,  25.19122118273868150,
         27.89927138384089157,  30.67186010608067280,  33.50507345013688888,
         36.39544520803305358,  39.33988418719949404,  42.33561646075348503,
         45.38013889847690803,  48.47118135183522388,  51.60667556776437357,
         54.78472939811231919,  58.00360522298051994,  61.26170176100200198,
         64.55753862700633106,  67.88974313718153498,  71.25703896716800901
  };                  /* ln(k!) for k=0,1,2,...,29 */
 static double my_prev = -1.0,p0,a,h,g,q;
 long k;
 double u,x,lf,ry,ry2,yk,ym,fk,fm,r;

 if (my_prev != my)                                  /* Set-up */
        {
         my_prev = my;
         if (my<5.5)
                 {
        p0=exp(-my);                            /* Set-up for Inversion */
        b=(long int)(my+10.0*(sqrt(my)+1.0));   /* safety-bound */
                 }
         else
                 {
        a = my + 0.5;                      /* Set-up for Ratio of Uniforms */
        r = sqrt(a + a);
        b = floor(a + 7.0 * (r + 1.5));    /* safety-bound */
        m = my;
        g  = log(my);
        ym = m;
        yk = floor(a - r);
        x = (a - yk - 1.0) / (a - yk);
        if (my * x * x > yk + 1) yk++;
        if (yk > 29)                      /* computing ln(k!) with    */
                {                              /* Stirling's Approximation */
                 ry= 1.0/yk;
                 ry2=ry*ry;
                 fk = (yk + 0.5)*log(yk)-yk+C0+ry*(C1+ry2*(C3+ry2*(C5+ry2*C7)));
                }
        else
                {
                 k = yk;                      /* ln(k!) out of table logfak */
                 fk = logfak[k];              /* applied for k<=29          */
                }
        if (ym > 29)                     /* ln(m!) analogous to ln(k!) */
                {
                 ry= 1.0/ym;
                 ry2=ry*ry;
                 fm = (ym + 0.5)*log(ym)-ym+C0+ry*(C1+ry2*(C3+ry2*(C5+ry2*C7)));
                }
        else fm = logfak[m];
        q = ym * g - fm;
        h = (a - yk) * exp(0.5 * ((yk - ym) * g + fm - fk) + logfak[2L]);
                 }
         }
 if (my<5.5)                              /* Inversion/Chop-down */
         {
          double pk;

          pk=p0;
          k=0;
          u=drand(seed);
          while (u>pk)
          {
                ++k;
                if (k>b)
                        {
                         u=drand(seed);
                         k=0;
                         pk=p0;
                        }
                else
                        {
                         u-=pk;
                         pk*=my/k;
                        }
         }
         return(k);
        }
 for(;;)                                 /* Ratio of Uniforms */
        {
         do
                {
                 u = drand(seed);
                 x = a + h * (drand(seed) - 0.5) / u;
                }
         while ((x < 0) || (x >= b));      /* check, if x is valid candidate */
         k = (long int)(x);
         yk = k;
         if (k > 29)
                 {
        ry=1.0/yk;                               /* lf - ln(f) */
        lf = g * yk - (0.5 + yk) * log(yk) + yk - C0 - C1*ry - q; /* lower bound
*/
        if (lf >= u * (4.0 - u) - 3.0) break;    /* quick acceptance */
        ry2=ry*ry;
        lf -=ry*ry2*(C3+ry2*(C5+ry2*C7));        /* more exact ln(f) */
                 }
         else
                 {
        lf = yk * g - logfak[k] - q;            /* exact ln(f) using table */
        if (lf >= u * (4.0 - u) - 3.0) break;   /* to compute ln(k!)       */
                 }
         if (u * (u - lf) <= 1.0)               /* u*(u-lf)>1 - quick rejection */
                 if (2.0 * log(u) <= lf) break;      /* final acceptance */
        }
return(k);
}