FFT

发表于 2019-12-25 分类于总结

FFT有关的知识点总结。

前置知识

多项式求导

根据 $(x^n)’=nx^{n-1}$ 可推出。

多项式不定积分

根据 $\int x^ndx=\frac{x^{n+1}}{n+1} +C$ 可推出。

多项式加法

对应系数相加。

多项式减法

对应系数相减。

多项式乘法

DFT与IDFT的原理

给定两个多项式 $A(x)=\sum_{i=0}^na_ix^i,B(x)=\sum_{i=0}^nb_ix^i$，则两个多项式相乘后是：$(A * B)(x)=(a_0b_0)+(a_1b_0+a_0b_1)x+\dots+(a_nb_n)x^{2n}$

一个 $n$ 次多项式 $f(x)=\sum_{i=0}^na_ix^i$ 可以由 $n+1$ 个点 $(x_i,f(x_i))$ 唯一表示，那么FFT就是把两个多项式 $A(x),B(x)$（假设系数都为 $n$，不够最高位补0即可）的系数表示法转成由大于 $2n$ 个（因为相乘后多项式的次数是 $2n$的，要用大于 $2n$ 个点值才能表示）相同的点值（即$(x_i,A(x_i)),(x_i,B(x_i))$）表示，然后对应函数值相乘，得到点值：$(x_i,A(x_i) \times B(x_i))$，最后再把点值表示换成系数表示的过程。

如果随便取一些 $x_i$，那么每次操作是 $O(n^2)$ 的，毫无卵用。。。

假设我们要取 $n$ 个点值（不妨设 $n$ 为偶数）。

考虑 $x_k$ 为方程 $x^n=1$ 的复数解，即 $x_k=\omega _ n ^ k = e ^ {\frac{2k\pi i}{n}}=\cos(\frac{2k\pi}{n})+i\sin(\frac{2k\pi}{n})$

取 $1$ 的 $n$ 次单位根的作用在于，这个 $\omega _ n ^ k$ 有如下性质：

1，$\omega _ n^k=\omega _ n^{k \bmod n}$

2，若 $n,k$ 为偶数，则 $\omega _ n^k=\omega _ {n/2}^{k/2}$

证明都十分显然。

我们先回到原问题上面，现在要求的是 $f(\omega _ n^0),f(\omega _ n^1)…f(\omega _ n^{n-1})$，考虑分治，把多项式中次数为偶数的系数放在 $f_0$ 中，奇数的放在 $f_1$ 中，即：
$$f_0(x)=\sum_{i=0}^{2i\leq n}a_{2i}x^i\f_1(x)=\sum_{i=0}^{2i+1\leq n}a_{2i+1}x^i$$
那么有：$f(\omega _ n^k)=f_0((\omega _ n^k)^2)+(\omega _ n^k)f_1((\omega _ n^k)^2)$

又因为 $(\omega _ n^k)^2=\omega _ n^{2k}=\omega _ {n/2}^k$

所以：$$f(\omega _ n^k)=f_0(\omega _ {n/2}^k)+(\omega _ n^k)f_1((\omega _ {n/2}^k))\ (0\leq k < n)$$

当 $0\leq k < \frac{n}{2} $ 时，有：
$$f(\omega _ n^k)=f_0(\omega _ {n/2}^k)+(\omega _ n^k)f_1((\omega _ {n/2}^k))\f(\omega _ n^{k+\frac{n}{2}})=f_0(\omega _ {n/2}^k)-(\omega _ n^k)f_1((\omega _ {n/2}^k))$$

这样就变成了两个子问题，递归求解即可。时间复杂度 $O(n\log n)$

现在我们已经将多项式用系数表示转成点值表示(dft)了，那么这个操作的逆变换(idft)怎么做呢？

考虑dft的过程，可以变成矩阵相乘的形式：

$\begin{bmatrix}
f(\omega _ n^0)\
f(\omega 1)\
f(\omega 2) \
\vdots\
f(\omega _ n^{n-1})
\end{bmatrix}

\begin{bmatrix}
1 & 1 & 1 & \dots & 1\
1 & \omega _ n^1 & \omega _ n^2 & \dots &\omega _ n^{n-1} \
1 & \omega _ n^2 & \omega _ n^4 & \dots & \omega _ n^{2(n-1)}\
\vdots & \vdots & \vdots & \ddots & \vdots\
1 & \omega _ n^{n-1} & \omega _ n^{2(n-1)} & \dots & \omega _ n^{(n-1)(n-1)}
\end{bmatrix}
\begin{bmatrix}
a_0\
a_1\
a_2\
\vdots\
a_{n-1}
\end{bmatrix}\$

