#include<math.h>

#define min(x,y) (x < y) ? x : y

extern double flogfak(long k);
extern double drand(unsigned long *seed);

/******************************************************************
 *                                                                *
 *      Binomial Distribution - Ratio of Uniforms/Inversion       *
 *                                                                *
 ******************************************************************
 *                                                                *
 * Ratio of Uniforms method combined with Inversion for sampling  *
 * from Binomial distributions with parameters n (number of       *
 * trials) and p (probability of success).                        *
 * For  min(n*p,n*(1-p)) < 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  min(n*p,n*(1-p)) >=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. *
 * If the candidate k is near the mode and k>29 (respectively     *
 * n-k>29) f(k) is computed recursively starting at the mode m.   *
 * The algorithm is valid for np > 0, n > =1, p <= 0.5.           *
 * For p > 0.5, p is replaced by 1-p and k is replaced by n-k.    *
 *                                                                *
 ******************************************************************
 *                                                                *
 * FUNCTION:    - bruec samples a random number from the Binomial *
 *                distribution with parameters n and p . It is    *
 *                valid for    n * min(p,1-p) > 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: - flogfak(k)  ... log(k!) with long integer k     *
 *              - drand(seed) ... (0,1)-Uniform generator with    *
 *                unsigned long integer *seed.                    *
 *                                                                *
 * Implemented by R.Kremer 1990, revised by P.Busswald, July 1992 *
 ******************************************************************/

long bruec(unsigned long *seed, long n, double p)

{
 static long b,m,n_prev = -1;
 static double par,q1,p_prev=-1.0,np,a,h,g,r,t,r1,p0,ln2=0.69314718055994531;
 long i,k,k1,bh;
 double u,c,f,g1,g2,x,lf;


 if(n_prev!=n || p_prev!=p)                 /* Set_up */
        {
         n_prev=n;
         p_prev=p;
         par=min(p,1.0-p);
         q1=1.0-par;
         np=n*par;                                /*np=min(n*p,n*(1-p))*/
         if (np<5)
        {
         p0=exp(n*log(q1));               /* Set-up for Inversion */
         bh=(long int)(np+10.0*sqrt(np*q1));
         b=min(n,bh);                     /* safety-bound */
        }
         else
        {                                 /* Set-up for Ratio of Uniforms */
         m=(long int)(np+par);            /* mode */
         a=np+0.5;                        /* shift parameter */
         c = sqrt(2.0*a*q1);
         r=par/q1;
         t=(n+1)*r;
         r1=log(r);
         bh=(long int)(a+7.0*c);
         b=min(n,bh);                     /* safety-bound */
         g=flogfak(m) + flogfak(n-m);     /* binomial const. */
         k1=(long int)(a-c);
         x=(a-k1-1.0)/(a-k1);
         if((n-k1)*par*x*x > (k1+1)*q1) k1++;           /* h=2*s */
         h=(a-k1)*exp(.5*((k1-m)*r1+g-flogfak(k1)-flogfak(n-k1))+ln2);
        }
         }
 if (np<5)                                /* Inversion/Chop-down */
         {
          double pk;

          k=0;
          pk=p0;
          u=drand(seed);
          while (u>pk)
                {
                 ++k;
                 if (k>b)
                         {
                u=drand(seed);
                k=0;
                pk=p0;
                         }
                 else
                         {
                u-=pk;
                pk=(double)(((n-k+1)*par*pk)/(k*q1));
                         }
                }
          return ((p>0.5) ? n-k:k);
         }
 for (;;)                                  /* Ratio of Uniforms */
         {
          do
                 {
        u=drand(seed);
        x=a+h*(drand(seed)-0.5)/u;
                 }
          while (x < 0 || ((k=(long int)x) > b)); /* check, if k is valid
candidate */
          if ((labs(m-k)<=15) && ((k>29)||(n-k>29)) )
                 {
        f=1.0;                       /* compute f(k) recursively */
        if (m<k)
          {
                for (i=m;i<k;) f*=t/(double)++i-r;           /* f - f(k) */
                if (u*u <= f) break;                         /* u^2<=f   */
          }
        else
          {
                for (i=k;i<m;) f*=t/(double)++i-r;          /* f - 1/f(k) */
                if (u*u*f <= 1.0) break;                    /* u^2<=f     */
          }
        }
          else
                 {
        lf=(k-m)*r1+g-flogfak(k)-flogfak(n-k);       /* lf - ln(f(k)) */
        if ( u * (4.0 - u) - 3.0 <= lf) break;       /* lower squeeze */
        if (u*(u-lf) <= 1.0)                         /* upper squeeze */
                if (2.0*log(u) <= lf) break;           /* final acceptance */
                 }
        }
 return((p > 0.5) ? n-k : k);
}