泰勒（麦克劳林）级数的高效生成答案

【问题标题】：Efficient generation of Taylor (Maclaurin) series泰勒（麦克劳林）级数的高效生成
【发布时间】：2014-07-12 09:06:26
【问题描述】：

考虑函数
y=1/((1-x^5)(1-x^7)(1-x^11))

WolframAlpha 在几秒钟内计算出 MacLaurin 级数展开式的前 1000 个元素：
https://www.wolframalpha.com/input/?i=maclaurin+series+1%2F%28%281-x%5E5%29%281-x%5E7%29%281-x%5E11%29%29

出于好奇，我编写了一个非常幼稚的 java 程序，使用 BigInteger 来处理多项式系数。在伪代码中它会是这样的：

BigInt next=1;
BigInt factorial=1;
while(true){
   function=function.differentiate();
   factorial*=++next;
   print("Next coefficient is: " + function(0)/factorial);
}

该程序在计算前七个左右的系数后因 java.lang.outofmemory 异常而崩溃，因为分数的分子和分母变成了非常长的多项式。假设我的代码效率低下，但看起来 Wolfram 并没有使用他们在第一年微积分课上向您展示的相同技术。
问题是：Wolfram 使用什么？

作为比较，Wolfram 仅计算同一函数的十次导数所花费的时间比获得多项式的前 1000 项所需的时间要多得多，如果简单地完成，则需要对函数进行 1000 次微分。
https://www.wolframalpha.com/input/?i=tenth+derivative+1%2F%28%281-x%5E5%29%281-x%5E7%29%281-x%5E11%29%29

【问题讨论】：

您能否对伪代码这一步提供更详细的说明：function=function.differentiate();您是否使用 lambda 生成导数函数？
不，没有 lambda。我使用了一个非常便宜的实现：函数被定义为两个多项式：分子+分母；微分运算符使用为多项式定义的加法和乘法的quitient 规则。多项式定义为两个数组列表：coeffecients+exponents
生成 MacLaurin 系列只需要您评估函数并且它的导数为 0。这可以通过微分函数并插入值 0 来完成。或者您可以跳过所有这些，只需 compute the answers numerically.

标签： algorithm calculus taylor-series

【解决方案1】：

tl;dr：x^N 的系数是仅使用 5、7 和 11 可以划分 N 的方式数。

我不确定 Wolfram 是如何做到的，但对于这个函数，可以更有效地计算系数（使用微积分第一年结束时会看到的技术）。作为幂级数，1/(1-x)=∑_k=0^∞ x^k。但是我们可以将 x 替换为 xⁿ，这个关系仍然成立。这意味着
1/((1-x⁵)(1-x⁷)(1-x¹¹)) = (∑_{k =0}^∞x^5k)(∑_k=0^∞x^{7k)(∑_k=0^∞x^11k)}

乘以这个会很痛苦。但是所有的系数都是 1，所以我们只需要看看指数，它们加在一起。例如，Wolfram 显示 x⁴⁰ 的系数为 4，即来自 (x^5·1)(x^7·5) (x^0·11)+(x^5·0)(x^7·1)(x^{11·3 sup>)+(x^5·3)(x^7·2)(x^11·1)+(x^{5 ·8})(x^7·0)(x^11·0).}

但是如果我们只需要添加指数，那么我们就不需要关心系数或变量 x。最后，x^N的系数就是N可以写成5s、7s、11s之和的方式数。这是分区问题的受限版本，但同样的想法仍然成立。特别是，动态规划方法将能够计算线性时间和空间中的系数。

【讨论】：

【解决方案2】：

不确定分数的分子，但我知道为什么它的分母增长太快了：

factorial*=factorial+1;

不是你计算阶乘的方式。每次迭代平方分母上的“阶乘”值！所以你会得到 1, 2, 6, 42, 1806, 3263442... 相比之下，阶乘是 1, 2, 6, 24, 120, 720...

要增量计算阶乘，请维护一个循环计数器，并每次将factorial 乘以那个。

【讨论】：

对不起，这是我的问题中的一个错字。在代码中，我正确计算了阶乘。

【解决方案3】：

有理函数（以及这个特定的函数）既不需要微分也不需要阶乘。计算序列的一种方法是将每个因子扩展为自己的序列（例如1/(1 - x^5) = sum(n=[0,inf] x^(5n))，然后将结果乘以多项式。

【讨论】：

【解决方案4】：

您可以对正式的幂级数进行每个操作。给定 f,g 的幂级数，您可以找到 f(z)+g(z)、f(z)g(z)、f(z)/g(z)、f(g( z))，甚至 f^-1(z)。使用这些方法，您可以在多项式时间内计算几乎任何函数的幂级数。

在特殊情况下有更有效的方法。如果 f(z) 具有幂级数，则 f(z)/(1 - z) 的幂级数的系数只是 f(z) 的幂级数的部分和。因此，如果 f_n 是 f 的级数，则 g(z) = f(z)/(1 - z) 的级数由 g_n = f_n + g_(n-1) 给出。

您可以将此扩展到除以任何多项式。该算法与多项式的长除法基本相同。例如，让我们计算 1/(1 - z^2)。我们在分子上加上和减去 z^2 得到 (1 - z^2 + z^2)/(1 - z^2) = 1 + z^2/(1 - z^2)。然后我们加减 z^4 得到 (z^2 - z^4 + z^4)/(1 - z^2) = z^2 + z^4/(1 - z^2)。继续这样你会发现 1/(1 - z^2) = 1 + z^2 + z^4 + z^6 等等。

当您对 n 次的一般多项式执行此操作时，您总是有一个少于 n 项的分子。您可以将这些项的系数存储在一个数组中，并将其用作您的状态。从一个状态，您可以计算幂级数中的下一项和时间 O(n) 中的下一个状态。这为您提供了一个 O(nk) 时间算法来查找 1/p(z) 幂级数中的前 k 个项。

请注意，在 z=z0 点计算幂级数与在 z=z0 处求所有导数相同，因此这两个问题是等价的。您可以计算符号变量点处的幂级数以找到导数公式，因此理论上没有理由说明 Wolfram 在求 n 次导数时要慢得多。

【讨论】：