第一次用博客园写学习笔记，写的不好请见谅~

欢迎各位OIer在评论区说一下自己对这篇博客的建议！

有什么问题也欢迎在下方评论区中提出！顺便点赞

有关多项式的基本概念

对于求和式 \(\sum a_nx^n\) ，如果是有限项相加，称为多项式，记作 \(f(x)=\sum_{k=0}^n a_kx^k\)。

一个多项式的度（也称次数）就是这个多项式最高次项的次数，记作 \(\deg f\) 。

多项式的表示一般有两种方法：

系数表示，就是形如 \(f(x)=\sum_{k=0}^n a_kx^k\) 的形式。

点值表示，即给出 \(n+1\) 个点 \((x_{0}, y_{0}), (x_{1}, y_{1}), \dots, (x_{n}, y_{n})\) ，求一个 \(n\) 次多项式使得这 \(n+1\) 个点都在 \(f(x)\) 上。

定义加法和乘法以后，多项式的基本运算和整数基本一致。（这个大家都会吧）

多项式的导数公式：

\[\left(\sum_{k=0}^{+\infty}a_kx^k\right)'=\sum_{k=1}^{+\infty}ka_kx^{k-1}

\]

多项式的积分公式：

\[\int_0^x\left(\sum_{k=0}^{+\infty}a_kx^k\right)\text{d}x=\sum_{k=0}^{+\infty}\frac{a_kx^{k+1}}{k+1}

\]

一个长度为 \(n+1\) 的序列 \(a_i\) 的生成函数定义为 \(a(x)=\sum_{k=0}^n a_kx^k\) 。（以下全文中 \(f_i\) 的生成函数记为 \(f(x)\) ）

多项式乘法

单位根（前置知识）

前置知识：复数

顾名思义，单位根就是满足 \(x^n=1\) 的根。

众所周知，在复数域内， \(n\) 次多项式有 \(n\) 个根，因此 \(x^n=1\) 也有 \(n\) 个根，一般把第 \(k\) 个根记作 \(\omega_n^k\) 。

由于复数乘法的几何意义是模长相乘，辐角相加，因此 \(\omega_n^k\) 的模长一定为 \(1\) ，辐角为 \(\frac{2\pi}{n}\) 的倍数。

为方便，我们设 \(\omega_n^k\) 的辐角为 \(\frac{2k\pi}{n}\) ，或者 \(\frac{360^\circ k}{n}\) 。

单位根有三个重要的性质。对于任意正整数 \(n\) 和整数 \(k\)：

\[\begin{aligned}

\omega_n^n&=1\\

\omega_n^k&=\omega_{xn}^{xk}\\

\omega_{2n}^{k+n}&=-\omega_{2n}^k\\

\end{aligned}

\]

这三个性质在快速傅里叶变换中会被用到。

顺便一提，单位根还有一个比较重要的性质，是： \(\sum_{i=0}^{n-1}(\omega_n^k)^i=0\) 。

DFT

又称离散傅里叶变换、Discrete Fourier Transform。

DFT的思想是将两个多项式的系数表示与点值表示互相转化。

因为多项式的系数表示直接相乘是 \(O(n^2)\) 的，而点值表示相乘是 \(O(n)\) 的。（只需要把每个点处的值相乘）

先考虑将系数表示转化为点值。

如果暴力去求，时间复杂度还是 \(O(n^2)\) 的。（这个过程就被称为DFT）

那这个时候我们就要用到一个东西叫做：

FFT

又称快速傅里叶变换、Fast (Discrete) Fourier Transform。

我们先把多项式的次数 \(n\) 增加至 \(2^k-1\) 的形式。（一会儿我就会说明为什么要这样做了）

接下来，我们先尝试将系数表示转化为点值表示，这一步就被称之为FFT。

我们考虑一个 \(n\) 次多项式有什么性质：

以 \(n=7\) 为例：

\[f(x) = a_0 + a_1x + a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7

\]

我们先按照次数的奇偶来分成两组，然后右边提出来一个 x：

\[\begin{aligned}

f(x) &= (a_0+a_2x^2+a_4x^4+a_6x^6) + (a_1x+a_3x^3+a_5x^5+a_7x^7)\\

&= (a_0+a_2x^2+a_4x^4+a_6x^6) + x(a_1+a_3x^2+a_5x^4+a_7x^6)

\end{aligned}

\]

然后再用奇偶次次项数建立新的函数：

\[\begin{aligned}

f_0(x) &= a_0+a_2x+a_4x^2+a_6x^3\\

f_1(x) &= a_1+a_3x+a_5x^2+a_7x^3

\end{aligned}

\]

那么原来的 \(f(x)\) 就可以表示为：

\[f(x)=f_0\left(x^2\right) + x \times f_1\left(x^2\right)

\]

于是我们可以发现一个性质：它可以分治！

也就是说，我们可以先将原来的式子的奇数项和偶数项分别处理（即分别计算FFT），求出原文中的 \(f_0\) 和 \(f_1\) ，再代入式子求解。

那怎么合并呢？

由于 \(f_1\) 和 \(f_2\) 都只有 \(n/2\) 个点值，那么我们如果用以上那个式子的话，也只能推出 \(f\) 的 \(n/2\) 个点值。

但是，我们可以再考虑以下这个式子：

\[f(-x)=f_0\left(x^2\right) - x \times f_1\left(x^2\right)

\]

这样的话，我们能用这个式子求出另外 \(n/2\) 个点值。

你以为这就完了吗？还没有！

观察上面两个式子，我们发现：它只能求互为相反数的 \(n/2\) 对 \(x\) ，而如果在实数范围内考虑，右边就要求 \(n/2\) 个正数的点值，显然行不通。

于是我们考虑单位根。

