深挖26年赣州一模压轴：用概率dp+多项式求通项，oier爱看的一集~

hans774882968
3 阅读4分钟

引言

本文视频版传送门

视频封面：

260320-封面.png

最近赣州一模的压轴题，大家都说是老题、套路题，只不过计算量超大！但如果我告诉你，这题居然可以用概率dp+多项式减少计算量呢？又是OIER最喜欢的一集，一起来领略概率dp和组合数学的魅力！

题干：

现有一种不断分裂的M细胞，在每个分裂周期中，一个M细胞以 23\frac{2}{3} 的概率分化成一个新的M细胞，以 13\frac{1}{3} 的概率分裂成两个新的M细胞，分裂后原来的M细胞消失，新的M细胞在下一个分裂周期里会继续分裂。设初始状态下有1个M细胞， nn 个分裂周期后，M细胞的数目为 XnX_n

  1. 求 X2X_2 的分布列和数学期望
  2. 求概率 P(Xn=2)P(X_n = 2)
  3. 证明： P(Xn=3)<65(23)n−1P(X_n = 3) < \frac{6}{5}(\frac{2}{3})^{n-1}

本文 52pojie： www.52pojie.cn/thread-2101…

本文 博客园： www.cnblogs.com/hans77/arti…

本文 juejin： juejin.cn/post/762446…

作者：hans774882968以及hans774882968以及hans774882968以及hans77

（1）求 X2X_2 的分布列和数学期望

X2X_2 是分裂两个周期后的细胞个数，故 X2X_2 的可能取值为 1, 2, 3, 4

  • P(X2=1)=P(X1=1)⋅23=23⋅23=49P(X_2 = 1) = P(X_1 = 1) \cdot \frac{2}{3} = \frac{2}{3} \cdot \frac{2}{3} = \frac{4}{9}
  • P(X2=2)=P(X1=1)⋅13+P(X1=2)⋅C21⋅23⋅13=23⋅13+13⋅49=1027P(X_2 = 2) = P(X_1 = 1) \cdot \frac{1}{3} + P(X_1 = 2) \cdot C_2^1 \cdot \frac{2}{3} \cdot \frac{1}{3} = \frac{2}{3} \cdot \frac{1}{3} + \frac{1}{3} \cdot \frac{4}{9} = \frac{10}{27}
  • P(X2=3)=P(X1=2)⋅C21⋅23⋅13=13⋅49=427P(X_2 = 3) = P(X_1 = 2) \cdot C_2^1 \cdot \frac{2}{3} \cdot \frac{1}{3} = \frac{1}{3} \cdot \frac{4}{9} = \frac{4}{27}
  • P(X2=4)=P(X1=2)⋅(13)2=13⋅19=127P(X_2 = 4) = P(X_1 = 2) \cdot (\frac{1}{3})^2 = \frac{1}{3} \cdot \frac{1}{9} = \frac{1}{27}

分布列：

X2X_21234
PP49\frac{4}{9}1027\frac{10}{27}427\frac{4}{27}127\frac{1}{27}

数学期望 E[X2]=1⋅49+2⋅1027+3⋅427+4⋅127=169E[X_2] = 1 \cdot \frac{4}{9} + 2 \cdot \frac{10}{27} + 3 \cdot \frac{4}{27} + 4 \cdot \frac{1}{27} = \frac{16}{9}

（2）求概率 P(Xn=2)P(X_n = 2) - 法1：不用递推式

不懂概率dp+全概率公式（高中通常称为马尔可夫链）能做吗？能的

Xn=2X_n = 2 表示n个细胞周期中，某个周期有1个M细胞分裂成两个，其他周期都分化成1个。故枚举所有可能的一分为二的时刻 k=1…nk=1 \dots n ：

  1. 前k-1个周期，只有一个细胞且都分化为1个，概率 (23)k−1(\frac{2}{3})^{k-1}
  2. 第k个周期，细胞一分为二，概率 13\frac{1}{3}
  3. 后n-k个周期，两个细胞每次都分化为1个，概率 (23)2(n−k)(\frac{2}{3})^{2(n-k)}

