涉及了拉插，FFT，NTT，生成函数入门等多项式儿童套餐。

〇. 基本概念：

  • 多项式： 形如 f(x)=a0+a1⋅x+a2⋅x2+...+an⋅xnf(x)=a_0+a_1\cdot x+a_2\cdot x^2+...+a_n\cdot x^n 的 n 次函数。

  • 度（Degree）：对于一个多项式 f(x)f(x)，称其最高次项的次数为该多项式的度，记作 deg⁡x\deg x.

多项式四则运算： 加减乘很好记，对于除法运算我们同样引入两个概念：

  • 多项式的逆元（Inverse Element）： 对于多项式 f(x)f(x)，若存在 g(x)g(x) 满足：f(x)⋅g(x)≡1 (mod xn)f(x)\cdot g(x)\equiv 1~(mod~x^n)，
    则称 g(x)g(x) 为 f(x)f(x) 在模 xnx^n 意义下的逆元 ，记作 f(x)−1f(x)^{-1}。

  • 余数和商：对于 f(x),g(x)f(x),g(x) 存在 Q(x),R(x) (deg⁡R<deg⁡g)Q(x),R(x)~(\deg R<\deg g) 使得：

f(x)=Q(x)⋅g(x)+R(x)f(x)=Q(x)\cdot g(x)+R(x)

亦可表示为：

f(x)≡R(x) (mod g(x))f(x)\equiv R(x)~(mod~g(x))

  • 幂级数：指在级数的每一项均为与级数项序号 n 相对应的以常数倍的（x-a）的n次方。

说人话就是有无穷项的多项式。而我们主要探讨的是 形式多项式，即关注其 系数序列，主元 x 只是走个形式。

  • 生成函数：可简单理解为 ∑iai∗fi(x)\sum_i a_i*f_i(x) 的幂级数

其中fi(x)=xif_i(x)=x^i 时被称作 普通生成函数 ，fi(x)=xix!f_i(x)=\frac{x^i}{x!} 被称作指数生成函数，具体知识会在第四章介绍。

  • x 的 n 次下降幂：xn‾=∏0n−1(x−i)=x!(x−n)!x^{\underline n}=\prod_0^{n-1}(x-i) = \frac{x!}{(x-n)!}

壹. 拉格朗日插值：

1. 插值：

给定一系列点（共 n+1 个），确定一个 n 次多项式使得这些点在该函数图像上。
In other words, 就是将 点值表示 转换为 系数表示 的过程。

2. 拉插：

实质是对于每个点 (xi,yi)(x_i,y_i) 都凑一个式子，在联立求出所有的系数。
那怎么凑这个式子呢？由多项式 余数和商 的定义，我们不难给出一个引理：

f(x)≡f(xi)≡yi (mod (x−xi))f(x)\equiv f(x_i)\equiv y_i~(mod~(x-x_i))

上式可以作差证明。对每个点都列出该式，再套上 中国剩余定理 的公式，我们便可以得到拉插的一般形式：

f(x)=∑i=1n( yi⋅∏i≠jx−xjxi−xj )(mod M)f(x)=\sum_{i=1}^n (~y_i \cdot\prod_{i\neq j} \frac{x-x_j}{x_i-x_j}~) (mod~ M)

其中 MM 为 ∏(x−xi)\prod (x-x_i) ，注意到其中的变量可能有负数，取模注意一下。时间复杂度 O(n2)O(n^2).

3. 拉插特例：

当 xix_i 连续时，我们可以通过预处理阶乘，将时间复杂度降到 O(n)O(n)，尤其是可自选点值的题，往往通过选择 1→k+11\to k+1 来快速求解。

如例题： 求∑i=1nik\sum_{i=1}^n i^k ，已知 n≤1015,k≤106n\leq 10^{15},k\leq 10^6.

经推导可知 f(n)=∑i=1nikf(n)=\sum_{i=1}^n i^k 是一个 k+1 次多项式，我们取 1→k+21\to k+2作为点值带入，则原式变形为:

f(x)=∑i=1k+2( f(i)⋅∏i≠jn−ji−j )f(x)=\sum_{i=1}^{k+2} (~f(i) \cdot\prod_{i\neq j} \frac{n-j}{i-j}~)

对于其中各项：

  1. f(i)f(i) 可以边求值边计算得出；

  2. 后面一坨的分子可以预处理 前缀和 与 后缀和 ，把 i=ji=j 那一点挖掉求得；

  3. 而对于分母，我们可以注意到它对于任意一个 i 都是连续的：从负值 k+2-i 开始，到正值 i-1 结束，我们可以通过预处理阶乘将其拆解为： fac[k+2−i]⋅fac[i−1]fac[k+2-i]\cdot fac[i-1]，注意一下 k+2-i 若为奇数，就要多一个负号。

写的时候好像爆了 ll，套了一个龟速乘。

貳.快速傅里叶变换：

1. 概念：

给定两个不高于 n 次的多项式 f(x),g(x)f(x),g(x)，可以在 O(n⋅logn)O(n\cdot logn) 内求出 f(x)⋅g(x)f(x)\cdot g(x) 各个系数的算法。

其实质是利用 复数单位根 的种种性质，将两函数的系数表示转化为点值表示，在通过 逆变换 求出新函数的点值表示.

2. DFT 的实现：

对于两函数的系数表示，利用复数单位根性质将 ωn1,ωn2...ωnn\omega_n^1,\omega_n^2...\omega_n^n 带入，进行奇偶项分治，并最终求出整个函数对应点值的过程。

设 F(x)F(x) 为当前分治总函数的点值序列，O(x)O(x) 为奇数项分治后的序列，E(x)E(x) 为偶数项分治后的序列，根据性质：

F(ωnk)=O(ωn/2k)+ωnk×E(ωn/2k)F(ω^k_n)=O(ω^k_{n/2})+ω^k_n\times E(ω^k_{n/2})

F(ωnk+n/2)=O(ωn/2k)−ωnk×E(ωn/2k)F(ω^{k+n/2}_n)=O(ω^{k}_{n/2})-ω^k_n\times E(ω^{k}_{n/2})

不断向上合并答案，最终即可得到 f(x),g(x)f(x),g(x) 的点值序列，对应相乘即可。

3. IDFT 的实现：

即 DFT 的逆运算，将求出的 点值表示 转换为 系数表示。

可以通过 单位根求逆 或者 构造法 证明：把 DFT 中的 ωniω^i_n 换成 ωn−1ω_n^{-1}，做完之后除以 n 即可实现 IDFT。

#include <bits/stdc++.h>
#define FO(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
#define cle(x) memset(x,0,sizeof(x))
#define inf 2147483647
#define ll long long
using namespace std;

const int maxn=2100000;
const double pi=acos(-1);
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;
}

struct cp{
	cp (double X=0,double Y=0) {x=X; y=Y;}
	double x,y;
	cp operator + (cp const &b) const 
	{return cp(x+b.x,y+b.y);}
	cp operator - (cp const &b) const 
	{return cp(x-b.x,y-b.y);}
	cp operator * (cp const &b) const 
	{return cp(x*b.x-y*b.y,y*b.x+x*b.y);}
}f[maxn],g[maxn],y[maxn];
int n,m; 

void fft(cp *f,int l,int c)
{
	if (l==1) return;
	cp *L=f,*R=f+l/2;   
    //指针可直接 引用 原数组 ，从而快速对原函数进行奇偶分项 
	for (int i=0;i<l;i++) y[i]=f[i];
	for (int i=0;i<l/2;i++) L[i]=y[i<<1],R[i]=y[i<<1|1];   

	fft(L,l/2,c); fft(R,l/2,c);
	cp w0( cos(2*pi/l),c*sin(2*pi/l) ), b(1,0);
	for (int k=0;k<l/2;k++) {
		cp t=b*R[k]; b=b*w0;
		y[k]=L[k]+t; y[k+l/2]=L[k]-t;  
        //求出对应的点值表示 
        //用 Y，不能用 F 存！否则会覆盖L,R数组导致 F 的改变！
	}
	for (int i=0;i<l;i++) f[i]=y[i]; 
    //重新覆盖，原函数转化为点值，且不影响后续操作
}