将上述表达式从左到右的矩阵记为 $Y,V,A$，则有 $Y=VA$。

因为是逆过程，所以要找到 $V$ 的逆矩阵使得 $V^{-1}Y=A$.

可以发现~~我也不知道是怎么发现的~~ ，下面这个矩阵就是 $V$ 的逆矩阵：

$V^{-1}
=\frac{1}{n}
\begin{bmatrix}
1 & 1 & 1 & \dots & 1\
1 & \omega _ n^{-1} & \omega _ n^{-2} & \dots &\omega _ n^{-(n-1)} \
1 & \omega _ n^{-2} & \omega _ n^{-4} & \dots & \omega _ n^{-2(n-1)}\
\vdots & \vdots & \vdots & \dots & \vdots\
1 & \omega _ n^{-(n-1)} & \omega _ n^{-2(n-1)} & \dots & \omega _ n^{-(n-1)(n-1)}
\end{bmatrix}
$

证明的话…可以考虑暴力矩阵乘法：$C_{ij}=\sum_{k}V_{ik}V^{-1} _ {kj}$，把两个矩阵暴力乘起来后可以发现当 $i=j$ 时，$C_{ii}$ 的值为 $n$ ，否则为 $0$.

那么idft的过程就与dft的类似了。

所以只需要将序列分别dft一下，然后对应相乘，最后idft回去就可以啦~

到现在为止我们能做一个递归版的 fft 啦！

但时实测这样的速度非常慢，需要写非递归的。

考虑 fft 的过程，经过一堆分治后，系数 $a_i$最终会到达一个位置上，而这些位置跟 $i$ 有关，如果 $i$ 是偶数就会被分到 $f_0(x)$ 中（即左边），否则会被分到右边，因此 $i$ 最后的位置是将 $i$ 在二进制表示下将位翻转后的位置，预处理即可。

到现在为止我们就会做一个非递归版的 fft 啦！

NTT

那如果要对大质数取模呢？

我们直接跑FFT，最后取个模不就可以了？

但是FFT会有比较大的精度问题，当数字比较大的时候很容易挂掉。

考虑这个质数的一个原根 $g$，它的性质与 $n$ 次单位根相似，即 $g^{P-1} \equiv 1(\bmod p)$，但因为我们必须要保证每次做 dft 时次数 $n$ 必须是偶数，所以一开始的 $n$ 必须是2的次幂。所以要满足：$2^k | (P-1)$。

如当 $P=998244353$ 时， $k$ 最大是 $23$，即 $n$ 最大是 $8388608‬$。这就是快速数论变换（NTT）了。此时需满足 $P$ 为质数，且满足 $2^k|(P-1)$ 中最大的 $k$ 要满足 $n\leq 2^k$。

能跑NTT的一些模数有：$998244353,1004535809,469762049$ 等，这三个数的原根都有 $3$。

到现在为止我们就会做一个非递归版的对模数有限制的 ntt 啦！

一个小模板：

#define G 3//原根
#define mod 998244353ll//模数
#define ll long long
int R[N],L;
int pre_ntt(int n)
{
	int m;
	for(m=1,L=0;m<=n;m<<=1) L++;
	for(int i=1;i<m;i++) R[i] = (R[i>>1]>>1)|((i&1)<<(L-1));
	return m;
}
void ntt(LL *a,int n,int t)//t=1时为dft,-1时为idft
{
	for(int i=1;i<n;i++) if(i>R[i]) swap(a[i],a[R[i]]);
	for(int i=1;(i<<1)<=n;i<<=1)
	{
		ll wn=Pow(G,(mod-1)/(i<<1));
		if(t==-1) wn=Pow(wn,mod-2);
		for(int j=0;j<n;j+=(i<<1))
		{
			ll w=1;
			for(int k=0;k<i;k++,w=w*wn%mod)
			{
				ll x=a[j+k],y=a[i+j+k]*w%mod;
				a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
			}
		}
	}
	if(t==1) return;
	ll invn=Pow(n,mod-2);
	for(int i=0;i<=n;i++) a[i]=a[i]*invn%mod;
}

任意模数FFT

那如果模数不满足上面的限制呢？

