我有一个有理函数：f(x) = P(x)/Q(x)。
例如：
f(x) = (5x + 3)/(1-x^2)
因为 f(x) 是生成函数，所以可以写成：
f(x) = a0 + a1*x + a2*x² + ... + a_n*x^n + ... = P(x)/Q(x)
如何使用 sympy 找到生成函数 f(x)（即 a_n）的第 nth 项？
如果 Sympy 中没有这样的实现，我也很想知道这是否在其他包中实现，例如 Maxima。
感谢您的帮助。
【问题讨论】：
要得到有理形式的生成函数a_n的通式，可以使用SymPy的rational_algorithm。
例如：
from sympy import simplify
from sympy.abc import x, n
from sympy.series.formal import rational_algorithm
f = (5*x + 3)/(1-x**2)
func_n, independent_term, order = rational_algorithm(f, x, n, full=True)
print(f"The general formula for a_n is {func_n}")
for k in range(10):
print(f"a_{k} = {simplify(func_n.subs(n, k))}")
输出：
The general formula for a_n is (-1)**(-n - 1) + 4
a_0 = 3
a_1 = 5
a_2 = 3
a_3 = 5
a_4 = 3
a_5 = 5
a_6 = 3
a_7 = 5
a_8 = 3
a_9 = 5
【讨论】：
您可以取第 k 个导数并将 0 替换为 x 并除以 factorial(k)：
>>> f = (5*x + 3) / (1-x**2)
>>> f.diff(x, 20).subs(x, 0)/factorial(20)
3
参考here 谈到了有理生成函数。寻找重复出现时，您可以使用微分很快地看到模式：
[f.diff(x,i).subs(x,0)/factorial(i) for i in range(6)]
[3, 5, 3, 5, 3, 5]
【讨论】：
适应this post的方法，你可以尝试以下：
from sympy import *
from sympy.abc import x
f = (5*x + 3) / (1-x**2)
print(f.series(n=20))
k = 50
coeff50 = Poly(series(f, x, n=k + 1).removeO(), x).coeff_monomial(x ** k)
print(f"The coeffcient of x^{k} of the generating function of {f} is {coeff50}")
# to get the first 100 coeffcients (reversing the list to get a[0] the
# coefficient of x**0 etc.):
a = Poly(series(f, x, n=100).removeO(), x).all_coeffs()[::-1]
输出：
3 + 5*x + 3*x**2 + 5*x**3 + 3*x**4 + 5*x**5 + 3*x**6 + 5*x**7 + 3*x**8 + 5*x**9 + 3*x**10 + 5*x**11 + 3*x**12 + 5*x**13 + 3*x**14 + 5*x**15 + 3*x**16 + 5*x**17 + 3*x**18 + 5*x**19 + O(x**20)
The coeffcient of x^50 of the generating function of (5*x + 3)/(1 - x**2) is 3
按照 Cut The Knot 的 this 示例，该方法可用于找出使用 1 美分、5 美分、10 美分、25 美分和 50 美分硬币支付金额 n 的方式数量。
f = 1/((1 - x)*(1 - x**5)*(1 - x**10)*(1 - x**25)*(1 - x**50))
a = Poly(series(f, x, n=101).removeO(), x).all_coeffs()[::-1]
print(a[50]) # there are 50 ways to pay 50 cents
print(a[100]) # there are 292 ways to pay 100 cents
【讨论】：