int main()
{
	n=read(),m=read();
	for (int i=0;i<=n;i++) f[i].x=read();
	for (int i=0;i<=m;i++) g[i].x=read();
	for (m+=n,n=1;n<=m;n<<=1); 
	
	fft(f,n,1); 
	fft(g,n,1);  // 操作对象，长度，操作参数
	for (int i=0;i<=n;i++) f[i]=f[i]*g[i];
	fft(f,n,-1);
	for (int i=0;i<=m;i++) printf("%d ",(int)(f[i].x/n+0.49));
}

4. 位逆序置换优化：

我们注意到原版 FFT 进行了大量数组复制操作，加上是递归版本自带的小常数，很容易被卡。

通过分析递归底层的 fif_i 数组，我们惊讶地发现每个位置编号对应的正是原位置编号的 二进制翻转。

利用数位 DP 的思想，我们设 RxR_x 为 二进制数 x 进行二进制翻转后的大小，且长度是固定的 N=2nN=2^n （就像我们预先设计的那样）。那么得到如下转移方案：

  • 由于从小到大求解，我们已经得知了 Rx>>1R_{x>>1} 的大小，x>>1 此时代表着 二进制 x 除去个位后的值，而移位的同时导致最高位多了一个前导 0 。为了消除该影响，我们只需要将其再次移位即可。

  • 对于个位，若其为 1 我们要将其翻转至最高位，直接或上 2n−1=N/22^{n-1} = N/2 即可。

综上，转移方程为：

R(x)=R(x/2)>>1 ∣ (x mod 2)∗2n−1R(x)=R(x/2)>>1 ~|~ (x~\text{mod}~2)*2^{n-1}

同时我们也得到了非迭代版的FFT：

#include <cstring>
#include <cmath>
#include <cstdio>
#include <algorithm>
#define inf 2147483647
#define maxn 2100000
using namespace std;
const double pi=acos(-1);
struct cp{
	cp (double xx=0,double yy=0){x=xx,y=yy;}
	double x,y;
	cp operator + (cp const &B) const
	{return cp(x+B.x,y+B.y);}
	cp operator - (cp const &B) const
	{return cp(x-B.x,y-B.y);}
	cp operator * (cp const &B) const
	{return cp(x*B.x-y*B.y,x*B.y+y*B.x);}
}f[maxn],g[maxn],y[maxn];
int n,m,ro[maxn];

void fft(cp *f,int c) 
{
	for (int i=0;i<n;i++) if (i<ro[i]) swap(f[i],f[ro[i]]);  
	//直接得到递归至 1 的数组序列，非迭代向上转移 
	for (int h=2;h<=n;h<<=1)
	{
		cp w0( cos(2*pi/h), c*sin(2*pi/h) );
		for (int j=0;j<n;j+=h){
			cp buf(1,0);
			for (int i=j;i<(j+(h>>1));i++) {
				cp u=f[i]; cp t=buf*f[i+(h>>1)];
				f[i]=u+t; f[i+(h>>1)]=u-t;
				buf=buf*w0;
			}
		}
	}
} 

