有理函数级数展开的最佳算法答案

【问题标题】：Best algorithm for series expansion of Rational function有理函数级数展开的最佳算法
【发布时间】：2014-05-30 07:52:47
【问题描述】：

我需要用 C++ 编写函数，它可以有效地找到给定有理函数 (P(x) / Q(x)) 的泰勒级数的系数。

函数参数将是多项式的幂（在分母和分母上相等），具有多项式系数和展开项数的两个数组。

我的想法是这样的。考虑身份

P(x) / Q(x) = R(x) + ...

其中R(x) 是一个多项式，其项数等于我需要找到的系数数。然后我可以将两边乘以Q(x) 并得到

P(x) = R(x) * Q(x)

R(x) * Q(x) - P(x) = 0

因此，所有系数都应为零。这是具有 O(n^3) 算法求解的方程组。 O(n^3) 没有我想要的那么快。

有没有更快的算法？

我知道级数的系数满足线性递推关系。这让我觉得O(n)算法是可能的。

【问题讨论】：

@DavidEisenstat 谢谢，这简化了任务。但是我怎样才能找到带有长除法的1 / Q(x)？我以为除法只能找到商和余数。
@DavidEisenstat 谢谢！避免分离也会加快速度，因为乘法将是 O(nm)*，其中 n - 项数，m - 度。顺便说一句，我不明白为什么我的答案被标记为过于宽泛，它需要一种特定的算法。
@DavidEisenstat 你去吧 :) 也许我们应该开始一个关于元 [algorithm] 状态的话题。关于是否将这些问题转移到 CS 上，肯定需要做出一些一般性决策
@NiklasB。随意。我可能会参与，但 (i) 我认为 MSO 严重受损，并且 (ii) 我不喜欢那些 SO rep/MSO 代表比率为个位数（或更糟，小于 1）的人的回答。
@DavidEisenstat 有了它，你的第二个参数不再有效:)

标签： c++ algorithm taylor-series

【解决方案1】：

我将要描述的算法在数学上由formal power series 证明是合理的。每个具有泰勒级数的函数都有一个正式的幂级数。反之则不成立，但如果我们对泰勒级数的函数进行算术运算并得到一个泰勒级数的函数，那么我们可以对形式幂级数进行相同的算术运算并得到相同的答案。

形式幂级数的长除法算法类似于您可能在学校学过的long division 算法。我将在示例(1 + 2 x)/(1 - x - x^2) 上进行演示，其系数等于Lucas numbers。

分母必须有一个非零常数项。我们先写分子，也就是第一个残差。

             --------
1 - x - x^2 ) 1 + 2 x

[ 我们将残差的最低阶项 (1) 除以分母的常数项 (1)，然后将商放在首位。

              1
             --------
1 - x - x^2 ) 1 + 2 x

现在我们将1 - x - x^2 乘以1 并从当前残差中减去它。

              1
             --------
1 - x - x^2 ) 1 + 2 x
              1 -   x - x^2
              -------------
                  3 x + x^2

再做一次。

              1 + 3 x
             --------
1 - x - x^2 ) 1 + 2 x
              1 -   x -   x^2
              ---------------
                  3 x +   x^2
                  3 x - 3 x^2 - 3 x^3
                  -------------------
                        4 x^2 + 3 x^3

又一次。

              1 + 3 x + 4 x^2
             ----------------
1 - x - x^2 ) 1 + 2 x
              1 -   x -   x^2
              ---------------
                  3 x +   x^2
                  3 x - 3 x^2 - 3 x^3
                  -------------------
                        4 x^2 + 3 x^3
                        4 x^2 - 4 x^3 - 4 x^4
                        ---------------------
                                7 x^3 + 4 x^4

又一次。

              1 + 3 x + 4 x^2 + 7 x^3
             ------------------------
1 - x - x^2 ) 1 + 2 x
              1 -   x -   x^2
              ---------------
                  3 x +   x^2
                  3 x - 3 x^2 - 3 x^3
                  -------------------
                        4 x^2 + 3 x^3
                        4 x^2 - 4 x^3 - 4 x^4
                        ---------------------
                                7 x^3 + 4 x^4
                                7 x^3 - 7 x^4 - 7 x^4
                                ---------------------
                                       11 x^4 + 7 x^5

