编写之初

由于利用留数定理分解分式多项式的计算麻烦，所以决定用python做一个利用留数定理分解分式多项式程序，实现只要输入多项式就可以得到各种中间参数和最终拆分结果的目的。从本程序可以得到：分解后每项多项式分子值、计算分解后每项多项式分子值的过程展示、最终分解结果展示。

算数实现

利用留数定理分解分式多项式计算过程以下图为例所示

对于算数运算形象的理解：以该式为例，分解后各项分母为原多项式分母中的各项（s、(s+2)、(s+3)²）加上这些项中高次幂项去掉次幂为分母的项（(s+3)）;分解后各项分子为原多项式分母依次去掉各项后求极限所得值（-4、1/3、3）以及n(n>1)次幂项去掉该项后的n阶微分式求极限所得值（-10/3）。也依据这个思路编写程序。

依赖的包

1.sympy模块
sympy模块，可以进行符号计算，可以定义符号变量，进行代数运算，以及微分运算、积分运算等。

pip install sympy

2.re正则化模块
re模块是python独有的匹配字符串的模块，该模块中提供的很多功能是基于正则表达式实现的。

pip install re

完整代码及注释

import sympy as sp
import re

s = sp.symbols('s')     # 创建符号变量

print('以 6*(s+1)/(s*(s+2)*(s+3)**2) 形式为例')
Y = input('请输入待分解多项式：')   # 复制 例 6*(s+1)/(s*(s+2)*(s+3)**2)

# denominator = y.split('/')[0]   # 分子项
molecule = Y.split('/')[1]      # 分母项

moleculelist = molecule.strip('()').split('*')  # 分母每项
# print(moleculelist)

reg = re.compile(r"(?<=\*)\d+")     # 获取分式中的幂数
match=reg.search(Y)
pow = match.group(0)

molecule_number = 0     # 初始化分母项索引
# 计算拆分后各项值
while molecule_number < len(moleculelist):
    try:
        if moleculelist[molecule_number] != '':
            if moleculelist[molecule_number + 1] != '':
                y = Y.replace(moleculelist[molecule_number]+'*', '')    # 得到用来计算 该项分子值 的多项式
                if moleculelist[molecule_number] == 's':
                    value = sp.limit(y, s, 0)               # 求极限计算 该项分子值
                    print('\n计算分母为' + str(moleculelist[molecule_number]) +'项的分子：limit——>0 ', y)
                    print('计算得；', value)
                else:
                    limit_value = -int(re.findall('(\d+)', moleculelist[molecule_number])[0])   # 获取求极限时的 极限参数
                    value = sp.limit(y, s, limit_value)         # 求极限计算 该项分子值
                    print('计算分母为' + str(moleculelist[molecule_number]) +'项的分子：limit——>' + str(limit_value), y)
                    print('计算得：', value)
            else:
                y = Y.replace('*' + moleculelist[molecule_number] + '**' + pow, '')     # 得到用来计算 该项分子值 的多项式
                limit_value = -int(re.findall('(\d+)', moleculelist[molecule_number])[0])       # 获取求极限时的 极限参数
                value = sp.limit(y, s, limit_value)         # 求极限计算 该项分子值
                print('计算分母为' + str(moleculelist[molecule_number]) + '**' + pow + '项的分子：limit——>' + str(limit_value), y)
                print('计算得：', value)

                y_diff = sp.together(sp.diff(y, s, int(pow)-1))      # 求解 除n次幂项外分式 n阶导数
                limit_value = -int(re.findall('(\d+)', moleculelist[molecule_number])[0])       # 获取求极限时的 极限参数
                value = sp.limit(y_diff, s, limit_value)        # 求极限计算 该项分子值
                print('计算分母为' + str(moleculelist[molecule_number]) + '项的'+ str(int(pow)-1) + '阶导数的分子：limit——>' + str(limit_value), y_diff)
                print('计算得：', value)

    except: pass
    molecule_number += 1

# print(molecule+'\n'+molecule)
print('分解结果为：',sp.apart(Y, s))      # 拆分结果

运行结果

运行代码，首先以例示形式输入待分解多项式，得到分解过程及结果

感悟与不足

本来打算纯手写一个留数定理分解分式多项式的程序，但是由于其中涉及到极限与微分的计算，而这两个模块编写起来又要一定时间，于是借助了sympy模块辅助编程，但是我意外地发现，sympy模块居然自带多项式分解函数 sympy.apart()因此我也在呈现分解结果时用到了这个函数，免去了字符串处理当中的一系列麻烦。因此如果大家对多项式分解的过程不在意，只想得到分解结果，只需执行sympy.apart()函数即可，至于它其中的运算原理我也不得而知。

关于多项式

一个形如$f(x)={\sum }_{i=0}^n{a}_i{x}^i(a_n\ne 0)$的式子被称为一个$n$次多项式。例如，$3x^3+5x^2+x+2021$是一个$3$次多项式。

一个多项式有多种表示方法，其中最经典的是系数表示法。形式化地说，${\sum }_{i=0}^n{a}_i{n}^i$被表示为$(a_0,a_1,a_2\cdots a_n)$。一些多项式的题目中往往会直接或间接地给出这些系数。下面的题目中，“给定一个多项式”默认为“给出此多项式的各项系数”。

多项式还有点值表示法。假设这个多项式为$f$，那么在$f$的图像上任意取出$n+1$个点即可确定这个多项式的各项系数。例如，$3x+2$的一种点值表示为$(0,2)(1,5)$。点值表示法在多项式乘法中与多项式求逆中有很大的用处。

---

在多项式的学习中，可能会遇到各种各样的数学大坑，这一部分也是学习中的重难点，这也导致了笔者曾经几次三番放弃学习多项式。考虑到我的较差体验，笔者将会耐心地为大家去补充这些知识，相信你能够迅速理解这些妙妙的东西。

鸣谢: gh巨佬(他教会了我这个东西，还为本文贡献了许多图片)

多项式乘法

Description

给定一个$n$次多项式$A$和$m$次多项式$B$，请求出$A$与$B$的卷积(即$A\times B$的各项系数)。

$n,m\le 10^6$

Sample

Input
1 1
1 1
1 1

Output
1 2 1

Explanation
$(x+1)(x+1)=x^2+2x+1$

复数的表示、基本定义与一些性质

在初中阶段，我们已经学过了实数，它的特点是能在一条数轴上表示出来，比如$\sqrt{3},2021,\frac{998244353}{10^9+7}$。

复数是形如$a+bi$的数，其中$i=\sqrt{-1}$，它也被称为虚数单位。不难发现，当$b=0$时该数为实数，否则为虚数。

显然，虚数是不能使用一条数轴表示出来的。

那我们采用两条数轴呢？考虑一个“平面直角坐标系”: 横轴表示实部(即$a+bi$中的$a$)，纵轴表示虚部(即$a+bi$中的$bi$)，我们可以将$a+bi$表示在$(a,b)$处！如下图:

我们成功地将复数表示在了“平面直角坐标系”上。现在，我们再赋予复数几个定义：

①幅角$\phi$(虽然这是欧拉函数, 但是请您将就着看吧/kk): $(a,b)$与原点连线和$x$轴所形成的角的大小(注意区分夹角与角度)。
②模长$d$: $(a,b)$到原点的距离。

如下图:

好了，我们的定义到此为止。剩下的工作交给欧拉。

Euler: “$e^{i\phi }=\cos (\phi )+i\sin (\phi )$”

式子中的$e$是自然对数（即${\lim }_{x\to \infty }(x+\frac{1}{x}{)}^x$），$i$是虚数单位。

这是十分有用的一个公式。由于笔者太菜不会证明所以证明省略 结合我们上面对复数的定义与在坐标系上的表示，不难发现$e^{i\phi }$对应的就是一个模长为$1$，幅角为$\phi$的点！

为什么呢？考虑$\cos (\phi )+i\sin (\phi )$的坐标为$(\cos (\phi ),\sin (\phi ))$。根据勾股定理，此时的模长为$\sqrt{co{s}^2(\phi )+si{n}^2(\phi )}=1$。根据$sin$的定义，不难得到幅角为$\phi$。

更一般性地，一个模长为$a$，幅角为$\phi$的复数总可以用$a\times e^{i\phi }$来表达。为了方便叙述，我们将这个复数简记为$(a,\phi )$。

考虑两个复数相乘：$(a,{\phi }_1)\times (b,{\phi }_2)$，根据欧拉定理，这个式子等价于

$(a\times e^{i{\phi }_1})(b\times e^{i{\phi }_2})$

$=(ab)\times (e^{i({\phi }_1+{\phi }_2)})$

综上所述，复数相乘，模长相乘且幅角相加。

FFT预备知识

接下来的时间交给傅里叶。

他让我们画了一个以原点为圆心且半径为$1$的圆(是不是很可爱呢 )。然后我们把它$n(n=2^k,k\in N)$等分了：

于是我们在圆上得到了$8$个表示复数的点。从幅角为$0°$的点开始逆时针标号，分别标号为$0,1,2\cdots ,n-1,n,n+1\cdots$其中，标号为$k$的点简记为$w_n^k$。$w_n^k$被称为一个单位根。

单位根是有许多性质的。下面列出了部分重要性质：

①$n$个单位根两两不同；
②$w_n^k=w_{2n}^{2k}$
③$(w_n^k{)}^n=1$
④$w_n^k=-w_n^{k-\frac{n}{2}}$
⑤$w_n^a\times w_n^b=w_n^{a+b}$

证明十分简单，只需要瞅一瞅那个被$n$等分的圆即可。如果觉得不够严谨，也可以结合复数相乘“模长相乘，幅角相加”的性质来详细证明。说白了就是我懒不想证了

这些性质有什么用呢？马上你就知道了。

FFT

假设我们想要求出两个多项式$f,g$的乘积。

首先，一种方法是枚举$f,g$中的两项并向答案对应的项做贡献。这样子的时间复杂度是$O(nm)$的。

傅里叶又开始吊打我们了……祂提出了FFT的详细流程：

①带入$n$个单位根作为$x$，分别求出对应的$y=f(x)$(即转换为了点值表示法)；
②带入$n$个单位根作为$x$，分别求出对应的$y=g(x)$；
③将对应的点值相乘，即将$(x_1,y_1)(x_1,y_2)$转换为$(x_1,y_1{y}_2)$；
④将点值表示法还原为系数表示法，得到答案。

其中，①被称为DFT，④被称为IDFT，整个流程被称为FFT，即“傅里叶快速变换”。

考虑①②是基本相同的，③直接$O(n)$相乘，④是①、②的逆变换。现在，关键在于求出①与④的做法。

DFT

现在我们想要快速求出${\sum }_{i=0}^{n-1}{a}_i{w}_n^i$。考虑将它按照$i$的奇偶性拆开：

${\sum }_{i=0}^{n-1}{a}_i{w}_n^i$
$={\sum }_{i=0}^{\lfloor \frac{n}{2}\rfloor -1}{a}_{2i}{w}_n^{2ik}+{\sum }_{i=0}^{\lfloor \frac{n}{2}\rfloor -1}{a}_{2i+1}{w}_n^{(2i+1)k}$
$={\sum }_{i=0}^{\lfloor \frac{n}{2}\rfloor -1}{a}_{2i}{w}_n^{2ik}+w_n^k{\sum }_{i=0}^{\lfloor \frac{n}{2}\rfloor -1}{a}_{2i+1}{w}_n^{2ik}$

根据单位根的性质②$w_n^k=w_{2n}^{2k}$，不难得到