这时有两种方法：

三模数FFT

先拿几个能做NTT的模数跑一次 NTT，然后用中国剩余定理合并。

显然选三个比较大的模数就可以了。

拆系数FFT

设 $m=\left \lfloor \sqrt{P} \right \rfloor$，将FFT中每个数表示成 $bm+a$ 的形式，且 $a<m$。那么就可以将两个数的乘法拆成四部分：$(bm+a)(dm+c)=bdm^2+(bc+ad)m+ac$。然后分别做普通的 FFT ，这样就避免了精度问题。此时一共需要 $4$ 次DFT，$3$ 次IDFT。常数大到爆炸。

考虑优化DFT，将两次dft合并到一次。

假设现在要求 $dft(a)$ 与 $dft(b)$。

设 $f_k=a_k+ib_k,g_k=a_k-ib_k$。

将有 $dft(f)k=\text{conj}(dft(g){n-k})$。其中 $\text{conj}$ 表示共轭。

证明：$dft(f)k=\sum{j=0}^{n-1} w_n^{kj}f_j=\sum_{j=0}^{n-1} \text{conj}(w_n^{-kj}g_j)=\sum_{j=0}^{n-1} w_n^{kj}\text{conj}(g_j)=\text{conj}(dft(g)_{n-k})$。

于是 $dft(a)_k=\frac{dft(f)_k+dft(g)_k}{2}$ ，$dft(b)_k=-i\times \frac{dft(f)_k-dft(g)_k}{2}$。

这样做的前提条件显然是 $a,b$ 中的数均为实数。

考虑优化IDFT，设 $h=dft(a)+i\times dft(b)$。

那么 $idft(h)$ 的实部即为 $a$，虚部即为 $b$。

于是就优化成了分别 $2$ 次的DFT和IDFT了。

模板：

const ll mo=1e9+7;

db pi=acos(-1.);
struct P{
	db x,y;
	P(db _x=0,db _y=0) {x=_x,y=_y;}
	friend inline P operator +(const P&A,const P&B)
	{
		return (P){A.x+B.x,A.y+B.y};
	}
	friend inline P operator -(const P&A,const P&B)
	{
		return (P){A.x-B.x,A.y-B.y};
	}
	friend inline P operator *(const P&A,const P&B)
	{
		return (P){A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x};
	}
	friend inline P operator /(const P&A,const db &x)
	{
		return (P){A.x/x,A.y/x};
	}
};
P conj(P A) {return (P){A.x,-A.y};}

const int M=1<<20;
P W[M]; int R[M];
inline void PolyInit()
{
	for(int i=1;i<M;i<<=1)
		fo(j,0,i-1)
			W[i+j]=(P){cos(pi*j/i),sin(pi*j/i)};
}
typedef vector<P> Poly;
inline void ntt(P *a,int n,int t)
{
	fo(i,0,n-1)
	{
		R[i]=(R[i>>1]>>1)|((i&1)*(n>>1));
		if(i<R[i]) swap(a[i],a[R[i]]);
	}
	P w;
	for(int i=1;i<n;i<<=1)
		for(int j=0;j<n;j+=(i<<1))
			for(int k=0;k<i;k++)
				w=W[i+k]*a[i+j+k],
				a[i+j+k]=a[j+k]-w,
				a[j+k]=a[j+k]+w;
	if(t==1) return;
	reverse(a+1,a+n);
	fo(i,0,n-1) a[i]=a[i]/n;
}
inline void ntt(Poly &A,int n,int t){ntt(&A[0],n,t);}
inline Poly operator +(Poly A,Poly B)
{
	A.resize(max(A.size(),B.size()));
	fo(i,0,B.size()-1) A[i]=A[i]+B[i];
	return A;
}
#define Real(A) ((ll)floor(A.x+0.5))
#define Imag(A) ((ll)floor(A.y+0.5))
inline Poly operator *(Poly A,Poly B)
{
	static Poly C,D;
	int n=A.size(),m=B.size(),k=n+m-1,len=1;
	for(;len<k;len<<=1);
	C.resize(len); D.resize(len); A.resize(len); B.resize(len);
	ff(i,0,len) C[i]=(P){Real(A[i])&32767,Real(A[i])>>15},D[i]=(P){Real(B[i])&32767,Real(B[i])>>15};
	ntt(C,len,1); ntt(D,len,1);
	int j;
	ff(i,0,len)
	{
		P d4,d0,d1,d2,d3;
		j=(len-i)&(len-1);
		d4=conj(C[j]);
		d0=(d4+C[i])*P(0.5,0);
		d1=(d4-C[i])*P(0,0.5);
		d4=conj(D[j]);
		d2=(d4+D[i])*P(0.5,0);
		d3=(d4-D[i])*P(0,0.5);
		A[i]=d0*d2+d1*d3*P(0,1);
		B[i]=d0*d3+d1*d2;
	}
	ntt(A,len,-1); ntt(B,len,-1);
	ff(i,0,len) C[i]=(Real(A[i]) + (Imag(A[i]) % mo << 30) + (Real(B[i]) % mo << 15))%mo;
	C.resize(k);
	return C;
}