int main()
{
	scanf("%d%d",&n,&m);
	for (int i=0;i<=n;i++) scanf("%lf",&f[i].x);
	for (int i=0;i<=m;i++) scanf("%lf",&g[i].x);
	for (m+=n,n=1;n<=m;n<<=1);
	for(int i=0;i<n;i++) ro[i]=(ro[i>>1]>>1)|((i&1)?n>>1:0);
    //这里一定要带括号！！ 优先级问题已经挂了亿次了！
	fft(f,1); 
	for (int i=0;i<=m;i++) printf("%d ",(int)(f[i].x/n+0.49)); fft(g,1);
	for (int i=0;i<=n;i++) f[i]=f[i]*g[i];
	fft(f,-1);
	for (int i=0;i<=m;i++) printf("%d ",(int)(f[i].x/n+0.49));
}

叁. 快速数论变换：

一. 原根：

  1. 阶： 设 a,m∈N+a,m\in\mathbb{N^+}，且 a⊥ma\perp m，使 ax≡1(modm)a^x\equiv 1\pmod m 成立的 最小 正整数 xx，称为 aa 模 mm 的阶，记为 x=ordm ax=\text{ord}_m~a。

这篇博客给出了许多关于阶的重要性质。

  1. 原根的定义：设 g,m∈N+g,m\in N^+，且 g⊥mg\bot m，若 ordmg=φ(m)\text{ord}_mg=\varphi(m)，则称 gg 是模 mm 的原根。

不加证明地给出：膜 m 存在原根 的充要条件是m=2,4,pk,2pkm=2,4,p^k,2p^k.

很显然的是，如果 m 存在一个原根 g，那么对于φ(m)\varphi(m) 的任意质因数 pp 都有 gφ(m)p≠1(modm)g^{\frac{\varphi(m)}{p}}\neq 1\pmod m.

  1. 假设我们能快速地(?) 求出膜 m 下最小的原根 g，尝试着以此求出所有膜 m 的原根：

已知： gφ(m)≡1(modm),gφ(m)p≠1(modm)g^{\varphi(m)}\equiv 1\pmod m ,g^{\frac{\varphi(m)}{p}}\neq 1\pmod m.

那么对于 ∀s∈N+\forall s\in N^+，我们都可以来讨论 gsg^s 是原根的可能性：

  • 若 s>φ(m)s>\varphi(m)，由欧拉定理我们可将其转换为 gsmod  φ(m)g^{s\mod \varphi(m)}，也就轮不到 s 了。

  • 否则设 pp 为 φ(m)\varphi(m) 的质因子，若 s∣ps|p，则(gs)φ(m)p≡(gφ(m))sp≡1(modm)\left(g^s\right)^{\frac{\varphi(m)}{p}}\equiv \left(g^{\varphi(m)}\right)^{\frac{s}{p}}\equiv 1\pmod m，与二式矛盾。

综上，当且仅当 s∈[1,φ(m)],s⊥φ(m)s \in [1 ,\varphi(m) ],s\bot\varphi(m) 时，gsg^s 是模 mm 的原根。根据定义，模 mm 的原根有 φ(φ(m))\varphi(\varphi(m)) 个。

二. NTT 与 FFT 的转换

我们之所以选择 NTT 而非 FFT ，是因为后者由于使用单位根的关系，精度多多少少受到影响，加上部分数论题目要求取模的限制，找一个膜意义下的 “单位根”替身 未尝不可。

我们盯上了 原根，由于其阶的周期性与单位根的周期性匹配的非常好，所以对于一个质数 p=k⋅n+1p=k\cdot n+1 ，我们可以保证 其膜 p 意义下的原根 gkg^k 与 ωn1\omega_n^1 是等价的。

由于通常我们取 n=2qn=2^q，所以若一个质数减一后含有很多 2 的幂，那他就适合做 NTT 的模数 p。 （ 比如说 998244353=7∗17∗223+1998244353=7*17*2^{23}+1 ）

给出递归版的 NTT（因为非递归满地都是）：

#include <bits/stdc++.h>
#define FO(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
#define cle(x) memset(x,0,sizeof(x))
#define inf 2147483647
#define ll long long
using namespace std;

const ll maxn=2100000,mod=998244353;
const double pi=acos(-1);
ll read()
{  ll 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;
}

