📄 bnldev.cpp
字号:
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;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -