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

Question

我正在寻找一个计算 multinomial coefficients 的 Python 库函数。

我在任何标准库中都找不到任何这样的函数。 对于二项式系数（其中多项式系数是一种概括），有scipy.special.binom 和scipy.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 个不同的对象到不同的 bin 中，这样第一个 bin 包含 a 个对象等。我偶尔需要它们来解决 Project Euler 问题。

顺便说一句，其他语言确实提供此功能：Mathematica、MATLAB、Maple。

Answer 1

为了部分回答我自己的问题，这是我对多项式函数的简单且相当有效的实现：

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

到目前为止，从 cmets 看来，任何标准库中都不存在该函数的有效实现。

更新（2020 年 1 月）。 正如 Don Hatch 在 cmets 中指出的那样，这可以通过寻找最大的参数来进一步改进（尤其是在它支配所有其他参数的情况下）：

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

Answer 2

不，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].

注意：将多项式拆分为二项式的乘积通常也可以很好地防止溢出。

Answer 3

您写了 “sympy.ntheory.multinomial.multinomial_coefficients 返回与多项式系数相关的字典”，但从该评论中不清楚您是否知道如何从该字典中提取特定系数。使用维基百科链接中的符号，SymPy 函数为您提供 all 给定 m 和 n 的多项式系数。如果您只想要一个特定的系数，只需将其从字典中拉出即可：

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

忙碌的海狸给出了一个用scipy.special.binom 写的答案。该实现的一个潜在问题是binom(n, k) 返回一个浮点值。如果系数足够大，它将不准确，因此它可能无法帮助您解决 Project Euler 问题。除了binom，您可以使用scipy.special.comb，并带有参数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

Answer 4

这里有两种方法，一种使用阶乘，一种使用Stirling's approximation。

使用阶乘

您可以使用矢量化代码（而不是 for-loops）定义一个函数以在一行中返回多项式系数，如下所示：

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 自然地并行化矢量化代码。这也减少了程序中所需的行数，并且可以说更具可读性。如果有人有时间对这些不同的选项进行速度测试，那么我很想看看结果。

使用Stirling's approximation

事实上，多项式系数的对数计算起来更快（基于Stirling's approximation）并且允许计算更大的系数：

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

Answer 5

您自己的答案（被接受的）非常好，而且特别简单。但是，它确实有一个显着的低效率：您的外部循环 for a in lst 执行的次数超出了必要的次数。在第一次通过该循环时，i 和 j 的值始终相同，因此乘法和除法没有任何作用。在您的示例 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

Answer 6

开始Python 3.8，

• 因为标准库现在包含 math.comb 函数（二项式系数）
• 并且由于多项式系数可以计算为product of binomial coefficients

我们可以在没有外部库的情况下实现它：

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