ll n,m,f[maxn],g[maxn],y[maxn],inv3;

ll qpow(ll x,ll k)
{	
	ll X=1; for (;k;k>>=1,x=(x*x)%mod) if (k&1) X=(X*x)%mod; return X;
}

void ntt(ll *f,ll l,bool c)
{
	if (l==1) return; 
	ll *L=f,*R=f+l/2;
	for (ll i=0;i<l;i++) y[i]=f[i];
	for (ll i=0;i<(l>>1);i++) L[i]=y[i<<1],R[i]=y[i<<1|1];
	ntt(L,l/2,c); ntt(R,l/2,c);
	ll w0=qpow(c?3:inv3 , (mod-1)/l), w=1;
	for (ll i=0;i<(l>>1);i++) {
		ll t=w*R[i]%mod; w=(w*w0)%mod;
		y[i]=(L[i]+t)%mod, y[i+(l>>1)]=(L[i]-t+mod)%mod; 
		//用 Y，不能用 F 存！会覆盖L,R数组！ 
	}	
	for (ll i=0;i<l;i++) f[i]=y[i];
}

int main()
{
	n=read(), m=read(); inv3=qpow(3,mod-2);
	for (ll i=0;i<=n;i++) f[i]=read();
	for (ll i=0;i<=m;i++) g[i]=read();	

	for (m+=n,n=1;n<=m;n<<=1); 
//	for (ll i=0;i<n;i++) r[i]=(r[i>>1]>>1) | ((i&1)?n>>1:0); //括号括起来！ 
	ntt(f,n,1); ntt(g,n,1);	
	for (ll i=0;i<=n;i++) f[i]=(f[i]*g[i])%mod;
	ntt(f,n,0);  ll invn=qpow(n,mod-2);  
	for (ll i=0;i<=m;i++) printf("%lld ",(f[i]*invn)%mod);
}

这玩意比 FTT 还慢，仅供娱乐，请勿参考。

肆. 多项式儿童套餐

仅仅知道如何求出两个多项式的卷积没啥用，对于大多题目只有知道如何转化才能利用所学的知识解题。

但学的内容又比较肤浅，既没有 多项式计数 也没有 多项式指数对数计算，多项式快速幂 这些高阶算法。算不上全家桶，就当是儿童套餐吧。

一. 生成函数入门

对于一个序列 aia_i ，它的生成函数记做F(x)=∑ai∗xiF(x)=\sum a_i*x^i，注意到 aia_i 可以有限项，也可以无限项。

但不论是哪一种情况，如果我们能搞出生成函数的系数序列，我们也就推出了 aia_i 的通项，这是确定的。

1. 生成函数的四则运算

已知序列 a,b 的生成函数为f(x),g(x)f(x),g(x)，则对于序列 a+b 的生成函数自然是 f(x)+g(x)f(x)+g(x)，减法同理。

对于乘法，我们可以表示为 f(x)∗g(x)=∑nxi⋅∑i=0nai∗bn−if(x)*g(x)=\sum_n x^i \cdot \sum_{i=0}^n a_i*b_{n-i}，注意到 a,ba,b 的下标和始终为 n，这对应了多项式里的卷积形式。同时我们也称序列 ∑i=0nai∗bn−i\sum_{i=0}^n a_i*b_{n-i} 的生成函数为 f(x)∗g(x)f(x)*g(x).

除法暂且不谈。