一些技能

多项式求逆

给一个$n$次多项式$A(x)$，求$B(x)$满足$A(x)* B(x)\equiv 1\left ( \bmod x^n \right )$

不妨设$n$为2的次幂：

当$n=1$时，求一个整数的逆就好了。

当$n\geq 2$时，

假设已经求出了$A(x)* C(x)\equiv 1\left ( \bmod x^{\frac{n}{2}} \right )$

$\because A(x)* B(x)\equiv 1\left ( \bmod x^n \right )$

$\therefore A(x)* B(x)\equiv 1\left ( \bmod x^{\frac{n}{2}} \right )$

$\therefore A(x)* (B(x)-C(x))\equiv 0\left ( \bmod x^{\frac{n}{2}} \right )$

又$\because A(x)\not\equiv 0\left ( \bmod x^n \right )$

$\therefore (B(x)-C(x))\equiv 0\left ( \bmod x^{\frac{n}{2}} \right )$

仔细想一想可以发现$\left (B(x)-C(x)\right )^2\equiv 0\left ( \bmod x^{n} \right )$

化简得$B^{2}(x)+C^{2}(x)-2B(x)C(x)\equiv 0\left ( \bmod x^{n} \right )$

同时乘上$A(x)$得：$B(x)+A(x)* C^{2}(x)-2C(x)\equiv 0\left ( \bmod x^{n} \right )$

移项得：$B(x)\equiv C(x)\left ( 2-A(x)* C(x)\right )\left ( \bmod x^{n} \right )$

多项式乘法+倍增即可。

时间复杂度$T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)$

多项式除法

给两个次数分别为$n$和$m$的多项式$F(x),G(x)$，求$Q(x),R(x)$满足$F(x)=G(x)* Q(x)+R(x)$,且$Q(x)$的次数为$n-m$

$\because F(x)\equiv G(x)* Q(x)+R(x)\left ( \bmod x^{n+1}\right )$

$\therefore F(\frac{1}{x})\equiv G(\frac{1}{x})* Q(\frac{1}{x})+R(\frac{1}{x})\left ( \bmod x^{n+1}\right )$

同时乘以$x^n$得：$ F^{rev}(x)\equiv G^{rev}(x)* Q^{rev}(x)+R^{rev}(x)\times x^{n-dR}\left ( \bmod x^{n+1}\right )$(其中$rev$表示将多项式系数翻转，$dR$表示$R$这个多项式的次数，其余同理)

$\because dR < m $

$\therefore n-dR>n-m$

因此在$\bmod x^{n-m+1}$意义下，$R^{rev}(x)\times x^{n-dR}\equiv 0$

$\therefore F^{rev}(x)\equiv G^{rev}(x)* Q^{rev}(x) \left ( \bmod x^{n-m+1}\right )$

多项式求逆+翻转即可。

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

多项式取模

求出除法以后，取模变得非常简单。

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

多项式开方

给一个$n$次多项式$A(x)$，求$B(x)$满足$A(x)\equiv B^{2}(x)\left ( \bmod x^n \right )$

还是考虑倍增：

不妨设$n$为2的次幂：

当$n=1$时，求$x^2\equiv a\left ( \bmod p\right )$的解~~就好了~~。比较多的题目会保证常数项为1.

当$n\geq 2$时，

假设已经求出了$A(x)\equiv C^2(x)\left ( \bmod x^{\frac{n}{2}} \right )$，

易得$\left ( C(x)-B(x)\right ) * \left ( C(x)+B(x)\right )\equiv 0\left ( \bmod x^{\frac{n}{2}} \right )$