$={\sum }_{i=0}^{\lfloor \frac{n}{2}\rfloor -1}{a}_{2i}{w}_n^{2ik}+w_n^k{\sum }_{i=0}^{\lfloor \frac{n}{2}\rfloor -1}{a}_{2i+1}{w}_n^{2ik}$
$={\sum }_{i=0}^{\lfloor \frac{n}{2}\rfloor -1}{a}_{2i}{w}_{\frac{n}{2}}^{ik}+w_n^k{\sum }_{i=0}^{\lfloor \frac{n}{2}\rfloor -1}{a}_{2i+1}{w}_{\frac{n}{2}}^{ik}$

然后分治下去就做完了？

并不是。考虑直接分治的话，每次都要带入$w_n^0,w_n^1\dots \dots w_n^{n-1}$共$n$个单位根。虽然系数每次减半，但是带入的值并没有减半，这导致我们的时间复杂度依然不对。

考虑对于$k\le \frac{n}{2}$的时候按照上述方式去分治。当$\frac{n}{2}+1\le k$时，我们令$K=k-\frac{n}{2}$，将其带入得

${\sum }_{i=0}^{n-1}{a}_i{w}_n^{(K+\frac{n}{2})i}$
$={\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i}{w}_n^{(K+\frac{n}{2})2i}+{\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{w}_n^{(K+\frac{n}{2})(2i+1)}$

考虑有
(1)$w_n^{(K+\frac{n}{2})2i}$
$=w_n^{2iK}\times w_n^{ni}=w_{\frac{n}{2}}^{iK}$

(2)$w_n^{(K+\frac{n}{2})(2i+1)}$
$=w_n^{(2iK+K+ni+\frac{n}{2})}$
$=w_n^{2iK}\times w_n^K\times w_n^{ni}\times w_n^{\frac{n}{2}}$
$=-w_{\frac{n}{2}}^{iK}\times w_n^K$

往回带入得

${\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i}{w}_{\frac{n}{2}}^{Ki}-w_n^K{\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{w}_{\frac{n}{2}}^{Ki}$

至此，FFT的推导到此结束。我们总结一波。

用$A$表示一个用系数表示的多项式，$A_1$表示$A$的偶数项的系数形成的多项式，$A_2$表示$A$的奇数项的系数形成的多项式。则

$A(w_n^k)=A_1(w_{\frac{n}{2}}^k)+w_n^k{A}_2(w_{\frac{n}{2}}^k)$
$A(w_n^{k+\frac{n}{2}})=A_1(w_{\frac{n}{2}}^k)-w_n^k{A}_2(w_{\frac{n}{2}}^k)$

此时我们分治下去的时间复杂度是什么呢？

假设对于一个项数为$n$的多项式，做一遍DFT的时间复杂度为$T(n)$。每次拆开奇偶位并进行分治，这一部分的时间复杂度是$2T(\frac{n}{2})$；然后还要对于$w_n^0,w_n^1\dots \dots w_n^{n-1}$分别求出它们带进去当前多项式后得到的值，由于每次调用都是$O(1)$的，所以这一部分的时间复杂度是$O(n)$的。

即$T(n)=2T(\frac{n}{2})+O(n)$。通过主定理不难得到$T(n)=O(n\log n)$。

IDFT

学会了DFT你是不是很开心？可是如果不会IDFT，一切都是白搭/kk

不难发现，这一部分是一个解方程组的问题。

盗个图：

我们想要求出的就是$a_0,a_1\cdots a_{n-1}$。

写成矩阵的形式：

令上图中从左到右的三个矩阵分别是$A,X,B$。

考虑另一个矩阵$T$，其中$T_{ij}=w_n^{-ij}$。

令$F=A\times T$。考虑$F$长啥样，我们分类讨论如下：

①$i=j$，此时$F_{i,j}={\sum }_{k=1}^n{w}_n^{ik}{w}_n^{-jk}={\sum }_{k=1}^n{w}_n^0=n$。

②$i\ne j$

${\sum }_{k=0}^{n-1}{w}_n^{-ik}{w}_n^{kj}$

${\sum }_{k=0}^{n-1}{w}_n^{k(j-i)}$

这不就是一个等比数列求和的形式吗？

${\sum }_{k=0}^{n-1}{w}_n^{k(j-i)}$

$=\frac{1-w_n^{n(j-i)}}{1-w_n^{j-i}}$

由于$w_n^{n(j-i)}=1$，所以上面这个式子为$0$。

---

于是$F$是一个单位矩阵的$n$倍！

于是$A$的逆矩阵就是$T$除以$n$！

更形式化地说，

综上所述，我们只需要将$w_n^{-0},w_n^{-1}\cdots w_n^{-(n-1)}$逐一带入，用类似$DFT$的方法与各项系数为$A(w_n^0),A(w_n^1)\cdots A(w_n^{n-1})$的多项式相乘，最终再将各项系数除以$n$即可得到最终的答案。

显然这一部分的时间复杂度仍然是$O(n\log n)$的。

FFT之再优化

到这一步，我们似乎已经得到了$O(n\log n)$的解法。但是很可惜，这种做法的常数非常之大。gh尝曰: “此法尚不如暴力也。”

所以我们需要对这种做法的常数进行一定的优化。

我们观察一下做FFT的过程：

操作的位置的原序列与后序列如下：

发现了什么？将原序列的下标做二进制反转后就是后序列！

所以我们可以采用数位$dp$处理出每一个在$[0,n-1]$范围内的数反转后的值。每次在做DFT/IDFT前进行反转即可。

Attention

①要多次做复数相加，相减和相乘的运算，可以将一个复数拆成实部与虚部并放在一个结构体里面，然后重载运算符(比如$+,-,\times$)；

②单位根$w_n^k=cos(\frac{2\pi }{n}k)+isin(\frac{2\pi }{n}k)$——请注意C++里面的三角函数采用了弧度制；

③FFT有精度误差，最后输出的时候要精确到整数部分输出；

④$n,m$较大请使用快读；

⑤答案的最高项的次数可能会达到$n+m$，所以我们必须让$A,B$与答案的项数为一个大于$n+m$的$2$的次方数。

其他详见有注释的代码。

Code: FFT

#include <bits/stdc++.h>
#define int long long
#define double long double
using namespace std;
const int maxl=5000005;
const double PI=acos(-1.00);

int read(){//快读
	int s=0,w=1;char ch=getchar();
	while (ch<'0'||ch>'9'){if (ch=='-')  w=-w;ch=getchar();}
	while (ch>='0'&&ch<='9'){s=(s<<1)+(s<<3)+(ch^'0');ch=getchar();}
	return s*w;
}
int n,m,len=1,cnt=0;
int r[maxl];

struct Complex{
	double x,y;
	Complex(double a=0,double b=0){x=a,y=b;}
	Complex operator + (const Complex a){return Complex(a.x+x,a.y+y);}
	Complex operator - (const Complex a){return Complex(x-a.x,y-a.y);}
	Complex operator * (const Complex a){return Complex(x*a.x-y*a.y,x*a.y+y*a.x);}
	//重载运算符
}a[maxl],b[maxl];

void FFT(Complex *a,int len,int flag){//flag=1表示在做DFT，否则在做IDFT
	for (int i=0;i<len;i++){
		if (i<r[i])  swap(a[i],a[r[i]]);
		//如果不写i<r[i]那么会被交换两次，即交换回来
	}
	for (int i=2;i<=len;i*=2){
		Complex wn(cos(flag*2*PI/i),sin(flag*2*PI/i));//求出w_i^1
		for (int j=0;j<len;j+=i){
			Complex w(1,0);
			for (int k=j;k<j+(i/2);k++){
				Complex u=a[k],v=w*a[k+(i/2)];//千万别忘记乘上w
				a[k]=u+v,a[k+(i/2)]=u-v;//套公式
				w=w*wn;//实时维护w_i^(k-j)
			}
		}
	}
	if (flag==-1){
		for (int i=0;i<len;i++)  a[i].x/=len;//如果是在进行IDFT,每个复数均除以长度
	}
}

signed main(){
	n=read(),m=read();
	for (int i=0;i<=n;i++)  a[i].x=read();
	for (int i=0;i<=m;i++)  b[i].x=read();
	
	while (len<=n+m)  len*=2,cnt++;//n必须是2的整数次方，否则单位根就没有上面的性质了
	for (int i=0;i<len;i++)  r[i]=(r[i>>1]>>1)+((i&1)<<(cnt-1));//数位dp
	
	FFT(a,len,1),FFT(b,len,1);//进行两次DFT
	for (int i=0;i<=len;i++)  a[i]=a[i]*b[i];//点值相乘
	FFT(a,len,-1);//IDFT
	
	for (int i=0;i<=n+m;i++)  cout<<(int)(a[i].x+0.5)<<' ';//保留整数输出
	return 0;
}//by ducati

Summary

可以发现，FFT在优秀的$O(n\log n)$的复杂度内完成了多项式乘法。同时，它也能应对许多实际问题，这些东西留到多项式学习笔记[三]来详细讲解。

但是不难发现，FFT中用到了三角函数与double运算，会存在不小的精度误差。同时，对于一些需要取模且在不取模的情况下系数很大的题目，FFT将非常难运行出正确的结果。所以我们需要一个精度更高，能应对取模的多项式乘法算法。

前置芝士: 原根

何为原根？

定义如下: 正整数$g$是正整数$n$的原根，当且仅当$1\le g\le n-1$ ，且$g$模 $n$的阶为$\phi (n)$。

更简单地说，原根$g$在$[1,n)$范围内，且使得$g^t$模$n$余$1$的自然数$t$只有$\phi (n)$。

原根有许多优美的性质，这里简要地列出几条：
(1)$g^1,g^2\cdots g^{\phi (n)}$两两 不同，所以$g$可以作为在模$n$意义下$\log$的底数；
(2)原根的大小在$n^{0.25}$级别。
(3)令$n$的最小原根为$p$，那么所有其他的原根均为$p^r(\gcd (r,\phi (n))=1)$在模$n$意义下的值。

所以，我们可以采用如下方法计算出一个自然数的所有原根：

1. 枚举$i(1\le i\le n-1)$；
   (1)对于所有满足$d|\phi (n)$的自然数$d$，通过快速幂求出$i^d$模$n$的值$val$；
   ①如果$val=1$且$d\ne \phi (n)$，那么这个数不是原根；
   ②如果$val\ne 1$且$d=\phi (n)$，那么这个数不是原根；
   (2)如果这个数是原根，那么令$p=i$并退出循环。

2. 枚举$r$使得$1\le r\le \phi (n)-1$且$\gcd (\phi (n),r)=1$，求出$p^r$模$n$的值从而得到所有原根。

一些经典的数的原根如下:
①$998244353$的原根为$3$；
②$10^9+7$的原根为$5$。

现在我们已经搞懂了原根已经求法并记得了一些重要的数的原根。但是这东西对多项式乘法有什么用呢？

NTT

回顾一下FFT我们是如何定义$w_n^k$的。现在，我们将$w_n^k$的定义魔改一番：$$w_n^k=(x^{\frac{p-1}{n}}{)}^k$$

这里的$p$表示多项式乘法的模数，$x$表示$p$的一个原根。

这么定义之后，不难发现单位根拥有的所有性质都仍然满足。下面给出证明：

Lemma 1: $n$个“单位根”两两不同
Certification 1: 根据原根的定义显然。

Lemma 2: $w_n^k=w_{2n}^{2k}$
Certification 2: $(x^{\frac{p-1}{2n}}{)}^{2k}=(x^{\frac{p-1}{n}}{)}^k$

Lemma 3: $(w_n^k{)}^n=1$
Certification 3: $(x^{\frac{p-1}{n}}{)}^{kn}=(x^{p-1}{)}^k=1$

Lemma 4.1: $w_n^{\frac{n}{2}}=-1$
Certification 4.1:
∵$((x^{\frac{p-1}{n}}{)}^{\frac{n}{2}}{)}^2=1$
又∵$w_n^{\frac{n}{2}}\ne w_n^n$
∴$w_n^{\frac{n}{2}}=-1$

Lemma 4.2: $w_n^k=-w_n^{k-\frac{n}{2}}$
Certification 4.2:

$(x^{\frac{p-1}{n}}{)}^{k-\frac{n}{2}}=\frac{(x^{\frac{p-1}{n}}{)}^k}{(x^{\frac{p-1}{n}}{)}^{\frac{n}{2}}}=\frac{w_n^k}{-1}=-w_n^k$

Lemma 5: $w_n^a\times w_n^b=w_n^{a+b}$
Certification 5:
$(x^{\frac{p-1}{n}}{)}^a\times (x^{\frac{p-1}{n}}{)}^b=(x^{\frac{p-1}{n}}{)}^{a+b}=w_n^{a+b}$

综上所述，我们只需要把FFT中所有的$w_n^k$的定义改一改，其他的就基本相同了。

Code: NTT

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int maxl=5000005,mod=998244353;
const int inv3=(mod+1)/3;//998244353的一个原根是3, 而3在模998244353意义下的值为inv3

int read(){
	int s=0,w=1;char ch=getchar();
	while (ch<'0'||ch>'9'){if (ch=='-')  w=-w;ch=getchar();}
	while (ch>='0'&&ch<='9'){s=(s<<1)+(s<<3)+(ch^'0');ch=getchar();}
	return s*w;
}

int n,m,len=1,cnt;
int r[maxl],a[maxl],b[maxl];

int quick_power(int x,int y){
	int res=1;
	for (;y;y=y>>1,x=(x*x)%mod){
		if (y&1)  res=(res*x)%mod;
	}
	return res;
}

void NTT(int *a,int len,int flag){
	for (int i=0;i<=len;i++){
		if (i<r[i])  swap(a[i],a[r[i]]);
	}
	for (int i=2;i<=len;i*=2){
		int wn=1;
		if (flag==1)  wn=quick_power(3,(mod-1)/i);
		else wn=quick_power(inv3,(mod-1)/i);
		
		for (int j=0;j<len;j+=i){
			int w=1;
			for (int k=j;k<j+(i/2);k++){
				int u=a[k],v=(w*a[k+(i/2)])%mod;
				a[k]=(u+v)%mod,a[k+(i/2)]=(u+mod-v)%mod;//注意减法的取模
				w=(w*wn)%mod;
			}
		}
	}
	if (flag==-1){
		int p=quick_power(len,mod-2);//除以一个数相当于乘它的逆元
		for (int i=0;i<=len;i++)  a[i]=(a[i]*p)%mod;
	}
}

signed main(){
	n=read(),m=read();
	for (int i=0;i<=n;i++)  a[i]=read();
	for (int i=0;i<=m;i++)  b[i]=read();
	
	while (len<=n+m)  len*=2,cnt++;
	for (int i=0;i<=len;i++)  r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1));
	
	NTT(a,len,1),NTT(b,len,1);
	for (int i=0;i<=len;i++)  a[i]=(a[i]*b[i])%mod;
	
	NTT(a,len,-1);
	for (int i=0;i<=n+m;i++)  printf("%lld ",a[i]);
	return 0;
}

