double bnldev(double pp, int n, int &idum)
{
	int j,en,nold,bnl;
	double p,am,g,t,oldg,pc,plog,pclog,pold,sq,y,em,pi;
    pi = 3.141592654;
	nold=0;
	pold=0.0;
    if (pp <= 0.5)
	    p = pp;
    else
        p = 1.0 - pp;
    am = n * p;
    if (n < 25)
	{
        bnl = 0;
        for (j = 1; j<=n; j++)
            if (ran1(idum) < p) bnl = bnl + 1;
    }
    else
	{
		if (am < 1.0)
		{
        g = exp(-am);
        t = 1.0;
        for (j = 0;j<=n;j++)
			{
				t = t * ran1(idum);
				if (t < g) exit(1);
			}
        if (t >= g) j = n;
        bnl = j;
		}
        else
		{
             if (n!= nold)
					{
					en = n;
					oldg = gammln(en + 1.0);
					nold = n;
					}
			if (p != pold)
			{
	        pc = 1.0 - p;
            plog = log(p);
            pclog = log(pc);
            pold = p;
			}
			sq = sqrt(2.0 * am * pc);
			do
			{
				do
				{
				y = tan(pi * ran1(idum));
				em = sq * y + am;
				}
				while ((em < 0) || (em >= en + 1.0));
				em = int(em);
				t = en - em;
				t = exp(oldg - gammln(em + 1.0) - gammln(t + 1.0) + em * plog + t * pclog);
				t = 1.2 * sq * (1.0 + y *y) * t;
			}
			while (ran1(idum) > t);
        bnl = int(em);
		}
	}
    if (p != pp) bnl = n - bnl;
    t=bnl;
	return t;
}