user

1
2
3
4
5
6
7
8
9
电路瑕疵

在一个工厂生产的n块集成电路上随机分布有k个瑕疵（一块集成电路上可能有任意数量的瑕疵，瑕疵之间相互独立）。

记p(k,n)为存在一块集成电路上面有至少3个瑕疵的概率。
例如，p(3,7) ≈ 0.0204081633。

求p(20 000, 1 000 000)，并将你的答案四舍五入到10位小数，即形式为0.abcdefghij。
请给我python代码实现.

kimi

我来分析这个问题。这是一个经典的概率问题，可以用容斥原理或补集来解决。

问题分析

求补集更简单：计算”所有电路板上瑕疵数都 ≤ 2”的概率，然后用 1 减去它。

设 $n = 10^6$ 块电路板，$k = 20000$ 个瑕疵。

每个瑕疵独立随机落在某块板上。我们需要计算没有板子有 ≥3 个瑕疵的概率。

数学推导

设 $X_i$ 为第 $i$ 块板上的瑕疵数，则 $(X_1, …, X_n) \sim \text{Multinomial}(k; \frac{1}{n}, …, \frac{1}{n})$

我们需要：$P(\forall i: Xi \leq 2) = \sum{\substack{a_1+…+a_n=k \ 0\leq a_i\leq 2}} \frac{k!}{a_1!…a_n!} \cdot n^{-k}$

这等于：$\frac{k!}{n^k} \cdot [x^k] \left(\sum_{j=0}^{2} \frac{x^j}{j!}\right)^n = \frac{k!}{n^k} \cdot [x^k] \left(1 + x + \frac{x^2}{2}\right)^n$

Python 实现

使用高精度计算（Python 的 decimal 模块）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
from decimal import Decimal, getcontext
import math

getcontext().prec = 50  # 设置高精度

def solve():
    n = 1000000
    k = 20000
    
    # 使用 decimal 进行高精度计算
    # 计算 log(P(所有板子瑕疵数 <= 2))
    
    # 方法：直接计算概率
    # P = (k! / n^k) * [x^k] (1 + x + x^2/2)^n
    
    # 使用动态规划计算系数，但 k=20000, n=1000000 太大
    # 改用泊松近似或更好的方法
    
    # 实际上：当 n 很大，每个板子的瑕疵数 ~ Poisson(lambda=k/n=0.02)
    # P(X_i <= 2) 对于 Poisson(0.02) 很精确
    
    lam = Decimal(k) / Decimal(n)  # 0.02
    
    # Poisson 概率
    # P(X = j) = e^(-lambda) * lambda^j / j!
    
    # 由于瑕疵之间独立？不，总数固定为 k，是多项分布
    
    # 但对于大 n，可以用 Poisson 近似
    # 实际上我们需要精确值
    
    # 改用：直接数值计算
    # log P = log(k!) - k*log(n) + log(系数)
    
    # 系数 [x^k](1+x+x^2/2)^n 可以用组合数学
    # 设选 t2 个板子有2个瑕疵，t1 个板子有1个，t0=n-t1-t2 有0个
    # 2*t2 + t1 = k, 即 t1 = k - 2*t2
    # t2 从 0 到 k//2
    
    # 系数 = sum_{t2=0}^{k//2} C(n, t2) * C(n-t2, k-2*t2) * (1/2)^t2
    
    log_answer = Decimal(0)
    
    # 计算 log(k!)
    log_k_fact = Decimal(0)
    for i in range(1, k+1):
        log_k_fact += Decimal(i).ln()
    
    log_n = Decimal(n).ln()
    
    # 计算系数
    max_t2 = k // 2
    
    # 我们需要计算 sum_{t2} exp(log_term)
    # 为避免下溢，先找到最大项，然后求和
    
    log_terms = []
    
    for t2 in range(max_t2 + 1):
        t1 = k - 2 * t2  # 有1个瑕疵的板子数
        
        if t1 < 0 or t1 > n - t2:
            continue
            
        # C(n, t2) * C(n-t2, t1) * (1/2)^t2
        # = n! / (t2! * t1! * (n-t2-t1)!) * (1/2)^t2
        
        log_term = Decimal(0)
        
        # log(n!)
        log_n_fact = Decimal(0)
        for i in range(1, n+1):
            log_n_fact += Decimal(i).ln()
        
        # 减去分母
        log_t2_fact = Decimal(0)
        for i in range(1, t2+1):
            log_t2_fact += Decimal(i).ln()
            
        log_t1_fact = Decimal(0)
        for i in range(1, t1+1):
            log_t1_fact += Decimal(i).ln()
            
        log_t0_fact = Decimal(0)  # t0 = n - t2 - t1
        t0 = n - t2 - t1
        for i in range(1, t0+1):
            log_t0_fact += Decimal(i).ln()
        
        log_term = log_n_fact - log_t2_fact - log_t1_fact - log_t0_fact
        log_term -= t2 * Decimal(2).ln()  # (1/2)^t2
        
        # 加上前面的系数
        log_term += log_k_fact - k * log_n
        
        log_terms.append((t2, log_term))
    
    # 找到最大值用于数值稳定性
    max_log = max(t[1] for t in log_terms)
    
    sum_exp = Decimal(0)
    for t2, lt in log_terms:
        sum_exp += (lt - max_log).exp()
    
    log_prob_no_ge3 = max_log + sum_exp.ln()
    
    # 概率 = 1 - exp(log_prob_no_ge3)
    prob_no_ge3 = log_prob_no_ge3.exp()
    ans = 1 - prob_no_ge3
    
    return ans

# 优化版本，使用更高效的阶乘计算
from decimal import Decimal, getcontext
import math

getcontext().prec = 100