所以总概率

P(Xn=2)=∑k=1n(23)k−113(23)2(n−k)=∑k=1n13(23)2n−k−1=13(23)2n−1∑k=1n(32)k=(23)2n−1((32)n−1)=(23)n−1(1−(23)n)\begin{align} P(X_n = 2)&=\sum_{k=1}^{n} (\frac{2}{3})^{k-1}\frac{1}{3}(\frac{2}{3})^{2(n-k)}=\sum_{k=1}^{n} \frac{1}{3}(\frac{2}{3})^{2n-k-1} \\ &=\frac{1}{3}(\frac{2}{3})^{2n-1}\sum_{k=1}^{n}(\frac{3}{2})^{k}=(\frac{2}{3})^{2n-1}((\frac{3}{2})^n-1)=(\frac{2}{3})^{n-1}(1-(\frac{2}{3})^n) \end{align}

（3）证明： P(Xn=3)<65(23)n−1P(X_n = 3) < \frac{6}{5}(\frac{2}{3})^{n-1}

[!note] 想法 如果能求出通项公式，这题就很好解决。但这似乎不简单！那就考虑子问题的性质，求递推式吧！

记 P(Xn=m)P(X_{n}=m) 为 dp[m,n]dp[m,n] 。设第n-1个周期后有k个细胞，则：

  • 第n个周期需要在k个细胞中选m-k个来一分为二，剩下 2k−m2k-m 个细胞分化为一个，方案数 Ckm−kC_{k}^{m-k} 。而为了能做这个选取动作，需要 k≥m−k  ⟹  k≥⌈m2⌉k \geq m-k \implies k \geq \lceil \frac{m}{2} \rceil
  • 边界情况：0个周期后只有1个细胞，故 dp[m,0]=[m==1]dp[m,0]=[m==1] ，即n=0时值为1，其他情况值为0

于是由全概率公式：

