Python有计算多项式系数的函数吗?

27

我正在寻找一个Python库函数,用于计算多项式系数

我在任何标准库中都没有找到这样的函数。对于二项式系数(多项式系数的一般化),有 scipy.special.binomscipy.misc.comb。 此外,numpy.random.multinomial 从多项式分布中绘制样本,并且sympy.ntheory.multinomial.multinomial_coefficients返回与多项式系数相关的字典。

但是,我找不到适当的多项式系数函数,它可以给定 a,b,...,z 返回 (a+b+...+z)!/(a! b! ... z!)。我错过了吗?没有提供这样的函数有好的理由吗?

如果需要,我很乐意向SciPy捐献一个有效的实现(我必须弄清楚如何贡献,因为我从未这样做过)。

背景:当展开 (a+b+...+z)^n 时,它们会出现。此外,它们计算将 a+b+...+z 个不同对象放入不同容器中的方式,使得第一个容器包含 a 个对象等等。 我偶尔需要它们用于Project Euler问题。

顺便说一句,其他语言确实提供了这个函数:MathematicaMATLABMaple


1
由于这是我的第一个问题,我很好奇为什么这个问题会被踩。我的搜索没有得到我想要的答案。此外,我已经提供了一些背景信息。提前感谢任何澄清。 - Reiner Martin
2
请帮我理解,这样一个具体的技术问题会如何吸引充满个人观点的答案?要么该函数是可用的,但可能被隐藏得很好或使用了不寻常的名称,要么库设计者选择不实现它有一个很好的理由,或者它只是一个空缺(我很乐意填补)。请注意,我不是在寻求建议。 - Reiner Martin
1
这正是我所做的,这就是为什么我感到困惑的原因。 - Reiner Martin
@Reiner,你能举个这样的问题的例子吗? - Nick stands with Ukraine
是的,它出现在 https://projecteuler.net/problem=369,这正是我正在解决的问题(我的第401个)。 - Reiner Martin
显示剩余13条评论
6个回答

13
为了部分回答我的问题,这是我对多项式函数的简单且相当有效的实现:

def multinomial(lst):
    res, i = 1, 1
    for a in lst:
        for j in range(1,a+1):
            res *= i
            res //= j
            i += 1
    return res

从迄今为止的评论中看来,标准库中没有任何高效实现该函数的方法。

更新(2020年1月)。正如唐·哈奇(Don Hatch)在评论中指出的那样,这可以通过寻找最大参数来进一步改进(特别是对于主导所有其他参数的情况):

def multinomial(lst):
    res, i = 1, sum(lst)
    i0 = lst.index(max(lst))
    for a in lst[:i0] + lst[i0+1:]:
        for j in range(1,a+1):
            res *= i
            res //= j
            i -= 1
    return res

根据我的经验,使用sum()math.factorial()计算公式会更快。 - jarno
3
当然,这取决于你想要做什么。用那种方法,你无法解决任何(3位数)的欧拉计划问题。 - Reiner Martin
嗨Reiner,我在意识到这个答案其实不太好之前不小心点赞了它。 multinomial([10**7,2])需要0.8秒。multinomial([10**8,2])需要8秒。multinomial([10**9,2])会使我的机器崩溃(对于python2)或需要80秒(对于python3)。所有这些都可以瞬间完成。 - Don Hatch
你可以很容易地修复它。一种方法是:用res = 1i = sum(lst)index_of_max = lst.index(max(lst))替换函数中的前两行,然后用for a in lst[:index_of_max] + lst[index_of_max+1:]:替换i+=1,最后将i-=1替换掉。 - Don Hatch
事实上,微不足道的 multinomial([10**8,1]) 甚至更微不足道的 multinomial([10**8,0]) 都展示了您代码中的问题。 - Don Hatch
如果执行时间至关重要,您可以使用存储在内存中的“薄”二项式系数表;请参见如何高效计算多项式? \math.stackexchange - CopyPasteIt

10
没有Python内置的多项式库或函数。
无论如何,这次数学可能会帮助你。事实上,有一种简单的方法可以计算多项式。

保持关注性能是通过使用多项式系数的特征将其重写为二项式系数的乘积来实现的。

当然,这里

感谢 scipy.special.binom 和递归的神奇之处,您可以像这样解决问题:
from scipy.special import binom

def multinomial(params):
    if len(params) == 1:
        return 1
    return binom(sum(params), params[-1]) * multinomial(params[:-1])

其中params = [n1, n2, ..., nk]

注意:将多项式拆分为二项式的乘积也有助于通常情况下防止溢出。


谢谢,这基本上是我已经做的(将递归展开成循环以提高性能)。然而,这不是我的问题;相反,如果它在 SciPy 中不存在(似乎确实不存在),是否将其包含其中会是一个好主意,因为它已经存在于其他几种语言中?相信我,我不需要数学讲座。 - Reiner Martin
这真的是一个需要问“scipy”维护者的问题。 - chepner
好的,很好的建议,会这样做。可能这不是正确的听众。谢谢。 - Reiner Martin
话虽如此,我不确定除了 SciPy 之外是否还有其他包含它的库。 - Reiner Martin
1
一年多过去了,我来更新一下。对于我的问题解决编程(比如Project Euler),我已经完全转向Julia,我认为它更适合这样的工作。顺便说一句,Julia的组合数学库中包含了多项式函数,这是你所期望的。 - Reiner Martin
显示剩余4条评论

