[干饭现场] 多项式全家桶
- -1. 无关紧要的话
- 0. 泰勒展开与多项式牛顿迭代
- 0.1. 泰勒展开
- 0.1.1. 展开方式
- 0.1.2. 一些常见的展开
- 0.2. 多项式牛顿迭代
- 0.2.1. 问题描述
- 0.2.2. 牛顿迭代
- 0.1. 泰勒展开
- 1. 多项式求逆
- 1.1. 问题描述
- 1.2. 解法
- 1.3. 代码实现
- 2. 多项式开根
- 2.1. 问题描述
- 2.2. 解法
- 2.3. 代码实现
- 3. 多项式求 \(\rm ln\)
- 3.1. 问题描述
- 3.2. 解法
- 3.3. 代码实现
- 4. 多项式求 \(\rm exp\)
- 4.1. 问题描述
- 4.2. 解法
- 4.3. 代码实现
- 5. 多项式求幂
- 5.1. 问题描述
- 5.2. 解法
- 5.3. 代码实现
- To be continued...
-1. 无关紧要的话
\(\text{kfc}\) 全家桶都只够两个人吃，我这个全家桶不收录所有模板也不过分吧??????。
顺便推销一波自己的 和 .
0. 泰勒展开与多项式牛顿迭代
\(\color{red}{\text{Disclaimer}}\)：由于辣鸡博主数学很烂，所以可能会有很多不严谨的地方，如果哪些地方有问题请一定教教我 \(\text{sto}\).
0.1. 泰勒展开
0.1.1. 展开方式
泰勒展开是将一种在 \(x=x_0\) 处具有 \(n\) 阶导数的函数 \(f(x)\) 利用关于 \((x-x_0)\) 的 \(n\) 次来逼近函数多项式的方法。
若函数 \(f(x)\) 在包含 \(x_0\) 的某个闭区间 \([a,b]\) 上具有 \(n\) 阶导数，且在开区间 \((a,b)\) 上具有 \((n+1)\) 阶导数，则对 \([a,b]\) 上任意一点 \(x\)，成立下式
\[\begin{aligned} f(x)=\sum _{i=0}^{n}\frac{f^{(i)}(x_0)}{i!}\cdot (x-x_0)^i\end{aligned}+R_n(x) \]等号后的多项式称为函数 \(f(x)\) 在 \(x_0\) 处的泰勒展开式，剩余的 \(R_n(x)\) 是泰勒公式的余项，是 \((x-x_0)^n\) 的高阶无穷小，相当于泰勒展开式对原函数 \(f(x)\) 的一些相对而言极小的误差。每一次展开都是一次逼近。
当然，如果函数可以无限次求导的话，泰勒展开也可以无限次展开。但事实上对于我们所研究的有限次多项式，经过有限次求导就会变成 \(0\)，展开次数仍然是有限的。
0.1.2. 一些常见的展开
-
\(\sin x\) 的 \(0,1,2,3\) 阶导数为 \(\sin x,\cos x,-\sin x,-\cos x\)。这是一个循环。
取 \(x_0=0\)，泰勒展开有
\[\begin{align} \sin x&=\sum_{i=0}^{\infty}\frac{\sin^{(i)}(0)}{i!}\cdot x^i\\ &=\sum_{k=0}^{\infty}\frac{\sin^{(2k)}(0)}{2k!}\cdot x^{2k}+\sum_{k=0}^{\infty}\frac{\sin^{(2k+1)}(0)}{(2k+1)!}\cdot x^{2k+1}\\ &=\sum_{k=0}^{\infty}\frac{(-1)^k}{(2k+1)!}\cdot x^{2k+1} \end{align} \] -
\(\ln (1+x)= \displaystyle \sum_{n\geqslant 1}(-1)^{n+1}\cdot \dfrac{x^n}{n}\)，同时可以发现 \((1+x)\) 可以很好地消去 \(n=0\) 的情况。
同样地，\(\ln (1-x)=\ln (1+(-x))=- \displaystyle \sum_{n\geqslant 1}\dfrac{x^n}{n}\).
-
其实是 \(\ln(1+x)\) 展开的拓展，把系数移了一位，\(\dfrac{1}{1+x}= \displaystyle \sum_{n\geqslant 0} (-1)^nx^n\).
-
还可以用泰勒展开证明二项式定理：\((1+x)^a= \displaystyle \sum_{n\geqslant 0}\text{C}(a,n)\cdot x^n\).
-
另外还有：\(\dfrac{1}{(1-x)^{a+1}}= \displaystyle \sum_{n\geqslant 0}\dbinom{a+n}{n}\cdot x^n\)，特别地，当 \(a=1\) 时存在 \(\dfrac{1}{(1-x)^2}= \displaystyle \sum_{n\geqslant 0}(n+1)\cdot x^n\).
-
\(\text{e}^x= \displaystyle \sum_{n\geqslant 0} \dfrac{x^n}{n!}\)，另外 \({\rm e}^{px}= \displaystyle \sum_{n\geqslant 0}\dfrac{p^nx^n}{n!}\).
0.2. 多项式牛顿迭代
0.2.1. 问题描述
给定多项式 \(g(x)\)，已知有 \(f(x)\) 满足
\[g(f(x))\equiv 0\pmod{x^n} \]求出模 \(x^n\) 意义下的 \(f(x)\).
其中 \(g(f(x))=\sum_{i=0}^{\infty}g_i\cdot f^i(x)\)。需要特别注意的是，\(g\) 是关于 \(f(x)\) 的函数，而不是关于 \(x\) 的函数。
0.2.2. 牛顿迭代
考虑倍增。首先当 \(n=1\) 时，需要预先算出对应 \(f(x)\)。记 \(f_0(x)\) 为在模 \(x^{n/2}\) 意义下的解，那么有
\[\begin{align} g(f_0(x))&\equiv 0\pmod{x^{n/2}}\\ g(f(x))&\equiv 0\pmod{x^{n/2}} \end{align} \]移项得
\[f(x)-f_0(x)\equiv 0\pmod{x^{n/2}} \]后面有类似的证明，这里不再赘述，反正可以得到
\[(f(x)-f_0(x))^2\equiv 0\pmod{x^n} \]\(\rm ok\) 接下来将 \(g(f(x))\) 在 \(f_0(x)\) 处泰勒展开[1]
\[g(f(x))=\sum_{i=0}^{\infty}\frac{g^{(i)}(f_0(x))}{i!}\cdot (f(x)-f_0(x))^i \]由上文可知，\(i\ge 2\) 的项都可以被舍去[2]，解一下方程即可得到
\[f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\pmod{x^n} \]1. 多项式求逆
1.1. 问题描述
给定多项式 \(f\)，求多项式 \(g\) 满足 \(f*g\equiv 1\ (\text{mod}\ x^n)\).
1.2. 解法
先开始属实没看懂这是什么意思。其实翻译过来就是卷积后满足前 \(n\) 项系数除了常数项为 \(1\)，其余均为 \(0\).
可以考虑倍增求解。若已求得 \(g'\) 满足
\[f*g'\equiv 1\ (\text{mod}\ x^{\frac{n}{2}}) \]显然要求的 \(g\) 也要满足
\[f*g\equiv 1\ (\text{mod}\ x^{\frac{n}{2}}) \] \[\begin{align} f*g*g&\equiv f*g'*g\ (\text{mod}\ x^{\frac{n}{2}})\\ g&\equiv g'\ (\text{mod}\ x^{\frac{n}{2}})\\ g-g'&\equiv 0\ (\text{mod}\ x^{\frac{n}{2}}) \end{align} \]将多项式 \(g-g'\) 与其本身做卷积，由于其 \(0\) 到 \(n/2-1\) 项系数都为 \(0\)，容易发现卷积完 \(0\) 到 \(n-1\) 项系数都为 \(0\)（考虑 \(n-1\) 最平均的分解也是 \((n/2-1)+(n/2)\) 了）。另外，\(n/2\) 需要 向上取整。
\[(g-g')^2\equiv 0\ (\text{mod}\ x^n) \]拆项有
\[\begin{align} g^2-2gg'+(g')^2&\equiv 0\ (\text{mod}\ x^{n})\\ g-2g'+(g')^2*f&\equiv 0\ (\text{mod}\ x^{n})\\ g&\equiv 2g'-(g')^2*f \ (\text{mod}\ x^n) \end{align} \]递归求解即可，复杂度 \(\mathcal O(n\log n)\)。这是因为里面层数加起来的复杂度也不及最外一层的复杂度。
上面就是朴素地推导多项式求逆的过程，接下来以其为例使用牛顿迭代再次推导。
假设 \(h(x)\) 是需要求逆的函数，\(f(x)\) 为它的逆，记 \(g(f(x))=\frac{1}{f(x)}-h(x)\)，容易发现这个式子在模 \(x^n\) 意义下是同余于 \(0\) 的。直接套用牛顿迭代得
\[\begin{align} f(x)&\equiv f_0(x)-\frac{\frac{1}{f_0(x)}-h(x)}{-\frac{1}{f_0^2(x)}}\pmod{x^n}\\ &\equiv 2f_0(x)-f_0^2(x)\cdot h(x)\pmod{x^n} \end{align} \]这里再说一下为啥 \(g(f_0(x))\) 求导是 \(-\frac{1}{f_0^2(x)}\)，这是因为 \(g\) 是关于 \(f_0(x)\) 的函数！所以 \(h(x)\) 相当于常数项。
可以发现，我们只需要配凑一个关于原函数和目标函数同余零的式子，同时将目标函数设为 \(g\) 的自变量，就可以使用牛顿迭代求解多项式的各类运算。
1.3. 代码实现
先开始想着写递推版会快一些，结果跑出来发现比之前写的递归版慢了一倍??????。后来发现是因为递归版有个很大的优化 —— 从给定长度 \(n\) 开始除以二。而递推版是从 \(2\) 开始倍增，当长度不小于 \(n\) 时才开始卷，此时的上界就变成 \(8n\) 了！虽然改完还是要跑五百多毫秒。
另外要注意的代码细节是从 \(2\) 开始卷而不是从 \(1\) 开始卷，因为从 \(g'\) 递推到 \(g\) 是在模 \(x^n\) 意义下而不是在模 \(x^{n/2}\) 意义下完成的。
# include
# include
# define print(x,y) write(x), putchar(y)
template
inline T read(const T sample) {
T x=0; char s; bool f=0;
while(!isdigit(s=getchar())) f|=(s=='-');
for(; isdigit(s); s=getchar()) x=(x<<1)+(x<<3)+(s^48);
return f? -x: x;
}
template
inline void write(T x) {
static int writ[50], w_tp=0;
if(x<0) putchar('-'), x=-x;
do writ[++w_tp]=x-x/10*10, x/=10; while(x);
while(putchar(writ[w_tp--]^48), w_tp);
}
# include
# include
using namespace std;
const int maxn = 4e5+5;
const int mod = 998244353, g = 3;
inline int inv(int x,int y=mod-2,int r=1) {
for(; y; y>>=1, x=1ll*x*x%mod)
if(y&1) r=1ll*r*x%mod;
return r;
} const int ig = inv(g);
inline int dec(int x,int y) { return x-y<0?x-y+mod:x-y; }
inline int inc(int x,int y) { return x+y>=mod?x+y-mod:x+y; }
int rev[maxn],n,F[maxn],G[maxn];
inline void ntt(int* f,int lim,bool opt=0) {
static int w, wn, tmp;
for(int i=0;i>1; g[0] = inv(h[0]);
for(; b_len; --b_len) {
w=bit[b_len], B=0, lim=1;
while(lim<(w<<1)) lim<<=1, ++B;
memcpy(f,h,sizeof(int)*w);
memset(f+w,0,sizeof(int)*(lim-w));
for(int i=0;i>1]>>1)|((i&1)< 2. 多项式开根
2.1. 问题描述
给定多项式 \(f\)，求多项式 \(g\) 满足 \(g^2\equiv f\ (\text{mod}\ x^n)\).
2.2. 解法
可以考虑倍增求解。若已求得 \(g'\) 满足
\[(g')^2\equiv f\ (\text{mod}\ x^{\frac{n}{2}}) \]显然要求的 \(g\) 也要满足
\[g^2\equiv f\ (\text{mod}\ x^{\frac{n}{2}}) \]故（由平方差公式）
\[g-g'\equiv 0\ (\text{mod}\ x^{\frac{n}{2}}) \]同上可得
\[(g-g')^2\equiv 0\ (\text{mod}\ x^n) \]拆项有
\[g^2-2gg'+(g')^2\equiv 0\ (\text{mod}\ x^{n}) \]\[f-2gg'+(g')^2\equiv 0\ (\text{mod}\ x^{n}) \]\[g\equiv \frac{f+(g')^2}{2g'} \ (\text{mod}\ x^n) \]复杂度同样也是 \(\mathcal O(n\log n)\)，因为求逆与开根的运算实际上是并行的。
2.3. 代码实现
# include
# include
# define print(x,y) write(x), putchar(y)
template
inline T read(const T sample) {
T x=0; char s; bool f=0;
while(!isdigit(s=getchar())) f|=(s=='-');
for(; isdigit(s); s=getchar()) x=(x<<1)+(x<<3)+(s^48);
return f? -x: x;
}
template
inline void write(T x) {
static int writ[50], w_tp=0;
if(x<0) putchar('-'), x=-x;
do writ[++w_tp]=x-x/10*10, x/=10; while(x);
while(putchar(writ[w_tp--]^48), w_tp);
}
# include
# include
using namespace std;
typedef long long ll;
const int maxn = 4e5+5;
const int mod = 998244353, g = 3;
inline int inv(int x,int y=mod-2,int r=1) {
for(; y; y>>=1, x=1ll*x*x%mod)
if(y&1) r=1ll*r*x%mod;
return r;
} const int ig = inv(g), iv2 = inv(2);
inline int dec(int x,int y) { return x-y<0?x-y+mod:x-y; }
inline int inc(int x,int y) { return x+y>=mod?x+y-mod:x+y; }
int rev[maxn],n,F[maxn],G[maxn];
inline void ntt(int* f,int lim,bool opt=0) {
static int tmp, w, wn;
for(int i=0;i>1; g[0]=inv(h[0]);
for(; b_len; --b_len) {
w=bit[b_len], B=0, lim=1;
while(lim<(w<<1)) lim<<=1, ++B;
for(int i=0;i>1]>>1)|((i&1)<>1; g[0]=1;
for(; b_len; --b_len) {
w=bit[b_len], B=0, lim=1;
while(lim<(w<<1)) lim<<=1, ++B;
memcpy(f,h,sizeof(int)*w);
memset(f+w,0,sizeof(int)*(lim-w));
memset(_g,0,sizeof(int)*lim);
polyInv(_g,g,w); // 没有预处理 rev[] 的原因是在 polyInv() 中已经处理过了
ntt(f,lim), ntt(_g,lim);
for(int i=0;i 3. 多项式求 \(\rm ln\)
3.1. 问题描述
给定多项式 \(f\)，求多项式 \(g\) 满足 \(g\equiv \ln(f)\ (\text{mod}\ x^n)\).
3.2. 解法
记 \(h(x)=\ln x\)，那么其实就是求
\[g\equiv h(f)\pmod{x^n} \]两边同时求导得[3]
\[g'\equiv \frac{f'}{f}\pmod{x^n} \]在求出 \(g'\) 后，再求一次积分即可
\[\int x^a\ {\rm d}x=\frac{1}{a+1}\cdot x^{a+1} \]另外还有一个需要注意的地方：
在模意义下 当且仅当 \([x^0]f(x)=1\) 时，\(f(x)\) 才有对数多项式，且常数项为零；当且仅当 \([x^0]f(x)=0\) 时，\(f(x)\) 才有指数多项式，且常数项为 \(1\).
解释一下为啥是这样（当然是感性的）。考虑 \(\ln f(x)\) 在 \(x_0=0\) 处进行泰勒展开时，第一项是 \(\dfrac{\ln f_0}{0!}\cdot x^0\)，若 \(f(x)\) 的常数项不为 \(1\)，首先 \(<1\) 肯定是未定义的，对于 \(>1\) 的情况，可以利用 \(\ln(x+1)\) 来观察，当 \(x>0\) 时这个函数在模意义下是不收敛的，所以没法算。\(\exp\) 也是同样的道理，这里就不讲了。
3.3. 代码实现
# include
# include
# define print(x,y) write(x), putchar(y)
template
inline T read(const T sample) {
T x=0; char s; bool f=0;
while(!isdigit(s=getchar())) f|=(s=='-');
for(; isdigit(s); s=getchar()) x=(x<<1)+(x<<3)+(s^48);
return f? -x: x;
}
template
inline void write(T x) {
static int writ[50], w_tp=0;
if(x<0) putchar('-'), x=-x;
do writ[++w_tp]=x-x/10*10, x/=10; while(x);
while(putchar(writ[w_tp--]^48), w_tp);
}
# include
# include
using namespace std;
const int maxn = 4e5+5;
const int mod = 998244353, g = 3;
int inv(int x,int y=mod-2,int r=1) {
for(; y; y>>=1, x=1ll*x*x%mod)
if(y&1) r=1ll*r*x%mod; return r;
} const int ig = inv(g);
int dec(int x,int y) { return x-y<0?x-y+mod:x-y; }
int inc(int x,int y) { return x+y>=mod?x+y-mod:x+y; }
int rev[maxn],F[maxn],G[maxn],n;
void ntt(int* f,int lim,bool opt=0) {
int tmp, w, wn;
for(int i=0;i>1; g[0]=inv(h[0]);
for(; b_len; --b_len) {
w=bit[b_len], lim=1, B=0;
while(lim<(w<<1)) lim<<=1, ++B;
for(int i=0;i>1]>>1)|((i&1)<=0;--i)
f[i+1] = 1ll*f[i]*inv(i+1)%mod;
f[0]=0;
}
void polyLn(int* g,int* h,int n) {
static int lim=1, Inv;
while(lim<(n<<1)) lim<<=1;
polyInv(g,h,n); dera(h);
ntt(g,lim), ntt(h,lim);
for(int i=0;i 4. 多项式求 \(\rm exp\)
4.1. 问题描述
给定多项式 \(h(x)\)，求多项式 \(f(x)\) 满足 \(f(x)\equiv {\rm e}^{h(x)}\ (\text{mod}\ x^n)\).
4.2. 解法
由于之前已经解决了多项式求 \(\ln\)，所以不妨将这个玩意转化一下
\[\ln f(x)\equiv h(x)\pmod{x^n} \]使用牛顿迭代来解决这个问题，记 \(g(f(x))=\ln f(x)-h(x)\)，由于它在模 \(x^n\) 意义下是为零的，所以有
\[\begin{align} f(x)&=f_0(x)-\frac{\ln f_0(x)-h(x)}{g'(f_0(x))}\\ &=f_0(x)-\frac{\ln f_0(x)-h(x)}{\frac{1}{f_0(x)}}\\ &=f_0(x)\cdot (1-\ln f_0(x)+h(x)) \end{align} \]由于一切操作都是并行的，所以复杂度仍然是 \(\mathcal O(n\log n)\) 的！哦！让我们高声地赞美牛顿迭代法吧！常数大得飞起就是了。
4.3. 代码实现
# include
# include
# define print(x,y) write(x), putchar(y)
template
inline T read(const T sample) {
T x=0; char s; bool f=0;
while(!isdigit(s=getchar())) f|=(s=='-');
for(; isdigit(s); s=getchar()) x=(x<<1)+(x<<3)+(s^48);
return f? -x: x;
}
template
inline void write(T x) {
static int writ[50], w_tp=0;
if(x<0) putchar('-'), x=-x;
do writ[++w_tp]=x-x/10*10, x/=10; while(x);
while(putchar(writ[w_tp--]^48), w_tp);
}
# include
# include
using namespace std;
const int maxn = 4e5+5;
const int mod = 998244353, g = 3;
int inv(int x,int y=mod-2,int r=1) {
for(; y; y>>=1, x=1ll*x*x%mod)
if(y&1) r=1ll*r*x%mod; return r;
} const int ig = inv(g);
int dec(int x,int y) { return x-y<0?x-y+mod:x-y; }
int inc(int x,int y) { return x+y>=mod?x+y-mod:x+y; }
int rev[maxn],F[maxn],G[maxn],n;
void ntt(int* f,int lim,bool opt=0) {
int tmp, w, wn;
for(int i=0;i>1; g[0]=inv(h[0]);
for(; b_len; --b_len) {
w=bit[b_len], lim=1, B=0;
while(lim<(w<<1)) lim<<=1, ++B;
for(int i=0;i>1]>>1)|((i&1)<=0;--i)
f[i+1] = 1ll*f[i]*inv(i+1)%mod;
f[0]=0;
}
void polyLn(int* g,const int* h,int n) {
static int lim=1, Inv, f[maxn];
while(lim<(n<<1)) lim<<=1;
memset(g,0,sizeof(int)*lim);
memcpy(f,h,sizeof(int)*lim);
polyInv(g,h,n); dera(f,n);
ntt(g,lim), ntt(f,lim);
for(int i=0;i>1; g[0]=1;
for(; b_len; --b_len) {
w=bit[b_len], lim=1, B=0;
while(lim<(w<<1)) lim<<=1, ++B;
polyLn(lng,g,w); lng[0] = dec(lng[0],1);
for(int i=0;i 5. 多项式求幂
5.1. 问题描述
给定多项式 \(h(x)\)，求多项式 \(f(x)\) 满足 \(f(x)\equiv (h(x))^k\ (\text{mod}\ x^n)\).
5.2. 解法
先取对数
\[\ln f(x)\equiv k\ln h(x)\pmod{x^n} \]再取一下指数
\[f(x)\equiv {\rm e}^{k\ln h(x)}\pmod {x^n} \]5.3. 代码实现
Believe in yourself.
To be continued...
为了方便这里直接按无限展开来写，之后再说为什么。 ??
这就是为啥可以直接无限展开，因为保留的项数是很少的。同时也可以发现我们对 \(g(x)\) 的项数存在一定限制。 ??
这似乎是 \(\rm ln\) 的一个通用方法，将求 \(\rm ln\) 变成求逆元。 ??