Summary

NTT可以应对自如许多取模的多项式乘法题，并且没有精度误差。

但是，我们回顾一下在NTT中$w_n^k=$的定义$w_n^k=(x^{\frac{p-1}{n}}{)}^k$。

只有当$n|(p-1)$的时候NTT才能正确运行！

由于$n$总是$2$的整数次方，所以$p-1$必须可以写成一个$2^r\times d+R$的形式，其中$\lceil {\log }_{2}n\rceil +1\le r$。对于这一类的模数$p$我们称之为NTT模数，如：

①$167772161$
②$469762049$
③$998244353$
④$1004535809$

$\cdots \cdots$

很可惜，有的时候模数并不是NTT模数，这导致FFT与NTT都都较难运行出正确的答案。我们需要一个更牛逼的做法来解决这种毒瘤的问题。

它，有一个好听的名字，叫做任意模数NTT。

多项式学习笔记[二]（全网最详细！有图有代码有解释有例题有总结！)

仅做个人学习总结使用
参考的文章：
https://www.jianshu.com/p/1f6025995502
https://zhuanlan.zhihu.com/p/69471608

相信大家都不会陌生，经常遇见含有这些分式的积分类型
现在说说有哪些技巧可以简单应付

一个真分式，分子的次数 < 分母的次数
我们把一个真分式拆解为几个小分式，通常第一步会先把分母进行因式分解，然后按照那个因式分裂为小分式