2
你写道"sympy.ntheory.multinomial.multinomial_coefficients返回与多项式系数相关的字典",但从该评论中并不清楚你是否知道如何从该字典中提取特定的系数。使用维基百科链接中的符号,SymPy函数为给定的mn提供了所有的多项式系数。如果你只想要一个特定的系数,只需从字典中提取即可:
In [39]: from sympy import ntheory

In [40]: def sympy_multinomial(params):
    ...:     m = len(params)
    ...:     n = sum(params)
    ...:     return ntheory.multinomial_coefficients(m, n)[tuple(params)]
    ...: 

In [41]: sympy_multinomial([1, 2, 3])
Out[41]: 60

In [42]: sympy_multinomial([10, 20, 30])
Out[42]: 3553261127084984957001360

Busy Beaver提供了一个用scipy.special.binom编写的答案。这种实现的潜在问题在于binom(n, k)返回一个浮点数值。如果系数足够大,它将不是精确值,因此它可能无法帮助您解决Project Euler问题。您可以使用scipy.special.comb代替binom,并使用参数exact=True。这是Busy Beaver的函数,修改后使用comb

In [46]: from scipy.special import comb

In [47]: def scipy_multinomial(params):
    ...:     if len(params) == 1:
    ...:         return 1
    ...:     coeff = (comb(sum(params), params[-1], exact=True) *
    ...:              scipy_multinomial(params[:-1]))
    ...:     return coeff
    ...: 

In [48]: scipy_multinomial([1, 2, 3])
Out[48]: 60

In [49]: scipy_multinomial([10, 20, 30])
Out[49]: 3553261127084984957001360

谢谢,我确实没有研究如何从那个字典中获取特定系数。对我来说,生成完整的字典只为计算单个值似乎非常浪费,因此我对一些较大的输入进行了快速计时。例如,sympy_multinomial([123,134,145])所需的时间是我自己简单函数的1000多倍。 - Reiner Martin

2
这里有两种方法,一种使用阶乘,另一种使用斯特林公式

使用阶乘

您可以定义一个函数,使用矢量化代码(而不是for循环)在单行中返回多项式系数,如下所示:

from scipy.special import factorial

def multinomial_coeff(c):
    return factorial(c.sum()) / factorial(c).prod()

(其中c是包含每个不同对象计数的np.ndarray)。使用示例:
>>> import numpy as np
>>> coeffs = np.array([2, 3, 4])
>>> multinomial_coeff(coeffs)
1260.0

在某些情况下,这可能会更慢,因为您将多次计算某些阶乘表达式,但在其他情况下,这可能会更快,因为我相信numpy自然地并行化矢量化代码。此外,这减少了程序中所需的行数,并且可能更易于阅读。如果有人有时间对这些不同选项进行速度测试,那么我很想看到结果。
使用斯特林近似
事实上,多项式系数的对数根据斯特林近似计算更快,而且允许计算更大的系数。
from scipy.special import gammaln

def log_multinomial_coeff(c):
    return gammaln(c.sum()+1) - gammaln(c+1).sum()

使用例子:
>>> import numpy as np
>>> coeffs = np.array([2, 3, 4])
>>> np.exp(log_multinomial_coeff(coeffs))
1259.999999999999

3
实际上,基于 Stirling 近似,多项式系数的对数计算速度更快,并且可以计算更大的系数:from scipy.special import gammaln def log_multinomial_coeff(c): return gammaln(c.sum()+1) - gammaln(c+1).sum() - Jake Levi
我已经更新了答案,包括之前的评论。 - Jake Levi

1
你自己的答案(被接受的)相当不错,而且特别简单。但是,它确实存在一个重要的低效性:你的外部循环for a in lst比必要执行了一次。在第一次通过该循环时,ij的值总是相同的,因此乘法和除法没有任何作用。在你的例子multinomial([123, 134, 145])中,有123个不需要的乘法和除法,增加了代码的时间。
我建议找到参数中的最大值并将其删除,以便不进行这些不必要的操作。这增加了代码的复杂性,但减少了执行时间,特别是对于大数列表的短列表。我的代码在111微秒内执行multcoeff(123, 134, 145),而你的代码需要141微秒。虽然这不是一个很大的增加,但也可能很重要。所以这是我的代码。这还将单独的值作为参数,而不是列表,所以这是与你的代码的另一个区别。
def multcoeff(*args):
    """Return the multinomial coefficient
    (n1 + n2 + ...)! / n1! / n2! / ..."""
    if not args:  # no parameters
        return 1
    # Find and store the index of the largest parameter so we can skip
    #   it (for efficiency)
    skipndx = args.index(max(args))
    newargs = args[:skipndx] + args[skipndx + 1:]

    result = 1
    num = args[skipndx] + 1  # a factor in the numerator
    for n in newargs:
        for den in range(1, n + 1):  # a factor in the denominator
            result = result * num // den
            num += 1
    return result

1

Python 3.8 开始,

我们可以在不使用外部库的情况下实现它:

import math

def multinomial(*params):
  return math.prod(math.comb(sum(params[:i]), x) for i, x in enumerate(params, 1))

multinomial(10, 20, 30) # 3553261127084984957001360

1
如果性能是一个问题,这种方法会效率有所下降,因为很多项都会抵消掉,而这种方法没有利用它们。 - Reiner Martin

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接