{dp[m,0]=[m==1]dp[m,n]=∑k=⌈m/2⌉mdp[k,n−1]⋅C(k,m−k)⋅(13)m−k⋅(23)2k−m\begin{cases} dp[m,0]=[m == 1] \\ dp[m,n]=\sum_{k=\lceil m/2 \rceil}^{m} dp[k,n-1] \cdot C(k, m-k) \cdot \left(\frac{1}{3}\right)^{m-k} \cdot \left(\frac{2}{3}\right)^{2k-m} \end{cases}

以m=2为例：

dp[2,n]=C11dp[1,n−1]⋅(13)1⋅(23)0+C20dp[2,n−1]⋅(13)0⋅(23)2=13dp[1,n−1]+49dp[2,n−1]=13⋅(23)n−1+49dp[2,n−1]\begin{align} dp[2,n]&=C_{1}^1dp[1,n-1] \cdot (\frac{1}{3})^{1} \cdot (\frac{2}{3})^{0}+C_{2}^0dp[2,n-1] \cdot (\frac{1}{3})^{0} \cdot (\frac{2}{3})^{2} \\ &=\frac{1}{3}dp[1,n-1]+\frac{4}{9}dp[2,n-1] \\ &=\frac{1}{3} \cdot (\frac{2}{3})^{n-1}+\frac{4}{9}dp[2,n-1] \end{align}

再以m=3为例：

dp[3,n]=C21dp[2,n−1]⋅(13)1⋅(23)1+C30dp[3,n−1]⋅(13)0⋅(23)3=49dp[2,n−1]+827dp[3,n−1]\begin{align} dp[3,n]&=C_{2}^1dp[2,n-1] \cdot (\frac{1}{3})^{1} \cdot (\frac{2}{3})^{1}+C_{3}^0dp[3,n-1] \cdot (\frac{1}{3})^{0} \cdot (\frac{2}{3})^{3} \\ &=\frac{4}{9}dp[2,n-1]+\frac{8}{27}dp[3,n-1] \end{align}

代入上一问求出的 dp[2,n]=(23)n−1(1−(23)n)dp[2,n]=(\frac{2}{3})^{n-1}(1-(\frac{2}{3})^n) 得： dp[3,n]=(23)n(1−(23)n−1)+827dp[3,n−1]dp[3,n]=(\frac{2}{3})^{n}(1-(\frac{2}{3})^{n-1})+\frac{8}{27}dp[3,n-1]

问题转化为：用这条递推式求出通项公式

（3）求数列通项公式-法1：将递推式转化为累加形式

我们希望通过两边同乘 f(n)f(n) ，使得递推式能够累加，即

{f(n)dp[3,n]=⋯+827f(n)dp[3,n−1]f(n−1)dp[3,n−1]=⋯+827f(n−1)dp[3,n−2]…,f(1)dp[3,1]=0+827f(1)dp[3,0]\begin{cases} f(n)dp[3,n]=\dots+\frac{8}{27}f(n)dp[3,n-1] \\ f(n-1)dp[3,n-1]=\dots+\frac{8}{27}f(n-1)dp[3,n-2] \\ \dots, \quad f(1)dp[3,1]=0+\frac{8}{27}f(1)dp[3,0] \end{cases}

里满足 827f(n)=f(n−1), n≥2\frac{8}{27}f(n)=f(n-1),\ n \geq 2 。令 827f(1)=f(0)=1\frac{8}{27}f(1)=f(0)=1 得 f(n)=(32)3nf(n)=(\frac{3}{2})^{3n} 。于是累加得：

(32)3ndp[3,n]−dp[3,0]=∑k=1n(32)3k(23)k(1−(23)k−1)=∑k=1n(32)2k−(32)k+1=9454((32)2n−1)−9412((32)n−1)=95((32)2n−1)−92((32)n−1)=95(32)2n−92(32)n+2710\begin{align} (\frac{3}{2})^{3n}dp[3,n]-dp[3,0]&=\sum_{k=1}^{n} (\frac{3}{2})^{3k}(\frac{2}{3})^{k}(1-(\frac{2}{3})^{k-1})=\sum_{k=1}^{n} (\frac{3}{2})^{2k}-(\frac{3}{2})^{k+1} \\ &=\frac{\frac{9}{4}}{\frac{5}{4}} \left((\frac{3}{2})^{2n} - 1\right)-\frac{\frac{9}{4}}{\frac{1}{2}} \left((\frac{3}{2})^n - 1\right) \\ &= \frac{9}{5} \left((\frac{3}{2})^{2n} - 1\right)-\frac{9}{2} \left((\frac{3}{2})^n - 1\right)=\frac{9}{5} (\frac{3}{2})^{2n} -\frac{9}{2} (\frac{3}{2})^n + \frac{27}{10} \end{align}

于是 dp[3,n]=(23)3n(95(32)2n−92(32)n+2710)=2710(23)3n−92(23)2n+95(23)ndp[3,n]=(\frac{2}{3})^{3n}(\frac{9}{5} (\frac{3}{2})^{2n} -\frac{9}{2} (\frac{3}{2})^n + \frac{27}{10})=\frac{27}{10}(\frac{2}{3})^{3n}-\frac{9}{2}(\frac{2}{3})^{2n}+\frac{9}{5}(\frac{2}{3})^{n}

（3）求数列通项公式-法2：构造 dp[3,n]+f(n)dp[3,n]+f(n) 为等比数列

我们希望找到 f(n)f(n) ，使得 dp[3,n]+f(n)=827(dp[3,n−1]+f(n−1))dp[3,n]+f(n)=\frac{8}{27}(dp[3,n-1]+f(n-1)) ，展开并对比原递推式得：

f(n)−(23)3f(n−1)=(23)2n−1−(23)n\textcolor{orange}{ \boldsymbol{ f(n) - (\frac{2}{3})^3 f(n-1) = (\frac{2}{3})^{2n-1} - (\frac{2}{3})^n } }

右边是关于 x=(23)n∈(0,23]x=(\frac{2}{3})^n \in (0, \frac{2}{3}] 的二次函数 32x2−x\frac{3}{2}x^2-x ，故设 f(n)=A⋅(23)2n+B⋅(23)n+Cf(n) = A \cdot (\frac{2}{3})^{2n} + B \cdot (\frac{2}{3})^n+C ，于是方程左边

LHS=[A(23)2n+B(23)n+C]−(23)3[A(23)2n−2+B(23)n−1+C]=13A(23)2n+59B(23)n+1927C=A3x2+59Bx+1927C\begin{align} \text{LHS} &= \left[ A \left(\frac{2}{3}\right)^{2n} + B \left(\frac{2}{3}\right)^n+C \right] - \left(\frac{2}{3}\right)^3 \left[ A \left(\frac{2}{3}\right)^{2n-2} + B \left(\frac{2}{3}\right)^{n-1}+C \right] \\ &=\frac{1}{3} A \left(\frac{2}{3}\right)^{2n} + \frac{5}{9} B \left(\frac{2}{3}\right)^n+\frac{19}{27}C=\frac{A}{3}x^2+\frac{5}{9}Bx+\frac{19}{27}C \end{align}

对比系数得

{A3=3259B=−11927C=0  ⟹  {A=92B=−95C=0\begin{cases} \frac{A}{3}=\frac{3}{2} \\ \frac{5}{9}B=-1 \\ \frac{19}{27}C=0 \end{cases} \implies \begin{cases} A = \frac{9}{2} \\ B = -\frac{9}{5} \\ C=0 \end{cases}

设新数列 cn=dp[3,n]+f(n)=dp[3,n]+[92(23)2n−95(23)n]c_n = dp[3,n] + f(n)=dp[3,n]+\left[ \frac{9}{2} \left(\frac{2}{3}\right)^{2n} - \frac{9}{5} \left(\frac{2}{3}\right)^n \right] ，则它是首项 c0=2710c_{0}=\frac{27}{10} ，公比 827\frac{8}{27} 的等比数列。故

dp[3,n]=cn−f(n)=2710(23)3n−[92(23)2n−95(23)n]=2710(23)3n−92(23)2n+95(23)n\begin{align} dp[3,n]&=c_{n}-f(n) \\ &=\frac{27}{10} \left(\frac{2}{3}\right)^{3n} - \left[ \frac{9}{2} \left(\frac{2}{3}\right)^{2n} - \frac{9}{5} \left(\frac{2}{3}\right)^n \right] \\ &=\frac{27}{10}(\frac{2}{3})^{3n}-\frac{9}{2}(\frac{2}{3})^{2n}+\frac{9}{5}(\frac{2}{3})^{n} \end{align}

（3）给证明收尾（ P(Xn=3)<65(23)n−1P(X_n = 3) < \frac{6}{5}(\frac{2}{3})^{n-1} ）

前面已经求出通项公式

dp[3,n]=2710(23)3n−92(23)2n+95(23)ndp[3,n]=\frac{27}{10}(\frac{2}{3})^{3n}-\frac{9}{2}(\frac{2}{3})^{2n}+\frac{9}{5}(\frac{2}{3})^{n}

按原不等式想要的系数 65(23)n−1\frac{6}{5}(\frac{2}{3})^{n-1} 提取公因式：

P(Xn=3)=65(23)n−1(94(23)2n+1−154(23)n+1+1)=65(23)n−1(32(23)2n−52(23)n+1)\begin{align} P(X_{n}=3)&=\frac{6}{5}(\frac{2}{3})^{n-1}(\frac{9}{4}(\frac{2}{3})^{2n+1}-\frac{15}{4}(\frac{2}{3})^{n+1}+1) \\ &=\frac{6}{5}(\frac{2}{3})^{n-1}(\frac{3}{2}(\frac{2}{3})^{2n}-\frac{5}{2}(\frac{2}{3})^{n}+1) \end{align}

记 x=(23)n∈(0,23]x=(\frac{2}{3})^n \in (0, \frac{2}{3}] 则得 0≤32x2−52x+1<10 \leq \frac{3}{2}x^2-\frac{5}{2}x+1 < 1 ，证毕

扩展：求 P(Xn=m), m≥1P(X_{n}=m),\ m \geq 1 （即 dp[m,n]dp[m,n] ）的通项公式

题目做完了，但我们的探索才刚刚开始！这题只研究了m为1到3的情况。所以接下来我们推广到一般的m，求通项公式！再次回顾之前求出的递推式，我们知道：

{dp[m,0]=[m==1]dp[m,n]=∑k=⌈m/2⌉mdp[k,n−1]⋅C(k,m−k)⋅(13)m−k⋅(23)2k−m\begin{cases} dp[m,0]=[m == 1] \\ dp[m,n]=\sum_{k=\lceil m/2 \rceil}^{m} dp[k,n-1] \cdot C(k, m-k) \cdot \left(\frac{1}{3}\right)^{m-k} \cdot \left(\frac{2}{3}\right)^{2k-m} \end{cases}

用sympy符号计算的能力可以直接求解。于是有以下提示词：

大佬，你是一名数学科研工作者，精通概率论。有以下问题：（问题复述略）我希望写Python代码求解 P(Xn=m)P(X_n=m) 的表达式（m为常数）。容易知道这个表达式只与n有关，所以考虑路线：先求出递推式的各个系数，再代入 P(Xn=1)P(X_n=1) 的关于n的表达式，依次算出 P(Xn=2…m)P(X_n=2 \dots m) 。提示： P(Xn−1=t)P(X_{n-1}=t) （t为常数）也是关于n的式子，可以由 P(Xn=t)P(X_{n}=t) 的式子调用.subs()得到。注意：

  1. 需要使用Fraction表示精确结果
  2. 不要写和这个任务无关的代码
  3. 代码里不写注释，但要在外部说明这段代码做了什么

大佬，你代码的sol = rsolve(eq, y(n), {y(0): 0})这一行报错：AttributeError: 'NoneType' object has no attribute 'subs'。请帮我看看是哪里有bug

from sympy import symbols, Function, rsolve, simplify, binomial, Rational, expand


def verify_expression(ans, rhs_list, n_range=range(1, 10)):
    if not isinstance(rhs_list, list):
        rhs_list = [rhs_list]
    for rhs in rhs_list:
        delta_simp = simplify(ans - rhs)
        # 直接 delta_simp == 0 返回 False ，所以用数值解来验证
        for i in n_range:
            val = delta_simp.subs(n, i)
            if val != 0:
                raise ValueError(f"n={i} 时验证失败: ans={ans}, rhs={rhs}, diff={val}")


n = symbols('n', integer=True)
m = 5
probs = {}
probs[1] = Rational(2, 3)**n

for k in range(2, m + 1):
    y = Function('y')
    rhs = 0
    for j in range((k + 1) // 2, k):
        s = k - j
        coeff = binomial(j, s) * Rational(1, 3)**s * Rational(2, 3)**(j - s)
        term = probs[j].subs(n, n - 1) * coeff
        rhs += term

    rhs = simplify(expand(rhs))

    coeff_self = Rational(2, 3)**k

    eq = y(n) - coeff_self * y(n - 1) - rhs
    sol = rsolve(eq, y(n), {y(0): 0})

    if sol is None:
        raise ValueError(f"rsolve failed for k={k}")

    probs[k] = simplify(sol)

for k in range(1, m + 1):
    print(probs[k])

verify_expression(probs[2], Rational(2, 3)**(n - 1) * (1 - (Rational(2, 3))**n))
verify_expression(probs[3], [
    Rational(6, 5) * (Rational(2, 3))**(n - 1) * (1 - (Rational(2, 3))**(n - 1)) * (1 - (Rational(2, 3))**n),
    Rational(27, 10) * Rational(2, 3)**(3 * n) - Rational(9, 2) * Rational(2, 3)**(2 * n) + Rational(9, 5) * Rational(2, 3)**n
])
verify_expression(probs[4], 27 * Rational(2, 3)**n * (-153 * Rational(2, 3)**(n * 3) + 342 * Rational(2, 3)**(n * 2) - 247 * Rational(2, 3)**n + 58) / 760)

运行结果：

2: 3*(-4**n + 6**n)/(2*9**n)
3: 9*(-5*108**n + 2*162**n + 3*72**n)/(10*243**n)
4: 27*(58*118098**n - 153*34992**n + 342*52488**n - 247*78732**n)/(760*177147**n)
5: 81*(7209*1377495072**n - 19890*243**n*54**(4*n) + 20007*3099363912**n - 8710*4649045868**n + 1384*6973568802**n)/(49400*10460353203**n)

能求出更简洁的形式吗——概率dp求多项式系数

回顾之前求出的m=1到3情况的通项公式，我们知道

{dp[1,n]=(23)ndp[2,n]=(23)n−1(1−(23)n)dp[3,n]=2710(23)3n−92(23)2n+95(23)n\begin{cases} dp[1,n]=(\frac{2}{3})^n \\ dp[2,n]=(\frac{2}{3})^{n-1}(1-(\frac{2}{3})^n) \\ dp[3,n]=\frac{27}{10}(\frac{2}{3})^{3n}-\frac{9}{2}(\frac{2}{3})^{2n}+\frac{9}{5}(\frac{2}{3})^{n} \end{cases}

如果令 x=(23)nx=(\frac{2}{3})^n ，则它们都是关于x的多项式。并且 m=4m=4 时的通项公式27*(58*118098**n - 153*34992**n + 342*52488**n - 247*78732**n)/(760*177147**n)中，118098, 34992等数的质因数都只有2和3。我们猜测：

dp[m,n]dp[m,n] 是关于x的m次多项式，且常数项为0。

设c[m][i]为 dp[m,n]dp[m,n] 的 xix^i 项的系数，即 dp[m,n]=∑i=0mc[m,i]xidp[m,n]=\sum_{i=0}^{m} c[m,i]x^i 。于是我们的任务转化为求c数组。

为了求c数组，我们要重新观察递推式：

dp[m,n]=∑k=⌈m/2⌉mdp[k,n−1]⋅C(k,m−k)⋅(13)m−k⋅(23)2k−mdp[m,n]=\sum_{k=\lceil m/2 \rceil}^{m} dp[k,n-1] \cdot C(k, m-k) \cdot \left(\frac{1}{3}\right)^{m-k} \cdot \left(\frac{2}{3}\right)^{2k-m}

我们知道 dp[k,n−1]dp[k,n-1] 是关于 x′=(23)n−1x'=(\frac{2}{3})^{n-1} 的多项式，那 x′x' 和x有什么关系呢？显然是 x′=32xx'=\frac{3}{2}x 。又因为 dp[m,n]=∑i=0mc[m,i]xidp[m,n]=\sum_{i=0}^{m} c[m,i]x^i 里n从未单独出现，只蕴含在x里，所以如果记 f(x)=dp[m,n]f(x)=dp[m,n] ，那么 dp[m,n−1]=f(32x)dp[m,n-1]=f(\frac{3}{2}x)

所以 dp[k,n−1]dp[k,n-1] 对 c[m,i]c[m,i] 的贡献就是 c[k,i]⋅(32)ic[k,i] \cdot (\frac{3}{2})^i ，即：

c[m,i]=∑k=⌈m/2⌉m(c[k,i]⋅(32)i)⋅C(k,m−k)⋅(13)m−k⋅(23)2k−mc[m,i]=\sum_{k=\lceil m/2 \rceil}^{m} (c[k,i] \cdot (\frac{3}{2})^i) \cdot C(k, m-k) \cdot \left(\frac{1}{3}\right)^{m-k} \cdot \left(\frac{2}{3}\right)^{2k-m}

但当 k=mk=m 时，我们发现这里出现了自引用：为了求 c[m,i]c[m,i] ，需要先知道 c[m,i]c[m,i]

别担心！在概率dp里这很常见。出现自引用的情况，说明需要解方程。假设 k<mk<m 对 c[m,i]c[m,i] 的总贡献已经求出是 TT ，把递推式里 k=mk=m 的情况单独拿出来：

c[m,i]=T+c[m,i]⋅(32)i(23)m  ⟹  c[m,i]=T1−(23)m−ic[m,i]=T+c[m,i] \cdot (\frac{3}{2})^i (\frac{2}{3})^{m} \implies c[m,i]=\frac{T}{1-(\frac{2}{3})^{m-i}}

现在 c[m,i]c[m,i] 能求了，但我们看分母，当 m=im=i 时分母变成0了！那怎么求 c[m,m]c[m,m] 呢？

dp[m,n]dp[m,n] 的常数项为0，所以 x=1x=1 （即 n=1n=1 ）时 dp[m,n]=0dp[m,n]=0 ，又 dp[m,n]=∑i=0mc[m,i]xi  ⟹  ∑i=0mc[m,i]=0dp[m,n]=\sum_{i=0}^{m} c[m,i]x^i \implies \sum_{i=0}^{m} c[m,i]=0 。这意味着，求出 c[m,1…m−1]c[m,1 \dots m-1] 就能求出 c[m,m]=−∑i=0m−1c[m,i]c[m,m]=-\sum_{i=0}^{m-1} c[m,i] 了~

怀揣着把套路题推广到一般情况的热情，我们便欣赏到了不一样的风景！

思考

  1. 设 x=(23)nx=(\frac{2}{3})^n ，证明 dp[m,n]dp[m,n] 是关于x的多项式
  2. 如何证明 c[m,0]=0c[m,0]=0 ？
from sympy import symbols, Rational, binomial, simplify

x = symbols('x')
m = 10
p, q = Rational(1, 3), Rational(2, 3)

coeffs = {1: {1: Rational(1)}}
ans_list = [None] * (m + 1)

for k in range(2, m + 1):
    coeffs[k] = {}
    for i in range(1, k):
        d_ki = Rational(0)
        for j in range(max((k + 1) // 2, i), k):
            s = k - j
            d_ki += binomial(j, s) * p**s * q**(j - s) * coeffs[j][i] * Rational(3, 2)**i
        coeffs[k][i] = simplify(d_ki / (Rational(1) - q**(k - i)))
    coeffs[k][k] = simplify(-sum(coeffs[k][i] for i in range(1, k)))

for k in range(1, m + 1):
    ans_list[k] = simplify(sum(coeffs[k][i] * x**i for i in range(1, k + 1)))
for k in range(1, m + 1):
    print(f'{k}: {ans_list[k]}')

运行结果：

2: 3*x*(1 - x)/2
3: 9*x*(3*x**2 - 5*x + 2)/10
4: 27*x*(-153*x**3 + 342*x**2 - 247*x + 58)/760
5: 81*x*(7209*x**4 - 19890*x**3 + 20007*x**2 - 8710*x + 1384)/49400
avatar
文章
阅读
粉丝