对于小分式，分子的次数 总会 比分母的次数少1次方：deg(分子) = deg(分母) - 1
例如分母是二阶 $a{x}^2+bx+c$，则分子为$Ax+B$
若分母是一阶 $ax+b$，则分子为常数A

不过，对于高阶极点来说，小分式的个数 = 分母的因式个数
例如$(x+5{)}^3$，因式为$(x+5{)}^3$，$(x+5{)}^2$，$(x+5)$，共三个因式
$(x^2+4 )4，因式为(x2+4)4，(x2+4)3，(x2+4)2，(x^2+4)，共四个因式

常用的方法无非都是那几种：
添项减项法：这个方法对

$\frac{1}{(x+a)(x+b)}$

$1/[(x+a)(x+b)]$

型有效
待定系数法：即小分式通分后，把分子与原式的分子恒等，从而解出对应系数
留数法：即通过消去零因式来解出系数，分母要求为线性(ax+b)型因式，可以是高阶极点
这个方法其实跟z变换类似

添项减项法 和 待定系数法：

留数法：

留数法对于一次因式，一阶极点的因式时最好用的
例如：

而待定系数法，则需要对联立多元方程有很好的运算技巧
通常对于二次或以上的因式最好用
例如：

下面是练习，你们可以试试：

如果分子的次数 ≥ 分母的次数，这是假分式，设法自然会有些改变

