Part 0 泰勒展开的推广&多项式牛顿迭代
§ 0.1 记号和约定
为了不在下文引起混淆，这里简述一下这篇文章使用的记号：
- \(f(x)\) 或者 \(F(x)\) ：一个形式幂级数，此处的 \(x\) 是一个形式记号。
- \(F(f(x))\) ：两个多项式的复合。
- \(\mathcal{F}[f(x)]\) 或者：一个自变量为多项式的函数。此处的 \(\mathcal{F}\) 的自变量是多项式，而不是实数。而 \(x\) 仍然是一个表示多项式系数位置的形式符号，没有实际意义。为了体现出它与复合函数的区别，这里使用方括号和手写体。
- \(f'(x)\) ：多项式 \(f(x)\) 的形式导数，换言之可以看做对于多项式进行的一种类似加减乘除的运算。
- \(\int f(x)\) ：多项式 \(f(x)\) 的形式积分，含义同上。
- \(f(x^m)\) ：将多项式每一系数的下标翻 \(m\) 倍后得到的新多项式，并不是复合函数。
接下来我们会讨论如何求出 \(\mathcal{F}[f(x)]\) 的导数，并尝试将其泰勒展开。
§ 0.2 混账，你导了甚么！
第一个遇到的问题就是：我们该如何定义 \(R[[x]]\) 上的导函数。
要定义导函数，就需要定义微分，同时微分建立在无穷小量上，而要定义无穷小量，就需要先定义极限。
可惜的是，我等凡人实在是难以想象一个 \(R[[x]]\) 上的极限该怎么定义。所幸，我们可以耍耍小聪明，绕过这一个障碍。
观察一个以多项式为自变量的函数 \(\mathcal F[f(x)]\) ，在OI中，它的解析式一般也是人畜无害的多项式形式，长成这样：
\[\mathcal F[f(x)] = k_0 +k_1 \mathcal f(x)+k_2 \mathcal f^2(x) +\cdots +k_n \mathcal f^n(x) \]因为这里不存在针对 \(f(x)\) 的代入求值，多项式 \(f(x)\) 的存在更加接近于一个 “常量” 而非 “函数”，又考虑到多项式们有和实数几乎完全一致的加减乘运算，我们可以直接把一个多项式 \(f(x)\) 视作一个实数：
\[\mathcal F[\mathcal X] = k_0 +k_1 \mathcal X+k_2 \mathcal X^2 +\cdots +k_n \mathcal X^n \]现在就是我们的小聪明了：直接照抄实数域上的导函数定义，捏一个 \(\mathcal F\) 的 “导函数” 出来！
\[\mathcal F'[\mathcal X] = k_1 +2k_2 \mathcal X + 3k_3 \mathcal X^2 +\cdots +nk_n \mathcal X^{n-1} \]我们暂且把这个东西叫做 “形式导数”，尽管直接跳过了极限和微分，我们仍然可以由它开始一步步打通泰勒展开的定义。
§ 0.3 噫！好！我展了！
回顾高等数学，实数域中泰勒展开的由来无非是挑选一个点，强制令新函数在这一点上的各阶导数和原来的函数相等。而且泰勒展开还有一个特点：如果被展开的函数是一个多项式函数或者能用幂级数表达的函数，其正确性的证明可以完全不依靠微分和无穷小量来得出。
如果有人想看的话我就去补写一下证明qwq
现在实在是懒得搞了
那么，既然多项式环有着和实数几乎一致的加减乘（导数需要的三种运算）和数乘运算（泰勒展开需要的额外运算），那么就算把实数版本的证明中的 \(x\) 改为多项式而不是实数，证明也依然成立。因此没有理由相信泰勒展开在多项式环上不能用。
换句话说，只要有了导数的形式，对多项式函数和指对函数泰勒展开就是正确的。
以自然底数为底的指数、对数函数可以写成无穷项幂级数的形式，所以它们也可以包含在内。
如果被展开函数是无穷多项的多项式函数，那么的确会用到极限。
但是因为求极限的对象是无穷多个多项式的加和，取到极限的变量仍是一个实数，我们仍然不需要定义 \(R[[x]]\) 上的极限。
正好，我们的 \(\mathcal F\) 也有着多项式的形式，这就为我们创造了机会，可以再次跳过微分，直接捏出一个新的泰勒展开：
\[\mathcal F[\mathcal X] = \sum_{k=0}^{+\infty} \frac{ \mathcal F^{(k)} [\mathcal X_0] }{k!}(\mathcal X - \mathcal X_0)^k \]接下来可以把 \(\mathcal X\) 换成我们更熟悉的多项式 \(f(x)\) .
\[\mathcal F[ f(x) ] = \sum_{k=0}^{+\infty} \frac{ \mathcal F^{(k)} [ f_0(x) ] }{k!}( f(x) - f_0(x) )^k \]好了！我们做完了！
虽然可能不甚严谨，但是管他呢，用就完了。
这里有一篇博客提供了另一种理解：多项式牛顿迭代([email protected])
简单来说，我们挑选所有能让 \(f_0(x)\) 收敛的实数 \(x\) ，把它代入 \(f_0(x)\) 可以求出一些值，于是就可以把 \(f_0(x)\) 视作一个实数变量，同时自然把 \(\mathcal F[f(x)]\) 变成了一个以实数变量 \(f_0(x)\) 为自变量的函数，对其强行进行泰勒展开一定是正确的。
接下来再对所有的能让 \(f(x)\) 收敛的实数 \(x\) 代入其中，由于 \(f(x)\) 和 \(f_0(x)\) 的值域都为 \(R\) ，多项式表达式的差异可以大致等同于函数值（\(\mathcal F\) 函数的自变量）的差异。
综上，我们就得到了这种泰勒展开的正确性。
§ 0.4 开始迭代
多项式牛顿迭代可以解决这样的问题：
已知多项式 \(g(x)\) ，求满足这种条件的多项式 \(f(x)\)，系数对 \(998244353\) 取模。
\[g(f(x)) \equiv 0 \pmod{x^n} \]换句话说就是求零点。
这里多项式取模的意思就是只保留次数最低的 \(n\) 位，后面的东西不要了。
考虑用倍增实现。每一轮通过模 \(x^{\lceil \frac{n}{2} \rceil}\) 意义下的解 \(f_0(x)\) 得到模 \(x^n\) 意义下的解 \(f(x)\) .
第一轮中，我们需要求出 \(g(x)\) 在模意义下的根，这个单独处理。
如何迭代计算模 \(x^n\) 意义下的答案呢？我们可以用刚才推导出的泰勒展开，把 \(\mathcal G\)（把 \(g\) 视作一个自变量为多项式的函数）在 \(f_0(x)\) 处展开：
\[\mathcal G[f(x)] = \sum_{k=0}^{+\infty} \frac{ \mathcal G^{(k)}[f_0(x)] }{k!}(f(x) - f_0(x))^k \]然后，观察到 \((f(x) - f_0(x))^k\) 在 \(k>2\) 情况下前 \(n\) 位都是零，这样我们就可以舍弃掉后面的东西了
\[\begin{align*} &\mathcal G[f(x)] \equiv 0 &\pmod{x^n} \\ \Rightarrow & \sum_{k=0}^{+\infty} \frac{ \mathcal G^{(k)}[f_0(x)] }{k!}(f(x) - f_0(x))^k \equiv 0 &\pmod{x^n} \\ \Rightarrow & \mathcal G[f_0(x)] + \mathcal G'[f_0(x)] \times(f(x) - f_0(x))& \pmod{x^n} \\ \Rightarrow & \mathcal G[f_0(x)] + \mathcal G'[f_0(x)] f(x) - \mathcal G'[f_0(x)] f_0(x)& \pmod{x^n} \end{align*} \]当然还有另一种方式理解为什么要用倍增实现：
普通的牛顿迭代每次可以把精度翻倍，那么对应到多项式上，每次确定的项数也会翻倍。
化简式子
\[\mathcal G[f_0(x)] + \mathcal G'[f_0(x)] f(x) - \mathcal G'[f_0(x)] f_0(x)\equiv 0 \pmod{x^n} \]得到
\[f(x)\equiv f_0(x) - \frac{\mathcal G[f_0(x)]}{\mathcal G'[f_0(x)]}\pmod{x^n} \]诶嘿，正好是经典牛顿迭代法的形式。难道不应该是吗
这就做完了，时间复杂度分析如下：
已知多项式求逆的复杂为 \(O(n\log n)\) 求导和加减的复杂度是 \(O(n)\) ，那么就有
\[T(n) = T(\frac{n}{2}) + O(n) + O(n\log n) = O(n\log n) \]最终的复杂度是 \(O(n\log n)\) .
但是实际跑起来常数巨大无比。
Part 1 多项式乘法逆
§ 1.1 常规方法
在很多多项式算法（例如刚刚讲到的牛顿迭代）中都需要用一个多项式在模 \(x^n\) 意义下除以另一个多项式。这种操作等价于乘上另一个多项式的乘法逆元，现在我们来求这个东西。
根据刚才牛顿迭代的启发，我们用倍增解决这个问题，已知已经求得了模 \(x^{n/2}\) 意义下 \(f(x)\) 的逆元。
\[f_0^{-1}(x)f(x) \equiv 1 \pmod{x^{n/2}} \]又因为目标答案长这样：
\[f^{-1}(x)f(x) \equiv 1 \pmod{x^{n/2}} \]两式相减，得到：
\[f^{-1}(x) - f_0^{-1}(x) \equiv 0 \pmod{x^{n/2}} \]式子两边平方，模数也平方，依然能整除：
\[f^{-2}(x) - 2 f_0^{-1}(x)f^{-1}(x) + f_0^{-2}(x) \equiv 0 \pmod{x^{n}} \]然后式子把他自己收拾好了：
\[f^{-2}(x) \equiv f_0^{-1}(x) \big( 2f^{-1}(x) - f_0^{-1}(x) \big) \pmod{x^{n}} \]两边同时乘上 \(f(x)\) :
\[f^{-1}(x) \equiv f_0^{-1}(x) \big( 2 - f(x) f_0^{-1}(x) \big) \pmod{x^{n}} \]递归计算就做完了！采用本式计算乘法逆的话每轮会进行两次 NTT 和一次 INTT .
模板代码（含NTT）
namespace poly{
const LL gn = 3;
const LL gni = 332748118;
int rev[maxn<<2],now_rev;
void ReArrange(LL *F,int N)
{
if(now_rev != N)
{
rev[0] = 0; now_rev = N;
for(int i=1; i<N; i++) rev[i] = (rev[i/2]/2) + (i&1)*N/2;
}
for(int i=0; i<N; i++) if(rev[i] > i) swap(F[rev[i]], F[i]);
return;
}
void NTT(LL *F,int N,int inv)
{
ReArrange(F,N);
for(int len=1; len<N; len<<=1)
{
LL stp = qpow((inv == 1) ? gn : gni, (mod-1)/(len<<1));
for(int st=0; st<N; st+=(len<<1))
{
LL omega = 1;
for(int i=0; i<len; i++)
{
LL U = F[st+i], V = M(F[st+len+i], omega);
F[st+i] = madd(U,V); F[st+len+i] = msub(U,V);
omega = M(omega,stp);
}
}
}
if(inv == -1)
{
LL Ninv = qpow(N,mod-2);
for(int i=0; i<N; i++) F[i] = M(F[i],Ninv);
}
return;
}
LL pinv[maxn];
void Inverse(LL *f1,int n)
{
int K = 1; while(K<n) {K = K<<1;} n = K;
memset(tmp,0,sizeof(tmp));
pinv[0] = qpow(f1[0],mod-2);//求出乘法逆元的第一项
for(int len=2; len<=n; len<<=1)
{
K = len<<1;//两个多项式相乘后的长度
memcpy(tmp, f1, sizeof(LL) * len); //临时复制过来f(x)
memset(tmp+len, 0, sizeof(LL) * len); //后面补零
NTT(tmp,K,1); NTT(pinv,K,1);
for(int i=0; i<K; i++) pinv[i] = M(pinv[i], msub(2, M(tmp[i], pinv[i]))); //计算刚刚得出的式子
NTT(pinv,K,-1);//变换回来
memset(pinv+len, 0, sizeof(LL)*len); //取模，去掉多余的部分
//此轮循环结束后，pinv[]中存储的由 f_0^{-1}(x) 变为 f^{-1}(x)
}
memcpy(f1, pinv, sizeof(LL)*n);
return;
}
}
§ 1.2 哇塞，这居然可以
既然我们有了牛顿迭代，为什么不用用它呢？
好！那么我们先构造出这么一个函数：
\[\mathcal F[g(x)] = \frac{1}{g(x)} - f(x) \]可以看出来， \(\mathcal F\) 在模 \(x^n\) 意义下的零点就是我们要求的乘法逆元。
接下来列出牛顿迭代的式子：
\[g_{t+1}(x) \equiv g_{t}(x) - \frac{ \mathcal F[g_{t}(x)] }{ \mathcal F'[g_{t}(x)] } \pmod{x^n} \]带入解析式：
\[g_{t+1}(x) \equiv g_{t}(x) - \frac{ \frac{1}{g_t(x)} - f(x) }{ \mathcal F'[g_{t}(x)] } \pmod{x^n} \]Part 2 多项式除法（取模）
模板代码（需要刚才的NTT和Inverse）
LL copy1[maxn],copy2[maxn];
void Divide(LL *f1,LL *f2,int n1,int n2,LL *Q,LL *R)
{
memcpy(copy1,f1,sizeof(LL)*n1); memcpy(copy2,f2,sizeof(LL)*n2); //备份
reverse(f1,f1+n1); reverse(f2,f2+n2); //两个乘数翻转
int qlen = n1-n2+1; //模数的次数
if(qlen < n2) memset(f2+qlen, 0, sizeof(LL) * (n2-qlen));
Inverse(f2,qlen); //求除数的乘法逆元
Multiply(f1,f2,qlen,qlen); //f1乘上f2的逆元
memset(f1+qlen, 0, sizeof(LL) * qlen); reverse(f1, f1+qlen); //乘积反转后就是商
memcpy(Q, f1, sizeof(LL) * qlen); //因为下面要用到f1，这里先复制一遍
memcpy(f2, copy2, sizeof(LL) * (n2)); //用copy2之前先把f2复制回来
Multiply(copy2, f1, n2, qlen); //用除数乘上商
for(int i=0; i<n2-1; i++) R[i] = msub(copy1[i],copy2[i]); //用被除数减去上一步的得数得到余数
memcpy(f1, copy1, sizeof(LL) * n1); //最后把被除数复制回来
return;
}