因此会有两个解，设$C(x)-B(x)\equiv 0\left ( \bmod x^{\frac{n}{2}} \right )$，则有$ C(x)+B(x)\equiv 0\left ( \bmod x^{\frac{n}{2}} \right )$

所以$ \left ( C(x)+B(x)\right )^2\equiv 0\Rightarrow C^2(x)+B^2(x)+2B(x)* C(x)\equiv 0\Rightarrow B(x)\equiv-\frac{C^2(x)+A(x)}{2C(x)}\left ( \bmod x^{n} \right )$

多项式求逆+乘法+倍增即可。

时间复杂度$T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)$

多项式ln

为啥多项式可以求ln啊。。。无法理解。。。

先来看看泰勒展开：

$$e^x=\sum_{i=0}^{\infty}\frac{x^i}{i!}$$

那么多项式求ln实际上就是知道$G(x)=e^{F(x)}$，让你求$F(x)$的过程。

对等式两边同时求$ln$得$\ln(G(x))=F(x)$

求导得：$\frac{G’(x)}{G(x)}=F’(x)$

求积分得：$\int {\frac{G’(x)}{G(x)}}=F(x)$

没了。。。

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

多项式exp

能求ln那当然能求exp啦。

$F(x)=e^{A(x)}$

求ln得：$\ln F(x)=A(x)$

$\ln F(x)-A(x)=0$

考虑牛顿迭代：

$F_{i+1}(x)=F_i(x)-\frac{\ln F_{i}(x)-A(x)}{\left (\ln F_{i}(x)-A(x) \right )’}=F_{i}(x)-\frac{\ln F_{i}(x)-A(x)}{\frac{1}{F_{i}(x)}}=F_{i}(x)\left ( 1-\ln F_{i}(x)+A(x)\right )\left ( \bmod x^n \right )$

每次迭代 $n$ 会翻一倍，时间复杂度$T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)$

多项式求幂

$F(x)=A^{k}(x)$

两边同时取对数得：$\ln F(x)=k \ln A(x)$

同时取指数得：$F(x)=\exp\left ( k\ln A(x)\right )$

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

多项式三角函数

不是很清楚这种东西有什么用…

求 $\sin(A),\cos(A)$ ：

直接硬上欧拉公式：

$$e^{ix}=\cos(x)+i\sin(x)\ e^{-ix}=\cos(x)-i\sin(x)$$

解得：

$$\cos(x)=\frac{e^{ix}+e^{-ix}}{2}\\sin(x)=\frac{e^{ix}-e^{-ix}}{2}$$

但是在模意义这个 $i$ 该怎么办呢？

考虑 $i^2=-1$，所以在模 $mod$ 的意义下这个 $i$ 就是 $mod-1$ 的二次剩余。

多点求值

“详见-多点求值”

多点插值

常系数齐次线性递推

“详见”

模板

#include <bits/stdc++.h>
using namespace std;
#define FO(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
#define fo(i,j,k) for(int i=(j),end_i=(k);i<=end_i;i++)
#define ff(i,j,k) for(int i=(j),end_i=(k);i< end_i;i++)
#define fd(i,j,k) for(int i=(j),end_i=(k);i>=end_i;i--)
#define DEBUG(x) cerr<<#x<<"="<<x<<endl
#define all(x) (x).begin(),(x).end()
#define cle(x) memset(x,0,sizeof(x))
#define lowbit(x) ((x)&-(x))
#define ll long long
#define ull unsigned ll
#define db double
#define lb long db
#define pb push_back
#define mp make_pair
#define fi first
#define se second
inline int read()
{
	int x=0; char ch=getchar(); bool f=0;
	for(;ch<'0'||ch>'9';ch=getchar()) if(ch=='-') f=1;
	for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<3)+(x<<1)+(ch^48);
	return f?-x:x;
}
#define CASET fo(___,1,read())

const ll mod=998244353;
inline ll Add(ll x,ll y){x+=y; return (x<mod)?x:x-mod;}
inline ll Dec(ll x,ll y){x-=y; return (x<0)?x+mod:x;}
inline ll Mul(ll x,ll y){return x*y%mod;}
inline ll Pow(ll x,ll y)
{
	ll ans=1;for(;y;y>>=1,x=x*x%mod)if(y&1) ans=ans*x%mod;
	return ans;
}

