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

Does Python have a function which computes multinomial coefficients?

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

我在任何标准库中都找不到任何此类函数。 对于二项式系数(其中多项式系数是一个推广)有 scipy.special.binom and also scipy.misc.comb. Also, numpy.random.multinomial draws samples from a multinomial distribution, and sympy.ntheory.multinomial.multinomial_coefficients returns 一个与多项式系数相关的字典。

但是,我无法找到一个正确的多项式系数函数,给定 a,b,...,z returns (a+b+...+z)!/(a!b!...z!). 我错过了吗?有 none 可用的充分理由吗?

我很乐意为 SciPy say 贡献一个有效的实现。 (我必须弄清楚如何做出贡献,因为我从来没有这样做过)。

作为背景,他们在扩展(a+b+...+z)^n时确实出现了。另外,他们计算了存款方式 a+b+...+z 不同的对象放入不同的容器中,这样第一个容器包含 a 对象等。我偶尔需要它们来解决 Project Euler 问题。

顺便说一句,其他语言也提供此功能:Mathematica, MATLAB, Maple.

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

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

你写了 "sympy.ntheory.multinomial.multinomial_coefficients returns 一本与多项式系数相关的字典",但如果你知道如何提取,从评论中不清楚该字典中的特定系数。使用维基百科中的符号 link,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

Busy Beaver 给出了 scipy.special.binom 的答案。该实现的一个潜在问题是 binom(n, k) returns 是一个浮点值。如果系数足够大,它就不会精确,因此它可能无法帮助您解决 Project Euler 问题。您可以使用带有参数 exact=Truescipy.special.comb 而不是 binom。这是 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

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

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

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

使用阶乘

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

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

开始 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

你自己的答案(采纳的)很好,而且特别简单。但是,它确实有一个显着的低效率问题:您的外层循环 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