def log_factorial(n):
    """使用 Stirling 近似或精确计算"""
    if n < 10000:
        # 小数值精确计算
        res = Decimal(0)
        for i in range(2, n+1):
            res += Decimal(i).ln()
        return res
    else:
        # Stirling 近似
        n = Decimal(n)
        return n * n.ln() - n + (2 * Decimal(math.pi) * n).ln() / 2 + 1 / (Decimal(12) * n)

def solve_fast():
    n = 1000000
    k = 20000
    
    lam = Decimal(k) / Decimal(n)  # 0.02
    
    # 对于多项分布，用 Poisson 近似非常精确（当 n 大，lambda 小时）
    # 但我们还是尝试精确计算
    
    # 实际上，由于 k << n，可以用包含-排除或直接计算
    
    # 精确方法：枚举有2个瑕疵的板子数 t2
    # 有1个瑕疵的板子数 t1 = k - 2*t2
    
    max_t2 = k // 2
    
    # 预计算 log 阶乘
    # 使用 Gamma 函数或直接累加
    
    # 为了精度，我们分段计算
    
    log_terms = []
    
    for t2 in range(0, max_t2 + 1, 1):  # 步长可以优化，找主要贡献
        t1 = k - 2 * t2
        if t1 < 0 or t1 + t2 > n:
            continue
        
        # 多项式系数部分：k! / (2!^t2 * 1!^t1 * 0!^t0) 其中 t0 = n - t1 - t2
        # 但还要乘上选择哪些板子的组合数
        
        # 总数： multinomial coefficient * (1/n)^k
        # = k! / (t2! * t1! * t0! * 2!^t2) * n! / (n-t2-t1)! / n^k  ？不对
        
        # 正确：先选 t2 个板子放2个瑕疵：C(n, t2)
        # 再从剩下选 t1 个放1个：C(n-t2, t1)
        # 然后分配瑕疵：k! / (2!^t2 * 1!^t1) 
        # 概率：(1/n)^k
        
        t0 = n - t1 - t2
        
        # log(C(n, t2) * C(n-t2, t1) * k! / (2^t2))
        log_val = log_factorial(n) - log_factorial(t2) - log_factorial(t1) - log_factorial(t0)
        log_val += log_factorial(k) - t2 * Decimal(2).ln()  # 除以 2^t2
        log_val -= k * Decimal(n).ln()  # 除以 n^k
        
        log_terms.append(log_val)
        
        # 只保留主要贡献的项（优化）
        if t2 > 100 and log_val < log_terms[0] - 50:  # 后续贡献极小
            break
    
    # 数值稳定求和
    max_log = max(log_terms)
    total = sum((lt - max_log).exp() for lt in log_terms)
    log_prob = max_log + total.ln()
    
    prob_no_ge3 = log_prob.exp()
    result = 1 - prob_no_ge3
    
    return result

# 更精确的版本，使用 mpmath
try:
    from mpmath import mp, log, exp, factorial
    
    mp.dps = 50  # 50位小数精度
    
    def solve_mpmath():
        n = 1000000
        k = 20000
        
        # 使用 Poisson 近似：当 n 很大，p=k/n 很小，多项分布 ≈ 独立 Poisson
        # 但严格来说，我们需要多项分布的精确概率
        
        # 对于多项分布 (k; 1/n, ..., 1/n)
        # P(all X_i <= 2) = sum_{t2=0}^{floor(k/2)} C(n,t2) * C(n-t2, k-2*t2) * k! / (2^t2) / n^k
        
        max_t2 = k // 2
        
        terms = []
        
        for t2 in range(max_t2 + 1):
            t1 = k - 2 * t2
            if t1 < 0 or t1 > n - t2:
                continue
            
            t0 = n - t1 - t2
            
            # log term
            log_term = log(factorial(n)) - log(factorial(t2)) - log(factorial(t1)) - log(factorial(t0))
            log_term += log(factorial(k)) - t2 * log(2)
            log_term -= k * log(n)
            
            terms.append(log_term)
        
        # 找到最大项
        max_log = max(terms)
        
        # 稳定求和
        s = mp.mpf('0')
        for lt in terms:
            s += exp(lt - max_log)
        
        log_prob_no_ge3 = max_log + log(s)
        prob_no_ge3 = exp(log_prob_no_ge3)
        
        return 1 - prob_no_ge3
    
    ans = solve_mpmath()
    print(f"Result: {ans}")
    print(f"Rounded: {mp.nstr(ans, 12)}")
    
except ImportError:
    print("mpmath not available, using decimal...")
    # 使用 decimal 回退
    ans = solve_fast()
    print(f"Result: {ans}")

关键优化：由于 $\lambda = k/n = 0.02$ 很小，有瑕疵的板子很少，主要贡献来自 $t_2$ 很小的项（即很少有板子有2个瑕疵）。实际上 $t_2$ 的期望值约为 $k \cdot \lambda / 2 = 200$，标准差也小，所以只需计算 $t_2 \in [0, 400]$ 左右即可。

如果你需要，我可以运行这段代码给出精确数值，或者你可以直接复制运行（需要安装 mpmath：pip install mpmath）。

fun fact

  • 高峰时段算力不足，已切换至 K2.5 快速，升级会员畅用思考模型 去升级

kimi永远说给你用k2.5快速,凌晨三点使用都会说给你快速模型用的.

这个题要考虑正难则反,我没想到怎么做,kimi提醒我的.

看了大佬博客发现要枚举瑕疵数.
因为最多有两个瑕疵,那么枚举2个瑕疵的工件的个数,1和0的都可以被确定,那么工件排列能算出来.因为每个瑕疵独立,那么全排列再乘上一个公式就是最终概率.