单位根有很好的的对称性，并且满足平方还是单位根这一特性。（因为 \((\omega_n^k)^2 = \omega_{n/2}^k\) ）

如果我们把 \(n\) 个单位根代入，那么式子就会变成这样：

\[\begin{aligned}

f(\omega_n^k) &= f_0(\omega_{n/2}^k) + \omega_n^k f_1(\omega_{n/2}^k)\\

f(\omega_n^{k+n/2}) &= f_0(\omega_{n/2}^k) - \omega_n^k f_1(\omega_{n/2}^k)

\end{aligned}

\]

（注意到之前我们让 \(n\) 扩大到 \(2^k-1\) 的形式，就是将项数扩大到 \(2^k\) ，方便减半）

于是递归即可。

边界条件： \(n=1\) 时，\(f(x)=a_0\) 。

时间复杂度：\(T(n)=2T(\frac n2)+O(n)\) ，化简得 \(O(n\log n)\) 。

逆FFT(IFFT)

现在我们已经解决了系数表示转换为点值表示，那么我们怎么将点值表示转换为系数表示呢？

我们知道有一个东西叫做高斯消元逆FFT。

它的证明需要依赖矩阵，这里不详细说明。（主要是作者也不是很会）

简单来说，你只需要把之前FFT的过程变成这样：

\[\begin{aligned}

f(\omega_n^k) &= f_0(\omega_{n/2}^k) + \omega_n^{-k} f_1(\omega_{n/2}^k)\\

f(\omega_n^{k+n/2}) &= f_0(\omega_{n/2}^k) - \omega_n^{-k} f_1(\omega_{n/2}^k)

\end{aligned}

\]

（就是系数变成它的共轭复数）

然后最终得出来的结果除以 \(n\) 即可。

好的，到了这一步，我们已经会多项式乘法的基本操作了。

但是别着急，后面还有一个优化：

非递归写法

虽然 FFT 和 IFFT 都是 \(O(n\log n)\) 的，但是它们有一个致命的缺陷：常数太大。

那我们应该怎么解决这个问题呢？

我们观察到 FFT 需要大量的复制数组的操作，这使得它的常数非常大。

有没有通过不移动数就能求出 FFT 的方法呢？

当然有。

我们先把整个过程压到一个数组中，再考虑之前的 FFT 都把数移到了哪些位置：（还是以 \(n=7\) 举例，同一个数组的用中括号括起来）

下标

0

1

2

3

4

5

6

7

第一层

[0

1

2

3

4

5

6

7]

第二层

[0

2

4

6]

[1

3

5

7]

第三层

[0

4]

[2

6]

[1

5]

[3

7]

边界

[0]

[4]

[2]

[6]

[1]

[5]

[3]

[7]

可以看出最后一层每一位都是这个位置下标的二进制翻转过来。

而我们如果最开始把边界数组存下来，再将固定距离的数进行合并即可。

于是我们可以先预处理出每个数的二进制翻转后的结果 \(rev_i\) ，顺便求出边界情况的数组 \(b_i = a_{{rev}_i}\) 。

然后从最下面一层开始枚举，每层距离变为原来的 \(2\) 倍，并对同一个区域内这个距离的数进行合并。（类似后缀排序的归并）

这样直接合并，时间复杂度仍为 \(O(n\log n)\) ，但常数小了不少。

代码

#include

using namespace std;

#define int long long

const double Pi=acos(-1);//用这个式子算pi要好一些

inline int read()

{

char ch=getchar();int x=0,r=1;

while(ch<'0'||ch>'9'){if(ch=='-')r=0;ch=getchar();}

while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+ch-'0',ch=getchar();

return r?x:-x;

}

struct complexs

{

double x,y;

complexs(double xx=0,double yy=0){x=xx;y=yy;}

complexs operator +(double b){return {x+b,y};}

complexs operator -(double b){return {x-b,y};}

complexs operator *(double b){return {x*b,y*b};}

complexs operator /(double b){return {x/b,y/b};}

complexs operator +(complexs b){return {x+b.x,y+b.y};}

complexs operator -(complexs b){return {x-b.x,y-b.y};}

complexs operator *(complexs b){return {x*b.x-y*b.y,x*b.y+y*b.x};}

complexs operator /(complexs b){return (complexs){x*b.x+y*b.y,-x*b.y+y*b.x}/(b.x*b.x+b.y*b.y);}

};//虽然看上去很复杂但实际上就是简单的复数运算

typedef vector vc;//vc就是vector

int rev[2100000];//二进制翻转数组

void FFT(vc &x,int limit,bool type=1)//type=1是FFT，否则是IFFT