2. 生成函数的封闭形式

  • 简单的例子（1）： f(x)=∑nxnf(x)=\sum_n x^n ，系数列为 1,1,1..{1,1,1..} 。可构造出等式f(x)∗x+1=f(x)f(x)*x+1=f(x)，解得f(x)=11−xf(x)=\frac{1}{1-x}.

  • 简单的例子（2）： g(x)=∑1n!⋅xng(x)=\sum \frac{1}{n!} \cdot x^n，由泰勒展开可知 g(x)=exg(x) = e^x

  • 简单的例子（3）： g(x)=∑npn⋅xng(x)=\sum_n p^n\cdot x^n，系数列为1,p,p2..1,p,p^2..，可构造px⋅f(x)+1=f(x)px\cdot f(x)+1=f(x) ，解得g(x)=11−pxg(x)=\frac{1}{1-px}.

  • 复杂的例子： g(x)=∑nn∗xng(x)=\sum_n n*x^n，系数列为1,2,3..1,2,3..，由卷积性质可知： g(x)=f(x)∗f(x)=∑nxi⋅∑i=1n1g(x)=f(x)*f(x)=\sum_n x^i \cdot \sum_{i=1}^n 1 ，即g(x)=1(1−x)2g(x)=\frac{1}{(1-x)^2}.

  • 经典的例子（1）：斐波拉契数列的生成函数

已知 a0=0,a1=1 ... an=an−1+an−2a_0=0,a_1=1~...~a_n=a_{n-1}+a_{n-2}，生成函数为 F(x)F(x)，凑项可得： F(x)=x∗F(x)+x2⋅F(x)+(a1−a0)x+a0F(x)=x*F(x)+x^2\cdot F(x)+(a_1-a_0)x+a0

解得：

F(x)=x1−x−x2F(x)=\frac{x}{1-x-x^2}

利用配项法，A1−ax+B1−bx=x1−x−x2\frac{A}{1-ax}+\frac{B}{1-bx}=\frac{x}{1-x-x^2}，底下的 a,b用韦达，上面通分联立求解，可知：

{A=15,B=−15a,b=1±52\begin{cases}A=\frac{1}{\sqrt5},B=-\frac{1}{\sqrt5} \\a,b=\frac{1\pm\sqrt5}{2}\end{cases}

我们还可以注意到，提了一个 15\frac{1}{\sqrt5}后里面满足简单例子（2）的格式，说明其展开式是两个公比为 a,ba,b 等比数列，分别代入整理可解得最终结果为：

F(x)=15∑xn⋅((1+52)n−(1−52)n)F(x)=\frac{1}{\sqrt5} \sum x^n \cdot \left( (\frac{1+\sqrt5}{2})^n-(\frac{1-\sqrt5}{2})^n \right )

但实际并没有什么卵用。。。应该吧。

指数型母函数：

对于一个多重集，其中 a1 重复 n1 次，a2 重复 n2 次，…若从中取出 r 个 排列不同的方案 所对应的母函数即为：

G(x)=∏(1+x1!+x22!+...+xii!)G(x)=\prod (1+\frac{x}{1!}+\frac{x^2}{2!}+...+\frac{x^i}{i!})

每一个括号内都称作该项的 指数型母函数，最终对应的答案为 G(r)∗F[r]G(r)*F[r]。

二.多项式求逆

咕咕咕。

三. 例题：K 阶差分前缀和

题目大意： 给定一序列 aia_i 与 KK ，求其 K 阶差分 or 前缀和序列。

解析：

我们对于 aia_i 给出他的生成函数 G(x)=∑iai∗xiG(x)=\sum_i a_i*x^i ，并尝试构造一个函数 G(x)G(x)，使二者的卷积就是我们答案。

我们容易想到的是，对于 k=1k=1 的前缀和，我们可以构造 G(x)=∑1⋅xiG(x)=\sum 1\cdot x^i，由卷积定义我们可以得到 F(x)∗G(x)=∑xi⋅∑k≤iakF(x)*G(x)=\sum x^i\cdot \sum_{k\leq i} a_k，即其第 i 项对应了 aia_i 的一阶前缀和。

同理，一阶差分要求前后两项相减，不难得到 G(x)=1−xG(x)=1-x 便可满足条件。

那么我们应该如何把它拓展到 k 阶呢？由于卷积的结合律，若要做 k 次前缀和，就相当于 F(x)∗G(x)kF(x)*G(x)^k，也就是 G(x)=1(1−x)k=(1−x)−kG(x)=\frac{1}{(1-x)^k}=(1-x)^{-k}。对于 kk 阶差分，也同样可以得到G(x)=(1−x)kG(x)=(1-x)^k。

而两者的处理方式却又微妙的差别：

  1. 对于差分的情况 运用二项式定理，其系数表示为 bi=(−1)i⋅(ki)b_i=(-1)^i\cdot \binom{k}{i}，根据 n 的范围完全可以线性递推。

  2. 对于前缀和的情况，运用广义二项式定理，其系数表示为

bi=(−1)i⋅−ki‾i !=(−1)i∗(−k)(−k−1)...(−k−n+1)i !b_i=(-1)^i\cdot \frac{-k^{\underline i}}{i~!}=(-1)^i*\frac{(-k)(-k-1)...(-k-n+1)}{i~!}

把 (−1)i(-1)^i 放进去，得到 bi=k(k+1)..(k+n−1)i !=(k+n−1i)b_i=\frac{k(k+1)..(k+n-1)}{i~!}=\binom{k+n-1}{i}，发现也可以线性递推。

然后把 aia_i 和 bib_i 一卷就可以了。

#include <bits/stdc++.h>
#define FO(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
#define cle(x) memset(x,0,sizeof(x))
#define inf 2147483647
#define ll long long
using namespace std;

const ll maxn=2000100,mod=1004535809;

ll read(ll p)
{	ll 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) )%p+(ch^48))%p;
	return f?-x:x;
}

