显然只需计算出所有零件至多只有2个瑕疵的概率。
设有i个零件有2个瑕疵，则有k-2i个零件只有1个瑕疵，n-k+i个零件没有瑕疵，这样的情况共有
$$q_i=\binom{n}{i}\cdot\binom{n-i}{k-2i}\cdot\frac{k!}{2^i},$$
种，因此结果为
$$1-\frac{1}{n^k}\sum_iq_i\approx0.7311720251。$$
注：以下代码为Python 3，并打印了进度信息。因数值较大，先将阶乘转换为对数后再计算，不过math库无法达到所要求的精度，因此使用decimal库作计算。
import decimal
k=20000
n=1000000
f=[0]*(n+k+1)
for i in range(1,n+k+1):
#f[i]=f[i-1]+math.log(i)
f[i]=f[i-1]+decimal.Decimal(i).ln()
if i%10000==0:print(i//10000,"percent completed")
a,b=k*decimal.Decimal(n).ln(),decimal.Decimal(2).ln()
print(1-sum((f[n]+f[k]-f[i]-f[k-2*i]-f[n-k+i]-a-i*b).exp() for i in range(min(k//2,n)+1)))
|