{

for(int i=0;i

for(int len=1;len

{

complexs Wn(cos(Pi/len),(2*type-1)*sin(Pi/len));//单位根，注意IFFT的时候单位根要取反

for(int i=0;i

{

complexs W(1,0),X,Y;

for(int j=i;j

{

X=x[j];Y=W*x[j+len];

x[j]=X+Y;x[j+len]=X-Y;//合并

W=W*Wn;

}

}

}

if(!type)for(int i=0;i<=limit;++i)x[i]=x[i]/limit;//IFFT时每个数要除以n

}

int limit;

void operator *=(vc &a,vc b)

{

int x=a.size()+b.size()-1;

limit=1;

while(1ull*limit<=a.size()-1+b.size()-1)limit<<=1;//注意limit是要大于等于项数，不是大于等于次数，项数=次数+1

for(int i=0;i

while(a.size()<1ull*limit)a.push_back(0);

while(b.size()<1ull*limit)b.push_back(0);

FFT(a,limit);FFT(b,limit);//系数转换成点值

for(int i=0;i

FFT(a,limit,0);//点值转换成系数

while(a.size()>1ull*x)a.pop_back();

}

vc a,b;

int n,m;

signed main()

{

n=read();m=read();

for(int i=0;i<=n;++i)a.push_back({1.0*read(),0});

for(int i=0;i<=m;++i)b.push_back({1.0*read(),0});

a*=b;

for(complexs i:a)printf("%lld ",(int)(i.x+0.5));//精度丢失较大，建议四舍五入

return 0;

}

三次变两次优化

省流：一般不用。

假设我们要求 \(f(x)*g(x)\) ，考虑 \(p(x)=f(x)+ig(x)\) 。

于是 \(p^2(x)=f^2(x)-g^2(x)+2i*f(x)g(x)\) 。

发现它的虚部除以 \(2\) 就是答案，输出即可。

NTT

又称快速数论变换、Number-Theoretic Transform。（在有的地方也被称为FNTT？总之都是这个意思）

虽然FFT好用，但它也有一些致命的弱点，例如复数乘法和三角函数常数太大，精度不高，不支持取模等等。

那我们能解决这些问题吗？不能

当然可以，于是我们就有了 NTT。

回忆一下FFT之所以用单位元是用了它的哪些性质。

\(\forall 0\le i,j

\(\omega_{kn}^{km}=\omega_n^m\)（倍增时）

\((\omega_n^m)^2=(\omega_n^{m+n/2})^2=\omega_n^{2m}\)（合并时）

\(\omega_1^0 =1\)（写起来方便）

于是我们发现数论中的原根也有这个性质。

前置知识：原根（你只要知道有这么个东西就行了，细节可以不用看）

更确切地说，我们设原根为 \(g\) ，令 \(\omega_n^m=g^\frac{(p-1)m}{n}\) 即可。

利用数论，它可以解决精度和取模问题。

其它代码基本和 FFT 一样。

注意：NTT不能用三次变两次优化！

证明

第一条由原根的定义，只要原根的指数不同那么值也不同。

由于 \(\frac{p-1}n\) 总是整数（并且 \(n\) 一定时是常数），且 \(m\) 值不一致，故成立。

第二条：\(\omega_{kn}^{km}=g^\frac{km(p-1)}{kn}=g^\frac{m(p-1)}{n}=\omega_n^m\)

第三条：\((\omega_n^m)^2=(g^\frac{m(p-1)}{n})^2=g^\frac{2m(p-1)}{n}=\omega_n^{2m}\)

并且 \((\omega_n^{m+n/2})^2=(g^\frac{(m+n/2)(p-1)}{n})^2=g^\frac{2m(p-1)}{n}\cdot g^{p-1}=\omega_n^{2m}\cdot 1=\omega_n^{2m}\)

证毕。

代码（用vector实现）

#include

using namespace std;

#define int long long

typedef vector vi;//vi指vector

const int mod=998244353,g=3,ig=332748118;//ig是g的逆元

inline int read()

{

char ch=getchar();int x=0,r=1;

while(ch<'0'||ch>'9'){if(ch=='-')r=0;ch=getchar();}

while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+ch-'0',ch=getchar();

return r?x:-x;

}

int mul(int x,int b)//快速幂

{

int ans=1;

while(b)

{

if(b&1)ans=(ans*x)%mod;

x=(x*x)%mod;

b>>=1;

}

return ans;

}

int rev[2100000];

void NTT(vi &x,int limit,bool type=1)

{

for(int i=0;i

for(int len=1;len

{

int Wn=mul(type?g:ig,(mod-1)/(2*len));//利用原根代替FFT中的单位元

for(int i=0;i

{

int W=1,X,Y;

for(int j=i;j

{

X=x[j];Y=W*x[j+len]%mod;

x[j]=(X+Y)%mod;x[j+len]=(X-Y+mod)%mod;//合并同FFT，但是要加上取模运算

W=(W*Wn)%mod;

}

}

}

if(!type)

{

int invs=mul(limit,mod-2);

for(int i=0;i

}

}

int limit;

void operator *=(vi &a,vi b)

{

int x=a.size()+b.size()-1;

limit=1;

while(1ull*limit<=a.size()-1+b.size()-1)limit<<=1;

for(int i=0;i

while(a.size()<1ull*limit)a.push_back(0);

while(b.size()<1ull*limit)b.push_back(0);

NTT(a,limit);NTT(b,limit);//FFT变为NTT

for(int i=0;i

NTT(a,limit,0);

while(a.size()>1ull*x)a.pop_back();

}

vi a,b;

int n,m;

signed main()

{

n=read();m=read();

for(int i=0;i<=n;++i)a.push_back(read());

for(int i=0;i<=m;++i)b.push_back(read());

a*=b;

for(int i:a)printf("%lld ",i);

return 0;

}

MTT

省流：解决任意模数的问题，一般不会用，但出现时比较恶心，可以写多次卷积+CRT（中国剩余定理），就是常数容易爆。

我们发现NTT只适用于模数为 \(a*2^n+1\) 的形式。（2^n要大于多项式项数，否则也会炸）

那如果模数不是 \(998244353\) 这样的和谐的数字，是 \(10^9+7\) 之类的呢？

解法一

假设原来的多项式值域是 \(m=10^9\) ，个数是 \(n=10^5\) （一般都是这个范围），那么能产生的最大数也不过就是 \(mn^2=10^{26}\) 。

考虑对多个模数卷积（ \(3\) 个就够了），再用CRT（中国剩余定理）合并， \(9\) 次FFT，能过一般题，但容易炸。

解法二

考虑拆系数。

我们先顺手把每个系数先模上 \(mod\) ，于是系数的值域一定 \(

把系数拆成 \(c=a*2^{15}+b\) 的形式。

把 \(a\) 和 \(b\) 拆成两个多项式，相乘的值域为 \(n*\sqrt m*\sqrt m=10^{14}\)

于是 \(c_1*c_2=(a_1*2^{15}+b_1)(a_2*2^{15}+b_2)=a_1a_2*2^{30}+(a_1b_2+a_2b_1)*2^{15}+b_1b_2\) 。

对这个式子的 \(a_1,a_2,b_1,b_2\) 分别做一次FFT，\(a_1a_2,a_1b_2,a_2b_1,b_1b_2\) 分别做一次IFFT，总共要用 \(8\) 次，常数还是比较大。

并且我们不能用可爱的NTT了，这意味着我们要面对FFT的精度问题（毕竟上界是 \(10^14\) ），建议写 \(\text{long double}\) 。

以及，千万不要相信自己背 \(\pi\) 的能力，建议写 \(\arccos (-1)\) 。

解法三

我们考虑在解法二的基础上进行优化。

我们现在有四个多项式 \(a_1,a_2,b_1,b_2\) ，我们要求 \(a_1a_2,a_1b_2,a_2b_1,b_1b_2\) 的值。

上面的三次变两次优化告诉我们：有的时候塞一个复多项式进去可以把常数变得更小。

不妨设 \(p=a_1+ia_2\) ， \(q=b_1+ib_2\) 。

则 \(pq = a_1b_1-a_2b_2+i(a_1b_2+a_2b_1)\) 。

再取 \(p'=a_1-ia_2\) ， 则 \(p'q = a_1b_1+a_2b_2+i(a_1b_2-a_2b_1)\)。

因此，我们对 \(p,p',q\) 做FFT，对 \(pq,p'q\) 做IFFT，剩下的就是一些基本的多项式加减操作了，总共用 \(5\) 次。

还是要写 \(\text{long double}\) 。

多项式的各种运算

多项式乘法逆（倍增法）

多项式乘法逆的定义：对于一个 \(n\) 次多项式 \(f(x)\) 来说， \(f(x)\) 的逆是满足 \(f(x)\cdot g(x)\equiv 1\pmod {x^{n+1}}\) 的多项式 \(g(x)\) 。

我们还是考虑递归。

首先便于计算，我们还是把 \(n\) 扩大到 \(2^k-1\) 的形式。（有些题解里不是这样写的，你也可以直接递归，但这样方便一些）

由于高次项的系数不会影响低次项的逆元，我们这步“扩大”可以直接补 \(0\) 。

为了下文叙述方便，我们记 \(m=x^{n+1},m'=x^{\frac {n+1}2}\) 。

我们可以先求出 \(f(x)\) 对 \(m'\) 的逆，记作 \(f_0^{-1}(x)\) 。

于是我们有：

\[\begin{aligned}

f(x)f^{-1}_0(x)&\equiv 1 &\pmod{m'}\\

f(x)f^{-1}(x)&\equiv 1 &\pmod{m'}\\

f^{-1}(x)-f^{-1}_0(x)&\equiv 0 &\pmod{m'}

\end{aligned}

\]

将最后一个式子平方可得：

\[f^{-2}(x)-2f^{-1}(x)f^{-1}_0(x)+f^{-2}_0(x)\equiv 0 \pmod{m}

\]

两边同乘 \(f(x)\) ：

\[f^{-1}(x)-2f_0^{-1}(x)+f(x)f^{-2}_0(x)\equiv 0 \pmod{m}

\]

移项得：

\[f^{-1}\left(x\right)\equiv f^{-1}_{0}\left(x\right)\left(2-f\left(x\right)f^{-1}_{0}\left(x\right)\right) \pmod{m}

\]

倍增即可。

时间复杂度： \(T(n)=T(\frac n2)+O(n \log n)\) ，即 \(T(n)=O(n \log n)\) 。

代码

//没错我又更了一堆插件

void read(int n,vi &a){for(int i=0;i<=n;++i)a.push_back(read());}//输入，n为次数，或者项数+1

void print(vi &a){for(int i:a)printf("%lld ",i);puts("");}//输出

void cut(vi &a,int x){while(a.size()>1ull*x)a.pop_back();while(a.size()<1ull*x)a.push_back(0);}//将a的系数个数强制设为x，注意系数个数和deg是两个不同的东西

void copy(vi &a,int x,vi &b)//将a的前x个系数复制到b，不足补0

{

if(1ull*x<=a.size()){b.clear();for(int i=0;i

else {b=a;while(b.size()<1ull*x)b.push_back(0);}

}

void operator +=(vi &a,int x){a.front()=(a.front()+x)%mod;}//多项式+常数

void operator +=(vi &a,vi &b)//多项式相加

{

while(a.size()

for(int i=0;1ull*i

}

void operator --(vi &a){for(int i=0;1ull*i

void operator -=(vi &a,int x){a.front()=(a.front()-x+mod)%mod;}//多项式-常数

void operator -=(vi &a,vi &b)//多项式相减

{

while(a.size()

for(int i=0;1ull*i

}

vi inv(vi &a)//逆多项式部分

{

int x=1;

vi ans={mul(a.front(),mod-2)},b,now;//显然a对于x^1的逆就是a0的逆元

while(1ull*x

{

x<<=1;//x为系数个数

copy(a,x,b);//将a的后几位存下来，避免ans每次直接与a相乘，影响时间复杂度

now=ans;ans*=b;ans-=2;--ans;ans*=now;//按式子进行计算，注意代码中--的定义与整数不同

cut(ans,x);//去除多余位

}

cut(ans,a.size());//让ans次数和a一样

return ans;

}

多项式除法

顾名思义，就是一个多项式 \(F(x)\) 除以另一个多项式 \(G(x)\) 的答案。

我们直接对 \(G(x)\) 取逆，再乘上 \(F(x)\) 即可。

好吧，真正的多项式除法实际上是带鱼带余除法。

或者说得正式点：给定一个 \(n\) 次多项式 \(f(x)\) 和一个 \(m\) 次多项式 \(g(x)\) ，请求出多项式 \(q(x)\), \(r(x)\)，满足以下条件：

\(\deg q(x)=n-m\)，\(\deg r(x)=m-1\)

\(f(x)=q(x)*g(x)+r(x)\)

要解决这样一个问题，我们首先构造这样一个变换：

\[f^R(x)=x^{\deg f}f(\frac{1}{x})

\]

观察到这个变换的实质是反转 \(f(x)\) 的系数，即 \(\text{vector}\) 当中的 \(\text{reverse}\) 函数。

并且实际上有 \((f^R)^R(x)=f(x)\) ，因此我们只需要求出 \(f^R(x)\) 就能求出 \(f(x)\) 。

（注： \(f^R(x)\) 和 \(f^{-1}(x)\) 是两个不同的东西，一个是反转系数，一个是乘法逆）

由 \(f(x)=q(x)*g(x)+r(x)\) 可知：

\[f(\frac{1}{x})=q(\frac{1}{x})*g(\frac{1}{x})+r(\frac{1}{x})

\]

两边同时乘上 \(x^n\) 得：

\[f^R(x)=q^R(x)*g^R(x)+x^{n-m+1}r^R(x)

\]

注意到我们只需要求出 \(q(x)\) 和 \(r(x)\) 中的任意一个就能求出另一个。

我们发现上式中 \(\deg q(x)=n-m\) ，而 \(r^R(x)\) 的系数为 \(x^{n-m+1}\) 。

因此我们把式子两边模上 \(x^{n-m+1}\) :

\[f^R(x)\equiv g^R(x)q^R(x) \pmod {x^{n-m+1}}

\]

\[f^R(x)(g^R(x))^{-1}\equiv q^R(x) \pmod {x^{n-m+1}}

\]

于是多项式求逆解决 \(q(x)\) ，再代入原式解出 \(r(x)\) 即可。

注意要先把 \(g^R(x)\) 的次数设为 \(n-m\) 再求逆。

时间复杂度 \(O(n\log n)\) 。

代码

typedef pair pvv;//除法返回答案时要用pair，你也可以设几个全局变量

#define mk make_pair

#define F first

#define S second

void R(vi &a){reverse(a.begin(),a.end());}//简单的reverse

pvv operator /(vi a,vi b)

{

vi q,r;

q=b;R(q);cut(q,a.size()-b.size()+1);//把q设为bR，把q的次数设为n-m（size设为n-m+1）

r=a;R(r);//r作为临时变量存储aR

q=inv(q);q*=r;//计算qR

cut(q,a.size()-b.size()+1);R(q);//将qR多余位数截掉后还原q

r=q;r*=b;r-=a;--r;cut(r,b.size()-1);//计算r

return mk(q,r);

}

多项式开方

给定一个 \(n\) 次多项式 \(f(x)\) ，请求出多项式 \(g(x)\) ，满足 \(g^2(x)=f(x)\) 。

设 \(f(x)\) 的系数为 \(a_0,a_1,\cdots,a_n\) ， \(g(x)\) 的系数为 \(b_0,b_1,\cdots,b_n\) 。

首先我们先求出 \(g(x)\) 的常数项 \(b_0\) 。

代入 \(x=0\) 可得： \(b_0^2=a_0\) ，因此如果 \(a_0\) 在模意义下没有平方根就无解。

于是我们取 \(b_0=\sqrt{a_0}\) 。

（由于 \(a_0\) 的平方根可能不止一个，我们随便取一个值作为 \(b_0\) ；选择不同的 \(b_0\) 会求出不同的 \(b(x)\) 。）

以下我们先假定 \(a_0=1\) ，因为这样讨论起来比较方便。

老规矩，我们继续把 \(n\) 扩大到 \(2^k-1\) 的形式。

假设我们求出了 \(f(x)\) 在模 \(x^\frac{n}{2}\) 意义下的平方根 \(g_0(x)\) 。

记 \(m=x^n,m'=x^{\frac n2}\) 。则有：

\[f(x)-g_0^2(x)\equiv 0 \pmod {m'}

\]

两边同时平方，得：

\[\begin{aligned}

(f(x)-g_0^2(x))^2&\equiv 0 \pmod m\\

(f(x)+g_0^2(x))^2&\equiv 4f(x)g_0^2(x) \pmod m

\end{aligned}

\]

将 \(4g_0^2(x)\) 移到左边得：

\[\left(\frac{f(x)+g_0^2(x)}{2g_0(x)}\right)^2\equiv f(x) \pmod m

\]

两边同时开方：

\[\begin{aligned}

\frac{f(x)+g_0^2(x)}{2g_0(x)}&\equiv g(x) \pmod m\\

2^{-1}(f(x)g_0^{-1}(x)+g_0(x))&\equiv g(x) \pmod m

\end{aligned}

\]

倍增计算即可。

时间复杂度： \(T(n)=T(\frac n2)+O(n \log n)\) ，即 \(T(n)=O(n \log n)\) 。

代码

vi sqrt(vi &a)

{

int x=1;

vi ans={1},b,now;

while(1ull*x

{

x<<=1;

cut(ans,x);//这里求逆的时候要

copy(a,x,b);

now=ans;ans=inv(ans);ans*=b;ans+=now;ans*=i2;

cut(ans,x);

}

cut(ans,a.size());

return ans;

}

（为了直观，我写的常数比较大，建议读者写的时候不要每次都做一遍cut和copy，建议多传几个参数）

然后我们解决 \(a_0\) 为任意数时的问题。（一般不会这么毒瘤）

若 \(a_0\ne 0\) ，我们用二次剩余即可。

否则，我们考虑原式中因式 \(x\) 的次数 \(k\) 。（即：找到一个最大的 \(k\) ，使得 \(a_0,a_1,a_2,\cdots,a_k=0\) ）

如果 \(k\) 是奇数，则无解；如果 \(k\) 是偶数，输出 \(\sqrt{\frac{a}{x^k}}*x^\frac k2\) 即可。

分治FFT

经典例题：给你一个数组 \(g_i\) ，求数组 \(f_i\) ，使得 \(f_i=\sum_{j=1}^if_{i-j}g_j\) 且 \(f_0=1\)。

我们先来对比一下普通卷积和分治FFT的式子：

普通卷积：\(h_i=\sum_{j=1}^if_{i-j}g_j\)

分治FFT：\(f_i=\sum_{j=1}^if_{i-j}g_j\)

我们发现：分治FFT实际上和普通卷积比较像，就是输出的数组就是它本身。

但这个性质好像没啥卵用

我们发现，对于每个区间 \([l,r]\) ，记 \(mid=\lfloor\frac{l+r}2\rfloor\) ，考虑一种类似于CDQ分治的算法：

递归 \([l,mid]\) 算出递归区间内的答案。

把 \([l,mid]\) 的贡献加到 \([mid+1,r]\) 上。

递归 \([mid+1,r]\) 。

其中第二步我们再详细讲一下。

假设当前区间为 \([0,5]\) ，则我们需要

\[\begin{aligned}

f_3&+=f_2*g_1+f_1*g_2+f_0*g_3\\

f_4&+=f_2*g_2+f_1*g_3+f_0*g_4\\

f_5&+=f_2*g_3+f_1*g_4+f_0*g_5

\end{aligned}

\]

于是，\(f_n+=\sum_{i=l}^{mid}f_ig_{n-i}=[x^n]f'(x)g'(x)\) ，其中 \(f'(x)\) 代表的是只考虑 \([l,mid]\) 这一段时的 \(f(x)\) ， \(g'(x)\) 代表只考虑 \([1,r-l]\) 时的 \(g(x)\) 。

卷积即可。

时间复杂度： \(T(n)=2T(\frac n2)+O(n \log n)\) ，即 \(T(n)=O(n \log^2 n)\) 。

代码

void copy(vi &a,int l,int r,vi &b)//把a的第l位到第r位复制给b

{

b.clear();

for(int i=l;i<=r;++i)b.push_back(a[i]);

}

void div_FFT(vi &a,vi &b,int l=-1,int r=-1)

{

if(l==-1&&r==-1)

{

cut(a,b.size());

l=0;r=b.size()-1;//初始化

}

if(l==r)return;

int mid=(l+r)>>1;

div_FFT(a,b,l,mid);//分治左区间

vi c,d;

copy(b,1,r-l,c);copy(a,l,mid,d);

c*=d;

for(int i=mid+1;i<=r;++i)a[i]=(a[i]+c[i-l-1])%mod;//卷积处理对右半的贡献

div_FFT(a,b,mid+1,r);//分治右区间

}

但是先别走！这道题也可以用多项式求逆解决！！！

我们发现

\[f(x)g(x) \equiv \sum_{i=0}^{n-1} x^i \sum_{j=0}^if_{i-j}g_j \pmod{x^n}

\]

由于 \(g_0=0\)，所以

\[f(x)g(x) \equiv \sum_{i=1}^{n-1} x^i \sum_{j=1}^if_{i-j}g_j \pmod{x^n}

\]

并且

\[\begin{aligned}

f(x)&\equiv\sum_{i=0}^{n-1} x^if_i\\

&\equiv x^0f_0+\sum_{i=1}^{n-1} x^if_i\\

&\equiv f_0+\sum_{i=0}^{n-1} x^i \sum_{j=1}^if_{i-j}g_j \pmod{x^n}

\end{aligned}

\]

则有 \(f(x)g(x)\equiv f(x)-f_0\) ，于是就有 \(f(x)(1-g(x))\equiv f_0\)，从而 \(f(x)\equiv f_0(1-g(x))^{-1}\)，多项式求逆即可。

时间复杂度 \(O(n\log n)\)。

多项式 \(\ln\)

各位都知道整数的 \(\ln\) 和 \(\exp\) 吧。（ \(exp(x)\) 实际上就是 \(e^x\)）

后面的证明需要依赖微积分，不会的同学请自行跳过。

实际上， \(\ln\) 和 \(\exp\) 的完整定义在这：

\[\ln f(x)=-\sum_{i=1}^\infty \frac{(1-f(x))^i}{i}\\

\]

\[\exp f(x)=\sum_{i=0}^\infty \frac{f(x)^i}{i!}

\]

但是，直接用定义式时间复杂度会炸。

首先，对于多项式 \(f(x)\)，若 \(\ln(x)\) 存在，则有 \([x_0]f(x)=1\) 。（否则常数项不收敛）

我们知道 \(\ln x\) 的导数为 \(\frac 1x\) ，不妨设 \(\ln f(x)=g(x)\) ，两边同时对 \(x\) 求导可得：

\[\frac 1{f(x)}\cdot f'(x)=g'(x)

\]

\[\frac {f'(x)}{f(x)}=g'(x)

\]

再同时积分： $$g(x)=\int\left(\frac {f'(x)}{f(x)}\right)$$

因此我们只需要多项式求导、积分、求逆和乘积就可以了。

时间复杂度 \(O(n\log n)\)。

代码

int invv[2100000];

void init(int n=2099999)

{

invv[1]=1;

for(int i=2;i<=n;++i)invv[i]=(mod-mod/i)*invv[mod%i]%mod;

}//线性求逆元

void D(vi &a){for(int i=1;1ull*i

void I(vi &a){a.push_back(0);for(int i=a.size()-1;i>=1;--i)a[i]=invv[i]*a[i-1]%mod;a[0]=0;}//积分

void ln(vi &a)

{

vi b=inv(a);

D(a);

a*=b;I(a);

cut(a,b.size());

}

多项式 \(\exp\)

首先，对于多项式 \(f(x)\)，若 \(\exp(x)\) 存在，则有 \([x_0]f(x)=0\) 。（否则常数项不收敛）

我们知道 \(\exp x\) 的导数就是它本身，不妨设 \(\exp f(x)=g(x)\) 。

两边求导得： \(g(x)*f'(x)=g'(x)\)

提取系数后就变成了 \(g'_n=\sum_{i=0}^{n}g_{n-i}f'_i\)

由多项式的求导，有 \((n+1)g_{n+1}=\sum_{i=0}^{n}g_{n-i}f_{i+1}(i+1)\)

即 \(ng_n=\sum_{i=1}^{n}(i\cdot g_{n-i}f_i)\)，\(g_n=\frac 1n\sum_{i=1}^{n}(i\cdot g_{n-i}f_i)\)。

并且 \(g_0=1\) ，分治FFT即可。

虽然和分治FFT的模板不同，但实际上差不多，我们把 \(i\) 乘到 \(f_i\) 上，再在 \(l=r\) 处把常数 \(\frac 1n\) 乘上。

其实，把上述式子再变换也可以求出 \(\ln\) 的递推形式，这里不再阐述，感兴趣的同学可以自己推一下。

时间复杂度 \(O(n\log^2 n)\)。

代码

void div_FFT(vi &a,vi &b,bool isexp=0,int l=-1,int r=-1)//对之前的分治FFT做了点手脚

{

if(l==-1&&r==-1)

{

if(a.empty())a.push_back(1);

cut(a,b.size());

l=0;r=b.size()-1;

}

if(l==r){if(isexp&&l)a[l]=(a[l]*invv[l])%mod;return;}//边界条件略改一下

int mid=(l+r)>>1;

div_FFT(a,b,isexp,l,mid);

vi c,d;

copy(b,1,r-l,c);copy(a,l,mid,d);

c*=d;

for(int i=mid+1;i<=r;++i)a[i]=(a[i]+c[i-l-1])%mod;

div_FFT(a,b,isexp,mid+1,r);

}

vi exp(vi &a)

{

for(int i=0;1ull*i

vi b={};

div_FFT(b,a,1);//1代表分治FFT最终的系数要除以下标

return b;

}

多项式快速幂

我们考虑封装多项式乘法，再倍增即可。

好吧，正解是这样的：

我们先保证 \(a_0=1\) 。

众所周知，多项式的 \(\ln\) 有个性质： \(\ln f^k(x)=k\ln f(x)\) 。

于是我们输出 \(\exp(k\cdot\ln f(x))\) 即可。

时间复杂度 \(O(n\log n)\) 。

但这还没完！

我们发现洛谷上的模板题是这样的：

\(0

那怎么办呢？还要写个高精？

肯定不可能，我们需要知道一个性质：

\[f(x)^k\equiv f(x)^{k\bmod p}\pmod{x^n}

\]

这个式子在 \(n

证明

我们先证明一下这样一个结论：\(f(x)^p\equiv f(x^p)\)。

显然有个式子叫做 \((a+b)^p\equiv a^p+b^p\) 。（组合数拆分系数）

不妨设上述结论对所有 \(k-1\) 次多项式成立。

设 \(f(x)\) 的次数为 \(k\) ，并且 \(f(x)=a_kx^k+g(x)\)，

则 \(f(x)^p\equiv g(x)^p+(a_kx^k)^p\equiv g(x^p)+a_k^px^{pk}\equiv g(x^p)+a_kx^{pk}=f(x^p)\)。

有了这个式子以后，由于 \(n

有些同学可能到了这里还有一个疑问：为什么整数快速幂的时候指数要模 \(p-1\) ，但多项式快速幂就变成模 \(p\) 了呢？

换句话说，如果把 \(f(x)\) 带入具体的值，那么模数应该模上 \(p-1\) ，而不是 \(p\) 。

其实，这两个式子还是有一点差别的。

快速幂的式子是：

\[f(x)^k\equiv f(x)^{k\bmod (p-1)} \pmod p

\]

而在 \(n

\[f(x)^k\bmod {x^n}\equiv f(x)^{k\bmod p}\bmod x^n

\]

也就是说一个式子算完了，而另一个式子只计算了前 \(n\) 项，因此它们推导的方式也不同。

一定要注意多项式快速幂的证明要依赖 \(n

有了这个性质，我们就能很轻松地解决刚开始的问题了，只需要在读入的时候将 \(k\) 模上 \(p\) 即可，根本不需要高精。

啥？怎么在读入的时候将 \(k\) 模上 \(p\) ？你想想快读是怎么读入的就知道了。

啥？还不会？快读每次乘 \(10\) 都取模啊！

代码

inline int read(int mod=0)//修改后的快读，加入了取模

{

char ch=getchar();int x=0,r=1;

while(ch<'0'||ch>'9'){if(ch=='-')r=0;ch=getchar();}

while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';if(mod)x%=mod;ch=getchar();}//这里有改动

return r?x:-x;

}

void mul(vi &a,int k)

{

ln(a);

a*=k;

a=exp(a);

}

接下来我们考虑 \(a_0\ne1\) 的情况。

当 \(a_0\ne0\) 时，我们先把整个式子除以 \(a_0\) ，然后再进行快速幂，最后再把整个式子乘上 \(a_0^k\) 即可。

对于 \(a_0=0\) 的情况，我们就把原式中的因式 \(x^m\) 提出来，转化成上一种情况求解。

注意这里要同时保留 \(k\bmod p\) 和 \(k\bmod (p-1)\) 的值，因此不能直接快读，建议把结果存到字符串里。

多项式全家桶代码

由于之前的代码段太多了，函数也比较多，所以这里整合了一下。

建议大家先把所有的函数写完再一起卡常，防止自己看不懂

全家桶代码

#include

using namespace std;

#define int long long

typedef vector vi;

typedef pair pvv;

#define mk make_pair

#define F first

#define S second

const int mod=998244353,g=3,ig=332748118,i2=499122177;

inline int read(int mod=0)

{

char ch=getchar();int x=0,r=1;

while(ch<'0'||ch>'9'){if(ch=='-')r=0;ch=getchar();}

while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';if(mod)x%=mod;ch=getchar();}

return r?x:-x;

}

int mul(int x,int b)

{

int ans=1;

while(b)

{

if(b&1)ans=(ans*x)%mod;

x=(x*x)%mod;

b>>=1;

}

return ans;

}

int invv[2100000];

void init(int n=2099999)

{

invv[1]=1;

for(int i=2;i<=n;++i)invv[i]=(mod-mod/i)*invv[mod%i]%mod;

}

int rev[2100000];

void read(int n,vi &a,int l=0){for(int i=0;i<=l+n-1;++i)a.push_back(i

void print(vi &a){for(int i:a)printf("%lld ",i);puts("");}

void D(vi &a){for(int i=1;1ull*i

void I(vi &a){a.push_back(0);for(int i=a.size()-1;i>=1;--i)a[i]=invv[i]*a[i-1]%mod;a[0]=0;}

void NTT(vi &x,int limit,bool type=1)

{

for(int i=0;i

for(int len=1;len

{

int Wn=mul(type?g:ig,(mod-1)/(2*len));

for(int i=0;i

{

int W=1,X,Y;

for(int j=i;j

{

X=x[j];Y=W*x[j+len]%mod;

x[j]=(X+Y)%mod;x[j+len]=(X-Y+mod)%mod;

W=(W*Wn)%mod;

}

}

}

if(!type)

{

int invs=mul(limit,mod-2);

for(int i=0;i

}

}

int limit;

void cut(vi &a,int x){while(a.size()>1ull*x)a.pop_back();while(a.size()<1ull*x)a.push_back(0);}

void copy(vi &a,int x,vi &b)

{

if(1ull*x<=a.size()){b.clear();for(int i=0;i

else {b=a;while(b.size()<1ull*x)b.push_back(0);}

}

void copy(vi &a,int l,int r,vi &b)

{

b.clear();

for(int i=l;i<=r;++i)b.push_back(a[i]);

}

void operator +=(vi &a,int x){a.front()=(a.front()+x)%mod;}

void operator +=(vi &a,vi &b)

{

while(a.size()

for(int i=0;1ull*i

}

void operator -=(vi &a,int x){a.front()=(a.front()-x+mod)%mod;}

void operator -=(vi &a,vi &b)

{

while(a.size()

for(int i=0;1ull*i

}

void operator *=(vi &a,int x){for(int i=0;1ull*i

void operator *=(vi &a,vi b)

{

int x=a.size()+b.size()-1;

limit=1;

while(1ull*limit<=a.size()-1+b.size()-1)limit<<=1;

for(int i=0;i

while(a.size()<1ull*limit)a.push_back(0);

while(b.size()<1ull*limit)b.push_back(0);

NTT(a,limit);NTT(b,limit);

for(int i=0;i

NTT(a,limit,0);

cut(a,x);

}

vi inv(vi &a)

{

int x=1;

vi ans={mul(a.front(),mod-2)},b,now;

while(1ull*x

{

x<<=1;

copy(a,x,b);

now=ans;ans*=b;ans*=-1;ans+=2;ans*=now;

cut(ans,x);

}

cut(ans,a.size());

return ans;

}

void R(vi &a){reverse(a.begin(),a.end());}

pvv operator /(vi a,vi b)

{

vi q,r;

q=b;R(q);cut(q,a.size()-b.size()+1);

r=a;R(r);

q=inv(q);q*=r;

cut(q,a.size()-b.size()+1);R(q);

r=q;r*=b;r-=a;r*=-1;cut(r,b.size()-1);

return mk(q,r);

}

vi sqrt(vi &a)

{

int x=1;

vi ans={1},b,now;

while(1ull*x

{

x<<=1;

cut(ans,x);copy(a,x,b);

now=ans;ans=inv(ans);ans*=b;ans+=now;ans*=i2;

cut(ans,x);

}

cut(ans,a.size());

return ans;

}

void div_FFT(vi &a,vi &b,bool isexp=0,int l=-1,int r=-1)

{

if(l==-1&&r==-1)

{

if(a.empty())a.push_back(1);

cut(a,b.size());

l=0;r=b.size()-1;

}

if(l==r){if(isexp&&l)a[l]=(a[l]*invv[l])%mod;return;}

int mid=(l+r)>>1;

div_FFT(a,b,isexp,l,mid);

vi c,d;

copy(b,1,r-l,c);copy(a,l,mid,d);

c*=d;

for(int i=mid+1;i<=r;++i)a[i]=(a[i]+c[i-l-1])%mod;

div_FFT(a,b,isexp,mid+1,r);

}

void ln(vi &a)

{

vi b=inv(a);

D(a);

a*=b;I(a);

cut(a,b.size());

}

vi exp(vi &a)

{

for(int i=0;1ull*i

vi b={};

div_FFT(b,a,1);

return b;

}

void mul(vi &a,int k)

{

ln(a);

a*=k;

a=exp(a);

}

signed main()

{

return 0;

}

多项式牛顿迭代

开个坑，以后再写。