单个除法有点无聊，因为我使用了一个以1 开头的除数，但如果我使用了2 - 2 x - 2 x^2，那么商中的所有系数都将除以2。

【讨论】：

谢谢，这很有帮助！我想找到的所有展开都是整数系数，这可能保证前导系数为 1。

【解决方案2】：

这可以在O(n log n) 时间内完成，任意P 和Q 度为n。更准确地说，这可以在M(n) 中完成，其中M(n) 是多项式乘法的复杂度，它本身可以在O(n log n) 中完成。

首先，级数展开式的第一个n 项可以简单地视为次数为n-1 的多项式。

假设您对P(x)/Q(x) 的级数扩展的第一个n 术语感兴趣。存在一种算法，可以在上面定义的 M(n) 时间内计算 Q 的倒数。

T(x) 的逆 Q(x) 满足 T(x) * Q(x) = 1 + O(x^N)。 IE。 T(x) * Q(x) 正好是 1 加上一些误差项，其系数都在我们感兴趣的第一个 n 项之后，所以我们可以删除它们。

现在P(x) / Q(x) 只是P(x) * T(x)，这只是另一个多项式乘法。

您可以在我的开源库 Altruct 中找到计算上述逆的实现。请参阅series.h 文件。假设您已经有一个计算两个多项式乘积的方法，那么计算逆的代码大约有 10 行（分治法的一种变体）。

实际算法如下：假设Q(x) = 1 + a1*x + a2*x^2 + ...。如果a0 不是1，您可以简单地将Q(x) 和之后的倒数T(x) 与a0 相除。假设在每一步你都有L 的逆项，所以Q(x) * T_L(x) = 1 + x^L * E_L(x) 出现一些错误E_L(x)。最初T_1(X) = 1。如果你在上面插入这个，你会得到一些E_1(x) 的Q(x) * T_1(x) = Q(x) = 1 + x^1 * E_1(x)，这意味着这对L=1 成立。现在让我们在每一步加倍L。您可以从上一步获得E_L(x) 作为E_L(x) = (Q(x) * T_L(x) - 1) / x^L，或者在实现方面，只需删除产品的第一个L 系数。然后，您可以将上一步中的T_2L(x) 计算为T_2L(x) = T_L(x) - x^L * E_L(x) * T_L(x)。错误将是E_2L(x) = - E_L(x)^2。现在让我们检查归纳步骤是否成立。

Q(x) * T_2L(x)
= Q(x) * (T_L(x) - x^L * E_L(x) * T_L(x))
= Q(x) * T_L(x) * (1 - x^L * E_L(x))
= (1 + x^L * E_L(x)) * (1 - x^L * E_L(x))
= 1^2 - (x^L * E_L(x))^2
= 1 + x^2L * E_2L(x)

Q.E.D.

我很确定多项式除法的计算效率不可能高于乘法，如下表所示，该算法仅比单次乘法慢 3 倍：

   n      mul        inv      factor
10^4       24 ms      80 ms    3,33x
10^5      318 ms     950 ms    2,99x
10^6    4.162 ms  12.258 ms    2,95x
10^7  101.119 ms 294.894 ms    2,92x

【讨论】：

【解决方案3】：

如果你仔细观察你的计划得到的系统，你会发现它已经是对角线的，不需要 O(n^3) 来解决。它简单地退化为线性递归（P[]、Q[] 和 R[] 是相应多项式的系数）：

R[0] = P[0]/Q[0]
R[n] = (P[n] - sum{0..n-1}(R[i] * Q[n-i]))/Q[0]

由于 Q 是多项式，因此总和不超过 deg(Q) 项（因此需要恒定时间来计算），从而使整体复杂度渐近线性。您还可以查看递归的矩阵表示以获得（可能）更好的渐近性。

【讨论】：

对于一个 R(i) 值是线性的，但是对于所有值 R(0), ..., R(n) 它变成 O(n^2)
没有。如果 Q 有无限个系数，它将变为 O(n^2)。总和没有比 deg(Q) 更多的项，因此需要恒定时间来计算。我将编辑我的答案 - 目前还不清楚。谢谢。
然后变成 O(n*q)，其中 q - Q 的度数。看起来和接受的帖子一样。