const int M=1<<19;
ll W[M];
inline void PolyInit()
{
	ll w;
	for(int i=1;i<M;i<<=1)
	{
		W[i]=1; w=Pow(3,(mod-1)/2/i);
		fo(j,1,i-1) W[i+j]=W[i+j-1]*w%mod;
	}
}
typedef vector<ll> Poly;
inline void print(Poly A)
{
	ff(i,0,A.size()) printf("%d ",A[i]);
	printf("\n");
}
inline void ntt(ll *a,int n,int t)
{
	static int R[M];
	fo(i,0,n-1)
	{
		R[i]=(R[i>>1]>>1)|((i&1)*(n>>1));
		if(i<R[i]) swap(a[i],a[R[i]]);
	}
	ll w;
	for(int i=1;i<n;i<<=1)
		for(int j=0;j<n;j+=(i<<1))
			for(int k=0;k<i;k++)
				w=W[i+k]*a[i+j+k]%mod,
				a[i+j+k]=Dec(a[j+k],w),
				a[j+k]=Add(a[j+k],w);
	if(t^1)
	{
		reverse(a+1,a+n);
		w=Pow(n,mod-2);
		fo(i,0,n-1) a[i]=w*a[i]%mod;
	}
}
inline void ntt(Poly &A,int n,int t){ntt(&A[0],n,t);}
inline Poly operator +(Poly A,Poly B)
{
	A.resize(max(A.size(),B.size()));
	fo(i,0,B.size()-1) A[i]=Add(A[i],B[i]);
	return A;
}
inline Poly operator -(Poly A,Poly B)
{
	A.resize(max(A.size(),B.size()));
	fo(i,0,B.size()-1) A[i]=Dec(A[i],B[i]);
	return A;
}
inline Poly df(Poly A)
{
	fo(i,1,A.size()-1) A[i-1]=A[i]*i%mod;
	A.resize(A.size()-1);
	return A;
}
inline Poly jf(Poly A)
{
	A.pb(0);
	fd(i,A.size()-1,1) A[i]=A[i-1]*Pow(i,mod-2)%mod;
	A[0]=0; return A;
}
inline Poly operator *(Poly A,ll k)
{
	fo(i,0,A.size()-1) A[i]=A[i]*k%mod;
	return A;
}
inline Poly operator *(Poly A,Poly B)
{
	int n=A.size(),m=B.size(),k=n+m-1,len=1;
	for(;len<k;len<<=1);
	A.resize(len); ntt(A,len,1);
	B.resize(len); ntt(B,len,1);
	fo(i,0,len-1) A[i]=A[i]*B[i]%mod;
	ntt(A,len,-1);
	A.resize(k);
	return A;
}
inline Poly operator ~(Poly f)
{
	int n=f.size();
	Poly g,h;
	g.pb(Pow(f[0],mod-2));
	int m=2;
	for(;m<n;m<<=1)
	{
		h.resize(m<<1); g.resize(m<<1);
		fo(i,0,m-1) h[i]=f[i];
		ntt(h,m<<1,1); ntt(g,m<<1,1);
		fo(i,0,(m<<1)-1) g[i]=Mul(2+mod-Mul(g[i],h[i]),g[i]);
		ntt(g,m<<1,-1); g.resize(m);
	}
	g.resize(m<<1); f.resize(m<<1);
	ntt(f,m<<1,1); ntt(g,m<<1,1);
	fo(i,0,(m<<1)-1) g[i]=Mul(2+mod-Mul(g[i],f[i]),g[i]);
	ntt(g,m<<1,-1); g.resize(n);
	return g;
}
inline Poly Ln(Poly A)
{
	int n=A.size();
	A=jf((~A)*df(A));
	A.resize(n); return A;
}
inline Poly Exp(const Poly &A)
{
	int n=1; for(;n<A.size();n<<=1);
	Poly B,C,D; B.clear(); B.pb(1);
	for(int m=2;m<=n;m<<=1)
	{
		C=B; C.resize(m); D=A; D.resize(m);
		C=D-Ln(C); C[0]=Add(C[0],1);
		B=B*C; B.resize(m);
	}
	B.resize(A.size()); return B;
}
inline Poly operator ^(Poly A,ll k)
{
	if(!A.size()) return A;
	ll tmp=A[0],w=Pow(tmp,k);
	tmp=Pow(tmp,mod-2);
	fo(i,0,A.size()-1) A[i]=A[i]*tmp%mod;
	A=Exp(Ln(A)*k);
	fo(i,0,A.size()-1) A[i]=A[i]*w%mod;
	return A;
}
inline Poly Cos(Poly A)
{
	static ll w4=Pow(3,(mod-1)/4);
	return (Exp(A*w4)+Exp(A*(mod-w4)))*((mod+1)/2);
}
inline Poly Sin(Poly A)
{
	static ll w4=Pow(3,(mod-1)/4);
	return (Exp(A*w4)-Exp(A*(mod-w4)))*(Pow(w4,mod-2)*((mod+1)/2)%mod);
}
inline Poly Sqrt(Poly A)
{
	Poly C,D,B(1,1);
	C.clear(); D.clear();
	int n=A.size();
	for(int l=4;(l>>2)<n;l<<=1)
	{
		C=A; C.resize(l>>1);
		D=B; D.resize(l>>1); D=(~D);
		C.resize(l); D.resize(l);
		ntt(C,l,1); ntt(D,l,1);
		ff(i,0,l) C[i]=C[i]*D[i]%mod;
		ntt(C,l,-1);
		B.resize(l>>1);
		ff(i,0,(l>>1)) B[i]=Add(C[i],B[i])*((mod+1)>>1)%mod;
	}
	B.resize(n);
	return B;
}
inline Poly operator/(Poly A,Poly B)
{
	int len=1,deg=A.size()-B.size()+1;
	reverse(all(A)); reverse(all(B));
	for(;len<=deg;len<<=1);
	B.resize(len); B=~B; B.resize(deg);
	A=A*B; A.resize(deg);
	reverse(all(A));
	return A;
}
inline Poly operator%(const Poly &A,const Poly &B)
{
	if(A.size()<B.size()) return A;
	Poly C=A-(A/B)*B;
	C.resize(B.size()-1);
	return C;
}
inline Poly Pow(Poly A,ll n,Poly M)
{
	Poly B=A; n--;
	for(;n;n>>=1,A=(A*A)%M) if(n&1ll) B=(B*A)%M;
	return B;
}
const int N=64005*4;
#define lc (u<<1)
#define rc ((u<<1)|1)
Poly P[N<<2];
ll ans[N],a[N];
inline Poly MulT(Poly A,Poly B)
{
	int n=A.size(),m=B.size();
	reverse(all(B));
	int len=1;
	for(;len<n;len<<=1);
	A.resize(len); B.resize(len);
	ntt(A,len,1); ntt(B,len,1);
	ff(i,0,len) A[i]=A[i]*B[i]%mod;
	ntt(A,len,-1);
	B.clear();
	len--;
	fo(i,m-1,n+m-2) B.pb(A[i&len]);
	return B;
}
void solve(int u,int l,int r)
{
	if(l==r)
	{
		P[u].pb(1); P[u].pb(Dec(0,a[l]));
		return;
	}
	int mid=(l+r)>>1;
	solve(lc,l,mid); solve(rc,mid+1,r);
	P[u]=P[lc]*P[rc];
}
void solve(int u,int l,int r,Poly A)
{
	A.resize(r-l+1);
	if(l==r) return (void)(ans[l]=A[0]);
	int mid=(l+r)>>1;
	solve(lc,l,mid,MulT(A,P[rc])); solve(rc,mid+1,r,MulT(A,P[lc]));
}
int n,m,k;
Poly F,G;
int main()
{
	PolyInit();
	n=read(); m=read(); k=max(n,m);
	fo(i,0,n) F.pb(read());
	F.resize(n+k+1);
	fo(i,1,m) a[i]=read();
	solve(1,1,k);
	F=MulT(F,(~P[1]));
	solve(1,1,k,F);
	fo(i,1,m) printf("%lld\n",ans[i]);
	return 0;
}

前置知识

多项式求导

多项式不定积分

多项式加法

多项式减法

多项式乘法

DFT与IDFT的原理

$\begin{bmatrix}f(\omega _ n^0)\f(\omega 1)\f(\omega 2) \\vdots\f(\omega _ n^{n-1})\end{bmatrix}

NTT

任意模数FFT

三模数FFT

拆系数FFT

一些技能

多项式求逆

多项式除法

多项式取模

多项式开方

多项式ln

多项式exp

多项式求幂

多项式三角函数

多点求值

多点插值

常系数齐次线性递推

模板

$\begin{bmatrix}
f(\omega _ n^0)\
f(\omega 1)\
f(\omega 2) \
\vdots\
f(\omega _ n^{n-1})
\end{bmatrix}