这个可依旧运用待定系数法：

或者多项式除法：

图5留数法第二道题解答有点错误，应该对A那个式子求导。希望作者更正下。☺☺☺

https://zhuanlan.zhihu.com/p/69471608

分解步骤总览：

1. 判别真假分式.

2. 真分式分解出待定式.

3. 待定系数求解方法: 实根法(一次式), 复根法(二次式), 求导法(一次n重), 极限法(一、二次的二重)

4. 判别真假分式

形如 $\frac{P(x)}{Q(x)}$
的分式, 若分子指数等于或高于分母, 则要化为真分式.[1]

化简方法: 做多项式除法[2]

例如:

2. 真分式分解

[公式] 重一次因式

形如: $\frac{P(x)}{(Ax+b{)}^k}$

当 $k=1$时

其中$A_i$为待定系数.

当 $k>1$ 时

例如:

$A_i,B_i,C_i$ 为待定系数.

$k$ 二次因式

形如

$\frac{P(x)}{(a{x}^2+bx+c{)}^k}$

当$k=1$时,

当$k>1$时,

例如:

3. 待定系数求解[3]

   无特征——反解方程法

将各项通分合并, 将分子与原式的分子做系数比对, 写出关于待定系数的方程, 进行求解

多个不同的一次式, 且无重因式——实根代入法