我需要用 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)算法是可能的。
我将要描述的算法在数学上由形式幂级数证明是合理的。每个具有泰勒级数的函数都有一个正式的幂级数。反之则不成立，但如果我们对泰勒级数的函数进行算术运算并得到一个泰勒级数的函数，那么我们可以对形式幂级数进行相同的算术运算并得到相同的答案。
形式幂级数的长除法算法就像你在学校学过的长除法算法一样。我将在示例中演示它(1 + 2 x)/(1 - x - x^2)，它的系数等于卢卡斯数。
分母必须有一个非零常数项。我们从写分子开始，它是第一个残差。
--------
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。
这可以在O(n log n)任意时间P和Q度数上完成n。更准确地说，这可以在 中完成M(n)，其中M(n)多项式乘法的复杂性本身可以在 中完成O(n log n)。
首先，n级数展开式的第一项可以简单地看成是度的多项式n-1。
假设您对n的级数展开的第一项感兴趣P(x)/Q(x)。存在一种算法，可以计算上面定义Q的时间倒数。M(n)
的倒数T(x)满足。即正好加上一些误差项，其系数都在我们感兴趣的第一个项之后，所以我们可以去掉它们。Q(x)T(x) * Q(x) = 1 + O(x^N)T(x) * Q(x)1n
现在P(x) / Q(x)很简单P(x) * T(x)，这只是另一个多项式乘法。
您可以在我的开源库Altuct中找到计算上述逆的实现。请参阅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. 如果你在上面插入这个，你会得到Q(x) * T_1(x) = Q(x) = 1 + x^1 * E_1(x)一些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)
量子点
我很确定多项式除法的计算效率不可能比乘法更有效，如下表所示，该算法仅比单次乘法慢 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
如果你仔细观察你的计划得到的系统，你会发现它已经是对角线的，不需要 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)项（因此需要恒定时间来计算），从而使整体复杂度渐近线性。您还可以查看递归的矩阵表示以获得（可能）更好的渐近。