ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

SymPy 离散具体数学(Concrete):超几何项判定、Gosper 求和与序列猜想工具实战指南

SymPy 离散具体数学(Concrete):超几何项判定、Gosper 求和与序列猜想工具实战指南 SymPy 离散具体数学Concrete超几何项判定、Gosper 求和与序列猜想工具实战指南【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympySymPy 的concrete子模块源码位于 sympy/concrete是纯 Python 计算机代数系统中面向“具体数学”Concrete Mathematics的离散符号计算工具集它以超几何项hypergeometric term为理论核心向上提供Sum/Product求和求积、Gosper 超几何求和算法向下提供序列识别guess与生成函数猜想工具。读完本文你将掌握如何判定一个序列是否为超几何项、如何用hypersimp得到相邻项之比的最简多项式商、如何用gosper_sum求封闭形式和、以及如何从一串有理数序列反推其递推关系与生成函数。本文内容以 doc/src/modules/concrete.rst 为骨架结合仓库内源码与测试如 test_gosper.py、test_guess.py展开所有示例均可直接在 Python 中复制运行。一、超几何项离散求和与递推的中心舞台在递推求解recurrence solving与求和summation中超几何项hypergeometric term占据核心地位。形式化地说超几何项是被一阶线性递推算子零化的序列给定序列a(n)若其相邻项之比a(n1)/a(n)是n的有理函数则称a(n)为超几何项。直观理解多项式、阶乘、组合数、上升/下降阶乘、伽马函数、指数函数等具体数学中的常见物种其相邻项比值都能化为n的多项式商因此都是超几何项。1.1 用is_hypergeometric快速判定SymPy 在Basic基类上提供了is_hypergeometric(k)方法实现见 sympy/core/basic.py它内部调用hypersimp(self, k)并判断结果是否为Nonefrom sympy import * n, k symbols(n,k) # 多项式当然是超几何项 (n**2 1).is_hypergeometric(n) # True # 具体数学中的常见物种 factorial(n).is_hypergeometric(n) # True binomial(n, k).is_hypergeometric(n) # True rf(n, k).is_hypergeometric(n) # True 上升阶乘 ff(n, k).is_hypergeometric(n) # True 下降阶乘 gamma(n).is_hypergeometric(n) # True (2**n).is_hypergeometric(n) # True需要注意二项式系数以及上升、下降阶乘对两个参数都是超几何的在另一个参数上同样成立binomial(n, k).is_hypergeometric(k) # True rf(n, k).is_hypergeometric(k) # True ff(n, k).is_hypergeometric(k) # True1.2 整数线性参数依然成立上述所有示例对n的整数线性参数依然成立这是求和算法能够处理形如factorial(2n)、binomial(3n1, k)这类项的关键factorial(2*n).is_hypergeometric(n) # True binomial(3*n1, k).is_hypergeometric(n) # True rf(n1, k-1).is_hypergeometric(n) # True ff(n-1, k1).is_hypergeometric(n) # True gamma(5*n).is_hypergeometric(n) # True (2**(n-7)).is_hypergeometric(n) # True1.3 非线性参数会使序列失去超几何性一旦参数变为非线性如n**2、n**3相邻项之比不再是有理函数判定结果即变为Falsefactorial(n**2).is_hypergeometric(n) # False (2**(n**3 1)).is_hypergeometric(n) # False实现细节从源码看is_hypergeometric对Piecewise表达式直接返回None见 basic.py即分段定义序列不在超几何判定范围之内。二、hypersimp把相邻项之比化简为最小次数多项式商如果不仅想知道是否是超几何项还想得到相邻项之比化简后的最简形式应使用hypersimp(f, k)函数定义于 sympy/simplify/simplify.py。其工作流程在源码 docstring 中清晰给出共三步尽可能把所有函数改写为伽马函数gamma形式把所有 gamma 改写为 gamma 与整数或绝对常数指数的上升阶乘的乘积化简嵌套分式与幂若结果恰为多项式商则约分降低其总次数。若f(k)是超几何项函数返回最小次数的多项式商f(k1)/f(k)否则返回None表示该序列不是超几何项 from sympy import hypersimp, factorial from sympy.abc import n hypersimp(factorial(2*n), n) 2*(n 1)*(2*n 1) hypersimp(factorial(n**2), n) # 返回 None空行第一行结果2*(n1)*(2*n1)正是factorial(2n2)/factorial(2n)的约分结果第二行因为相邻项之比不是有理函数返回None。源码补充hypersimp先用g f.subs(k, k1) / f构造相邻项之比再经rewrite(gamma)、expand_func、powsimp等化简最后用is_rational_function(k)判断是否为有理函数是则返回simplify约分结果否则返回None。同一文件中的hypersimilar(f, g, k)还提供超相似判定——两个项的商为k的有理函数时返回True在求解递推关系时很有用。三、Gosper 算法超几何求和的封闭形式Gosper 算法是超几何求和的经典算法用于计算形如s_n sum(f(k), (k, 0, n-1))的和式封闭形式其中f为不依赖于n的超几何项。其核心思想是寻找另一个超几何项g_n满足g_{n1} - g_n f_n即f的不定和分从而将求和问题转化为简单的代值计算。算法参考了 Marko Petkovsek、Herbert S. Wilf 与 Doron Zeilberger 的著作《A B》AK Peters, 1997, pp. 73–100。sympy.concrete.gosper模块提供三个逐层递进的函数sympy/concrete/gosper.py3.1gosper_normalGosper 正规形gosper_normal(f, g, n, polysTrue)把互素单变量多项式f(n)/g(n)改写为如下正规形f(n)/g(n) Z · (A(n)·C(n1)) / (B(n)·C(n))其中Z为任意常数A、B、C是n的首一多项式且满足三条互素性条件gcd(A(n), B(nh)) 1对所有自然数h成立、gcd(B(n), C(n1)) 1、gcd(A(n), C(n)) 1。这种有理分解是 Gosper 算法与差分方程求解的关键步骤也可用于判定两个超几何项是否相似。 from sympy.concrete.gosper import gosper_normal from sympy.abc import n gosper_normal(4*n5, 2*(4*n1)*(2*n3), n, polysFalse) (1/4, n 3/2, n 1/4)返回三元组(Z*A, B, C)。测试 test_gosper.py 同时验证了polysTrue返回Poly对象与polysFalse返回普通表达式两种模式结果一致。3.2gosper_term寻找不定和分gosper_term(f, n)对给定的超几何项f返回满足g_{n1} - g_n f_n的超几何项g_n。其内部流程为先用hypersimp求相邻项之比再用gosper_normal做有理分解之后求解一个关于待定系数的线性方程组源码中通过构造H A*x.shift(1) - B*x - C并solve系数实现。若f不是超几何项、或不可 Gosper 求和则返回None。 from sympy.concrete.gosper import gosper_term from sympy import factorial from sympy.abc import n gosper_term((4*n 1)*factorial(n)/factorial(2*n 1), n) (-n - 1/2)/(n 1/4)3.3gosper_sum直接得到封闭和gosper_sum(f, k)是对用户最友好的入口给定超几何项f计算g_n - g(0)其中g_{n1} - g_n f_n若和式无法表示为超几何项的封闭形式返回None。它接受两种调用形式定和(k, a, b)与不定和仅传符号k。 from sympy.concrete.gosper import gosper_sum from sympy import factorial from sympy.abc import n, k f (4*k 1)*factorial(k)/factorial(2*k 1) gosper_sum(f, (k, 0, n)) (-factorial(n) 2*factorial(2*n 1))/factorial(2*n 1) _.subs(n, 2) sum(f.subs(k, i) for i in [0, 1, 2]) True gosper_sum(f, (k, 3, n)) (-60*factorial(n) factorial(2*n 1))/(60*factorial(2*n 1)) _.subs(n, 5) sum(f.subs(k, i) for i in [3, 4, 5]) True注意到 docstring 与测试都用subs(n, 2) sum(...)做数值回验这是验证封闭形式正确性的标准做法。测试文件 test_gosper.py 还覆盖了更多经典结果 from sympy.concrete.gosper import gosper_sum from sympy import factorial, binomial from sympy.abc import k, n gosper_sum(1, (k, 0, n)) # n 1 gosper_sum(k, (k, 0, n)) # n*(n 1)/2 gosper_sum(k**2, (k, 0, n)) # n*(n 1)*(2*n 1)/6 gosper_sum(k**3, (k, 0, n)) # n**2*(n 1)**2/4 gosper_sum(2**k, (k, 0, n)) # 2*2**n - 1 gosper_sum(factorial(k), (k, 0, n)) # None不可超几何求和 gosper_sum(binomial(n, k), (k, 0, n)) # None实践要点factorial(k)与binomial(n, k)的定和返回None说明它们不满足 Gosper 可和性这正是为什么sum(binomial(n,k), k)需要借助其他机制如二项式定理/超几何恒等式。遇到None时应转而使用Sum对象或其数值求值能力而不是强行期望封闭形式。四、Sum、Product 与求和/求积入口函数concrete模块的类参考见 concrete.rst包含三个核心类sympy.concrete.summations.Sum含ExprWithIntLimits基类表示未求值的符号和式可通过.doit()求值sympy.concrete.products.Product表示未求值的符号连乘可通过.doit()求值sympy.concrete.expr_with_intlimits.ExprWithIntLimitsSum/Product共用的整数上下限抽象基类承载了换元、下限平移等公共逻辑。对应的函数级入口为summation与product定义于 sympy/concrete/summations.py 与 sympy/concrete/products.py两者语法与Integral一致都是f后跟(i, a, b)这样的元组且支持多重求和/求积重复传入多个符号元组 from sympy import summation, product, symbols, oo, log i, n, m symbols(i n m, integerTrue) # 求和计算失败时返回未求值的 Sum 对象 summation(2*i - 1, (i, 1, n)) n**2 summation(1/2**i, (i, 0, oo)) 2 summation(1/log(n)**n, (n, 2, oo)) # 无法封闭求值 Sum(log(n)**(-n), (n, 2, oo)) summation(i, (i, 0, n), (n, 0, m)) # 多重求和 m**3/6 m**2/2 m/3 # 无穷级数 from sympy import factorial from sympy.abc import x summation(x**n/factorial(n), (n, 0, oo)) exp(x) # 求积与 Sum 对称 i, k, m symbols(i k m, integerTrue) product(i, (i, 1, k)) factorial(k) product(m, (i, 1, k)) m**k从源码看summation的实现就是return Sum(f, *symbols, **kwargs).doit(deepFalse)product类似地构造Product并求值若无法求值则返回未求值对象因此掌握Sum/Product的doit机制即可理解这两个入口函数。Sum内部还集成了多种求值策略telescopic裂项、eval_sum多项式求和等见 summations.py 中的telescopic_direct/telescopic辅助函数Gosper 算法也是其超几何项求和的候选策略之一。五、序列猜想工具从数列反推公式与生成函数sympy.concrete.guess模块sympy/concrete/guess.py提供一组从若干项猜公式的工具核心函数guess改编自 Christian Krattenthaler 的 Mathematica 软件包Rate.m。整套工具在测试文件 test_guess.py 中有完整验证。5.1find_simple_recurrence识别线性递推find_simple_recurrence(v, AFunction(a), NSymbol(n))从若干个整数或有理数项中检测并返回递推关系。返回表达式中函数名默认为a主变量默认为n且最小下标恒为n不会是n-1、n-2等。其底层函数find_simple_recurrence_vector(l)返回长度为n的系数向量当发现n阶递推时若只返回[0]则说明未找到关系——注意该函数对二次无理数等特殊实数需谨慎使用源码 docstring 有明确警告。 from sympy.concrete.guess import find_simple_recurrence from sympy import fibonacci find_simple_recurrence([fibonacci(k) for k in range(12)]) -a(n) - a(n 1) a(n 2) # 自定义函数名与主变量 from sympy import Function, Symbol a [1, 1, 1] for k in range(15): a.append(5*a[-1]-3*a[-2]8*a[-3]) find_simple_recurrence(a, AFunction(f), NSymbol(i)) -8*f(i) 3*f(i 1) - 5*f(i 2) f(i 3) from sympy.concrete.guess import find_simple_recurrence_vector find_simple_recurrence_vector([fibonacci(k) for k in range(12)]) [1, -1, -1]5.2rationalize从浮点数识别有理数rationalize(x, maxcoeff10000)通过连分数从浮点值或mpmath.mpf识别有理数。算法在检测到超过阈值默认 10000的大部分商partial quotient时停止。与Fraction.from_decimal、mpmath.identify、nsimplify等方法不同它关注的是部分商的量级而非全局近似精度——如果该实数已知是有理数即使分母很大也能在默认参数下正确识别。 from sympy.concrete.guess import rationalize from mpmath import cos, pi rationalize(cos(pi/3)) 1/2 from mpmath import mpf rationalize(mpf(0.333333333333333)) 1/3 # 提高 maxcoeff 阈值可用作近似 rationalize(pi, maxcoeff250) 355/1135.3guess_generating_function猜生成函数guess_generating_function(v, XSymbol(x), types[all], maxsqrtn2)尝试为有理数序列v猜出生成函数返回一个字典键为生成函数类型名。目前实现了六种类型type形式定义ogff(x) Sum( a_k * x^k, k: 0..infinity )普通生成函数egff(x) Sum( a_k * x^k / k!, k: 0..infinity )指数生成函数lgff(x) Sum( (-1)^(k1) * a_k * x^k / k, k: 1..infinity )对数生成函数初始下标为 1hlgff(x) Sum( a_k * x^k / k, k: 1..infinity )双曲对数生成函数初始下标为 1lgdogff(x) d/dx log( Sum( a_k * x^k, k: 0..infinity ) )普通生成函数的对数导数lgdegff(x) d/dx log( Sum( a_k * x^k / k!, k: 0..infinity ) )指数生成函数的对数导数参数说明与使用要点types默认为[all]只关心部分类型时可传入类型列表以节省时间。注意丢弃某类型只是不为其做额外计算结果字典中仍可能包含该类型因为可以从其他类型轻松转换而来lgdogf与lgdegf在序列首项为 0 时不会被计算此时可先去掉前导零再重试maxsqrtn默认 2指定要测试的有理函数的 n 次方根的最大阶数用于识别生成函数为某有理函数的平方根等情形。 from sympy.concrete.guess import guess_generating_function as ggf ggf([k1 for k in range(12)], types[ogf, lgf, hlgf]) {hlgf: 1/(1 - x), lgf: 1/(x 1), ogf: 1/(x**2 - 2*x 1)} from sympy import sympify l sympify([3/2, 11/2, 0, -121/2, -363/2, 121]) ggf(l) {ogf: (x 3/2)/(11*x**2 - 3*x 1)} from sympy import fibonacci ggf([fibonacci(k) for k in range(5, 15)], types[ogf]) {ogf: (3*x 5)/(-x**2 - x 1)} from sympy import factorial ggf([factorial(k) for k in range(12)], types[ogf, egf, lgf]) {egf: 1/(1 - x)} ggf([k1 for k in range(12)], types[egf]) {egf: (x 1)*exp(x), lgdegf: (x 2)/(x 1)} # n 次方根检测对应 OEIS A108626 序列 ggf([1, 2, 5, 14, 41, 124, 383, 1200, 3799, 12122, 38919])[ogf] sqrt(1/(x**4 2*x**2 - 4*x 1))其中guess_generating_function_rational(v, XSymbol(x))是只处理有理生成函数的低层版本先用find_simple_recurrence_vector求分母q再按卷积关系求分子p。其返回(3*x 5)/(-x**2 - x 1)之类的分式当未找到时返回None。源码 docstring 同时提示它与 sympy/series/approximants.py 中的approximantsPadé 近似型序列逼近功能相关可互为补充。5.4guess从序列猜组合公式guess(l, allFalse, evaluateTrue, niter2, variablesNone)从一串有理数序列猜出闭式公式返回一个公式列表可能是多个等价结果。参数语义allFalse默认一旦某次迭代出结果即停止计算加速流程设为True则继续更多迭代可能返回更多可能与前序等价的公式evaluateTrue默认对结果中的连乘进行求值设为False则保留未求值的Product对象便于观察结构niter2迭代次数最大可取len(l)-1。迭代阶数越高结果越复杂第一次迭代返回多项式或有理函数第二次迭代返回上升阶乘及其倒数的乘积第三次迭代返回上升阶乘乘积的乘积依此类推。variablesNone返回公式默认包含符号i0, i1, i2, ...其中主变量是i0辅助变量为i1, i2, ...也可传入自定义符号列表长度应不小于niter主变量取列表第一个符号。 from sympy.concrete.guess import guess guess([1,2,6,24,120], evaluateFalse) [Product(i1 1, (i1, 1, i0 - 1))] from sympy import symbols r guess([1,2,7,42,429,7436,218348,10850216], niter4) i0 symbols(i0) [r[0].subs(i0,n).doit() for n in range(1,10)] [1, 2, 7, 42, 429, 7436, 218348, 10850216, 911835460]注意事项来自源码 docstring 与实现若序列除最后一项外含有 0guess直接返回空列表[]源码第 448 行检查any(a0 for a in l[:-1])内部通过rational_interpolate多项式有理插值见 sympy/polys/polyfuncs.py与相邻项比值变换逐层推进。六、完整工作流从猜到证再到求将上述工具串联起来可形成一条经典的具体数学研究闭环——例如从斐波那契数列出发from sympy import fibonacci, simplify, summation from sympy.concrete.guess import guess_generating_function, find_simple_recurrence fib [fibonacci(k) for k in range(15)] # 第 1 步识别递推关系 print(find_simple_recurrence(fib)) # -a(n) - a(n1) a(n2) # 第 2 步猜生成函数并解析验证 print(guess_generating_function(fib, types[ogf])) # 第 3 步对超几何项使用 Gosper 求和得到封闭形式 from sympy.concrete.gosper import gosper_sum from sympy.abc import k, n print(gosper_sum(k**3, (k, 0, n))) # n**2*(n 1)**2/4先由find_simple_recurrence得到递推结构再由guess_generating_function获得生成函数解析式最后对超几何项用gosper_sum/summation求封闭和——这正是concrete模块设计上以超几何项为中心统一支撑递推、求和与序列识别的体现。七、常用资源与扩展阅读模块文档doc/src/modules/concrete.rst求和与求积实现sympy/concrete/summations.py、sympy/concrete/products.py、sympy/concrete/expr_with_intlimits.pyGosper 算法实现sympy/concrete/gosper.py配套测试 sympy/concrete/tests/test_gosper.py序列猜想实现sympy/concrete/guess.py配套测试 sympy/concrete/tests/test_guess.py超几何判定与化简is_hypergeometricsympy/core/basic.py、hypersimpsympy/simplify/simplify.py相关算法参考文献W. Koepf《Algorithms for m-fold Hypergeometric Summation》J. Symbolic Computation, 1995为hypersimp的算法依据Graham、Knuth、Patashnik《Concrete Mathematics》与 OEIS 生成函数词条为guess_generating_function的参考来源。总结concrete模块以超几何项为统一视角将离散求和、递推求解与序列识别整合为一套可操作的符号计算工具链。掌握is_hypergeometric/hypersimp的判定与化简、gosper_sum的封闭求和、summation/product的符号求值以及guess系列工具的逆向猜想能力即可在日常研究与工程实践中完成从数列观察到封闭公式的完整闭环。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考
返回列表