ll L,n,k,t,a[maxn],inv3;
ll r[maxn],f[maxn],inv[maxn];

ll qpow(ll x,ll k){ ll X=1ll; for (;k;k/=2,x=(x*x)%mod) if (k&1) X=(X*x)%mod; return X%mod; }

void ntt(ll *f,bool c)
{
	for (int i=0;i<n;i++) if (i<r[i]) swap(f[i],f[ r[i] ]);
	for (int l=2;l<=n;l<<=1) {  //等于号不能丢，不然不能合并最后一步 
		ll w0=qpow( c?3:inv3 , (mod-1)/l ); 
		for (int i=0;i<n;i+=l) {
			ll w=1;
			for (int k=i;k<(i+l/2);k++,w=(w*w0)%mod) {
				ll t = (w*f[k+l/2])%mod , u = f[k] ;
				f[k] = (u+t)%mod; f[k+(l>>1)] = (u-t+mod)%mod;
			} 
		}
	} 
	if (!c) { ll invn= qpow(n,mod-2); for (int i=0;i<n;i++) f[i]=(f[i]*invn)%mod; }
}

void cal()
{
	L=2*n; for (n=1;n<L;n<<=1); 
	for (int i=1;i<n;i++) r[i]= (r[i>>1]>>1) | ((i&1)?(n>>1):0); 
	ntt(a,1);  ntt(f,1);
	for (int i=0;i<n;i++) a[i] = (a[i]*f[i])%mod;  ntt(a,0); 
}

int main()
{
	n = read(mod); k = read(mod); t = read(mod);  inv3 = qpow(3,mod-2);
	for (ll i=0;i<n;i++) a[i] = read(mod); inv[1]=1ll;
	
	for (ll i=2;i<n;i++) inv[i]= (ll)(mod-mod/i)*inv[mod%i]%mod;
	
	if (t == 1) { f[0]=1; for (ll i=1;i<n;i++) f[i] = ( ((mod-f[i-1])*inv[i])%mod * (k-i+1+mod)%mod )%mod; }
	else { f[0]=1; for (ll i=1;i<n;i++) f[i] = ( (f[i-1]*inv[i])%mod * (k+i-1+mod)%mod )%mod; }
	
	cal(); for (int i=0;i<L/2;i++) printf("%lld ",a[i]);
}

终. 参考资料

Oi WIKI

command_block的NTT全家桶

Karry的多项式题单