多项式

又回来做多项式了？

基础工具

FFT

单独写了博客，在这里。

NTT

这里介绍 NTT。

FFT 容易丢精度，常数有些大，所以我们考虑在模意义下找到一个单位根的替代。

NTT 对模数有特殊要求，若模数 \(p=k\cdot 2^t +1\)，则只能支持 \(n\le2^t\) 的卷积，比 FFT 更快。（引自 九阳哥Nine_Suns的博客）

常用的 NTT 模数（质数）有：\(998244353,1004535809,469762049\)。（引自 阿尼哥ANIG的博客）

具体怎么实现呢？

（下面部分引自 寇门德command_block的博客）

我们取模数的原根（上述三个模数原根都是 \(3\)）。

我们找一个东西能够代替 \(\omega_n^1\)。

发现是这个东西：\(g^{\frac{p-1}{n}}\)。

具体为什么，不会整，反正知道它也满足 \(\omega\) 的几条性质即可。

如果 \(n\) 无法整除 \(p-1\)，则这个东西无法当做单位根使用，这也印证了为什么九阳哥如是说。

因为我们分治的过程中，我们的长度始终是 \(2\) 的幂，所以这几个模数很好！

我们写一下这三个模数减一后的质因数分解：

\[\begin{aligned}&998244353-1=2^{23}\cdot 7\cdot 17\\&1004535809-1=2^{21} \cdot479\\&469762049-1=2^{26}\cdot 7\end{aligned} \]

通过上面的论证，我们已经能体会到这几个模数的优良之处！

界

我们记 \(F(x)\bmod x^n\) 表示只保留 \(F(x)\) 的 \(0\sim n-1\) 次项。

多项式初等函数

多项式牛顿迭代、求导、积分、复合

求导/积分

这部分粘贴自九阳哥的博客。

先来几个常见导数：

\[\begin{aligned}&(x^n)'=nx^{n-1}\\&(\ln x)'=\frac{1}{x}\\&(e^x)'=e^x\\&(a^x)'=a^x\ln x\\&(\log_a x)'=\frac{1}{x\ln a}\end{aligned} \]

运算法则：

\[\begin{aligned}&\Big(f(x)+g(x)\Big)'=f'(x)+g'(x)\\&\Big(f(x)-g(x)\Big)'=f'(x)-g'(x)\\&\Big(f(x)g(x)\Big)'=f(x)g'(x)+f'(x)g(x)\\&\Bigg(f\Big(g(x)\Big)\Bigg)'=f'\Big(g(x)\Big)\cdot g'(x)\end{aligned} \]

前两个是加减法，第三个是乘法，第四个是复合。

对于多项式 \(F(x)=\sum\limits_{i=0}^{n-1}f_ix^i\)，我们定义多项式的求导，积分和复合。

求导：

\[F'(x)=\sum\limits_{i=1}^{n-1}f_iix^{i-1} \]

积分：

\[\sum\limits_{i=1}^{n}\frac{f_{i-1}x^i}{i} \]

复合：

\[F\Big(G(x)\Big)=\sum\limits_{i=0}^{n}f_iG(x)^i \]

只能看懂求导是个啥，积分是个啥不明白，这个 \(C\) 更不明白……

咱手玩一下，发现不管那个 \(C\)，这个求导和积分就是一对逆运算。

也就是对 \(F(x)'\) 积分，能够得到 \(F(x)\)。

代码：

void qiudao(ll *f,ll n){
	for(ll i=1;i<n;i++){
		f[i-1] = f[i] * i %mod;
	}f[n-1]=0;//求导少一次 
}
ll invnum[2*N];
void getinv(int lim){
	invnum[1]=1;
	for (ll i=2;i<=lim;i++){
		invnum[i]=invnum[mod%i]*(mod-mod/i)%mod;
	}
}
void jifen(ll *f,ll n){
	for(int i=n;i>=1;i--){
		f[i] = f[i-1] * invnum[i] %mod;
	}f[0]=0;//积分多一次 
}

泰勒展开

\[F(x)=\sum\limits_{i=0}\frac{F^{(i)}(x_0)}{i!}(x-x_0)^i \]

其中 \(F^{(i)}(x)\) 为 \(F(x)\) 的 \(i\) 阶导数。

理解为在 \(x_0\) 处模拟 \(F(x)\) 的值，值的变化趋势，变化趋势的变化趋势……

多项式牛顿迭代

若 \(f\Big(G_0(x)\Big)\equiv 0\pmod {x^\frac{n}{2}}\)，\(f\Big(G(x)\Big)\equiv 0\pmod {x^n}\)

则有：

\[G(x)\equiv G_0(x)-\frac{f\Big(G_0(x)\Big)}{f'\Big(G_0(x)\Big)}\pmod {x^n} \]

证明：（引自九阳哥博客）

将 \(f\Big(G(x)\Big)\) 在 \(G_0(x)\) 处泰勒展开，有：

\[f\Big(G(x)\Big)=\sum\limits_{i=0}\frac{f^{(i)}\Big(G_0(x)\Big)}{i!}\Big(G(x)-G_0(x)\Big)^i \]

比较显然，有 \(G_0(x)\equiv G(x)\pmod {x^\frac{n}{2}}\)，所以 \(G(x)-G_0(x)\equiv 0\pmod {x^\frac{n}{2}}\)，即 \(G(x)-G_0(x)\) 的前 \(\dfrac{n}{2}\) 项都是 \(0\)。

于是有在 \(k>1\) 时，\(\Big(G(x)-G_0(x)\Big)^k\equiv0\pmod {x^n}\)，即这个式子从 \(2\) 次方开始，它的前 \(n\) 项必定为 \(0\)。于是上式化为：

\[f\Big(G(x)\Big)\equiv f\Big(G_0(x)\Big) + f'\Big(G_0(x)\Big)\Big(G(x)-G_0(x)\Big)\pmod {x^n} \]

由于 \(f\Big(G(x)\Big)\equiv 0\pmod {x^n}\)，拆个括号，有：

\[0\equiv f\Big(G_0(x)\Big) + f'\Big(G_0(x)\Big)G(x)-f'\Big(G_0(x)\Big)G_0(x)\pmod {x^n} \]

移项，有：

\[f'\Big(G_0(x)\Big)G(x)\equiv f'\Big(G_0(x)\Big)G_0(x) - f\Big(G_0(x)\Big)\pmod {x^n} \]

两边同时除以 \(f'\Big(G_0(x)\Big)\)，有：

\[G(x)\equiv G_0(x) - \frac{f\Big(G_0(x)\Big)}{f'\Big(G_0(x)\Big)}\pmod {x^n} \]

证毕。

多项式求逆

P4238 【模板】多项式乘法逆

已知 \(F(x)\)，求一个 \(G(x)\)，使得：

\[F(x)G(x)\equiv 1\pmod {x^n} \]

考虑一个迭代倍增的过程：先求出模 \(x^t\) 意义下的答案，用这个答案求出模 \(x^{2t}\) 意义下的答案，然后令 \(t\leftarrow2t\)，在求下一步，一直倍增到 \(t\ge n\) 结束。（引自九阳哥博客）

当 \(t=1\) 时，即多项式只有常数项。这时候直接求数字逆元即可。

假如我们已求出模 \(x^{\frac{t}{2}}\) 意义下的答案，记为 \(G_0(x)\)，即：

\[G_0(x)F(x)\equiv 1\pmod {x^{\frac{t}{2}}} \]

\[G_0(x)F(x) -1\equiv 0\pmod {x^{\frac{t}{2}}} \]

两边平方，有：（正确性是有的，可以自己找两个例子试一下）

\[(G_0(x)F(x) -1)^2\equiv 0\pmod {x^t} \]

\[G_0(x)^2F(x)^2 -2G_0(x)F(x) +1\equiv 0\pmod {x^t} \]

\[\Big(G_0(x)^2F(x) -2G_0(x)\Big)F(x) \equiv -1\pmod {x^t} \]

\[\Big(2G_0(x)-G_0(x)^2F(x)\Big)F(x) \equiv 1\pmod {x^t} \]

即在模 \(x^t\) 意义下的答案 \(G(x)\)，满足以下递推式：

\[G(x)\equiv2G_0(x)-G_0(x)^2F(x)\pmod{x^t} \]

迭代计算即可。复杂度 \(O(n\log n)\)。

void invp(ll *f,ll n){//f->f^{-1}
	ll m=get2p(n);
	res[0]=qpow(f[0],mod-2);
	for(ll len=2;len<=m;len<<=1){
		for(int i=0;i<len;i++) tmpf[i]=f[i];
		for(int i=0;i<(len>>1);i++) tmpg[i]=(res[i] << 1) %mod;
		ntt(res ,(len<<1),0);ntt(tmpf,(len<<1),0);
		for(int i=0;i<(len<<1);i++) res[i] = (res[i] %mod) * (res[i] %mod) %mod * (tmpf[i]%mod)%mod; 
		ntt(res ,(len<<1),1);
		for(int i=len;i<(len<<1);i++) res[i]=0;
		for(int i=0;i<len;i++) res[i] = (tmpg[i] - res[i]%mod +mod)%mod;
	}
	for(int i=0;i<n;i++) f[i] = res[i];
	for(int i=n;i<2*m;i++) f[i]=0;
	for(int i=0;i<2*m;i++) tmpf[i]=tmpg[i]=res[i]=0;
}

有几个很恶心的点，我不得不提一嘴：

1. 我们一定要把 rev 数组的计算放进 NTT 函数里！！！我们长度不同时 rev 的值也不同！

2. 我们循环中，两个长度为 \(len\) 的多项式相乘，结果是 \(2len\) 的！！！尽管我们只取最低的 \(len\) 项，但我们还是要以长度为 \(2len\) 来 NTT！！

3. 还是上面的点，由于我们只要最低的 \(len\) 项，所以我们一定要清空 res[len]~res[(len<<1)-1]！！！不要忘记清空数组！！！！

多项式开根

P5205 【模板】多项式开根

给出 \(F(x)\)，求一个 \(G(x)\)，满足 \(G(x)^2\equiv F(x)\pmod {x^n}\)。

移个项，有：

\[G(x)^2-F(x)\equiv 0\pmod {x^n} \]

我们使用牛顿迭代推式子，设 \(f\Big(G(x)\Big)=G(x)^2-F(x)\equiv 0\pmod {x^n}\)

有：

\[f'\Big(G(x)\Big) = 2G(x) \]

这里我们是对 \(G(x)\) 求导，由于 \(F(x)\) 已知，我们将其视作常数。

于是套牛迭式子，我们记 \(G_0(x)\) 为模 \(x^{\frac{n}{2}}\) 意义下的答案，于是有：

\[G(x)\equiv G_0(x)-\frac{G_0(x)^2-F(x)}{2G_0(x)}\equiv \frac{G_0(x)^2+F(x)}{2G_0(x)}\pmod {x^n} \]

根据这个倍增即可。

多项式只有常数时，开方的结果是求个二次同余方程就行。

• P5205 【模板】多项式开根 这个板子里常数项为 \(1\)，所以直接把 \(1\) 带进去就行了。

• P5277 【模板】多项式开根（加强版） 这个板子保证里常数项为二次剩余，所以套个奇波拉就行了。

void sqrtp(ll *f,ll *res,ll n){
	ll m=get2p(n);
	static ll tmp[2*N],invt[2*N];
	res[0]=sqrtf0;//常数项为 1 时这里就是 1，否则是 f[0] 的二次剩余。 
	for(ll len=2;len<=m;len<<=1){
		for(int i=0;i<(len>>1);i++) tmp[i] = (res[i] << 1) %mod;
		invp(tmp,invt,len); //这里必须是对 x^len 取模 
		ntt(res,len,0);for(int i=0;i<len;i++) res[i] = res[i] * res[i] %mod;//这里不用 len<<1 是因为 res 此时本身就是一个 len>>1 项的多项式，而且只和自己相乘，出来的结果肯定是 len 项的！
		ntt(res,len,1);
		for(int i=0;i<len;i++) res[i] = (res[i] + f[i]) %mod;
		mul(res,invt,len,len);
		for(int i=0;i<(len<<1);i++) invt[i]=tmp[i]=0;//记得清零 
	}
	for(int i=0;i<m;i++) tmp[i]=invt[i]=0;
	return;
}

多项式黛玉带余除法

P4512 【模板】多项式除法

给定 \(n\) 次多项式 \(F(x)\)，\(m\) 次多项式 \(G(x)\)，求多项式 \(Q(x)\)，\(R(x)\)，满足：

• \(Q(x)\) 次数为 \(n-m\)，\(R(x)\) 次数小于 \(m\)
• \(F(x)=G(x)Q(x)+R(x)\)

即带余除法，算 \(F(x)\div G(x)\)。

我们定义一种运算 \(F^T(x)\)，表示把多项式的系数翻转。具体的，有：

\[F^T(x)=x^nF(x^{-1}) \]

我们手玩一个数据，比如 \(F(x)=1+2x+3x^2\)，于是：

\[F^T(x)=x^2(1+2x^{-1}+3x^{-2})=x^2+2x+3 \]

发现确实是系数进行了翻转。

于是我们推式子：

\[F(x)=G(x)Q(x)+R(x) \]

两边同时把 \(x=x^{-1}\) 带入，有：

\[F(x^{-1})=G(x^{-1})Q(x^{-1})+R(x^{-1}) \]

两边同时乘以 \(x^n\)，有：

\[x^nF(x^{-1})=x^mG(x^{-1})\cdot x^{n-m}Q(x^{-1})+x^{n-m+1}\cdot x^{m-1}R(x^{-1}) \]

写成 \(F^T\) 的形式：

\[F^T(x)=G^T(x)Q^T(x)+x^{n-m+1}R^T(x) \]

两边同时模 \(x^{n-m+1}\)，有：

\[F^T(x)\equiv G^T(x)Q^T(x)\pmod {x^{n-m+1}} \]

于是有：

\[Q^T(x)\equiv F^T(x)G^T(x)^{-1}\pmod {x^{n-m+1}} \]

于是我们先把 \(F\) 与 \(G\) 系数翻转，对 \(G^T\) 求逆，就能算出 \(Q^T\) 了，再翻转一次就能得到 \(Q(x)\)。

因为我们有 \(F(x)=G(x)Q(x)+R(x)\)，所以：

\[R(x)=F(x)-G(x)Q(x) \]

直接带进去就能求 \(R\) 了。

void divp(ll *f,ll *g,ll n,ll m){//这里 n 和 m 是项数！ 
	static ll q[2*N],r[2*N],lim=n-m+1;
	reverse(f,f+n),reverse(g,g+m); 
	for(int i=0;i<lim;i++) q[i] = f[i];
	for(int i=0;i<lim;i++) r[i] = g[i];
	invp(r,lim);mul(q,r,lim,lim);
	reverse(q,q+lim);reverse(f,f+n),reverse(g,g+m);
	mul(g,q,n,n);
	for(int i=0;i<m-1;i++){
		g[i]=(f[i]-g[i]+mod)%mod;//f[i]=f[i]-g[i]*q[i]
	}for(int i=0;i<lim;i++) f[i]=q[i];
	for(int i=lim;i<n;i++) f[i]=0;
}//f=gq+r;f->q,g->r

多项式 \(\ln\)

\[\color{red}\text{我不理解它的实际含义。} \]

P4725 【模板】多项式对数函数（多项式 ln）

有一个多项式 \(\ln\) 的麦克劳林级数定义：

\[\ln F(x)=-\sum\limits_{i=1}^{\infty}\frac{\Big(1-F(x)\Big)^i}{i} \]

这啥？

给定多项式 \(F(x)\)，求一个多项式 \(G(x)\)，满足：

\[G(x)\equiv \ln F(x) \pmod {x^n} \]

我们先全当做不是同余式，虽然我不理解此时的意义是什么，但我们还是求：

\[G(x)=\ln F(x) \]

定义 \(f(x)=\ln x\)，于是有：

\[G(x)=f\Big(F(x)\Big) \]

两边同时对 \(x\) 求导，通过复合的求导运算法则，有：

\[G'(x)=f'\Big(F(x)\Big)\cdot F'(x) \]

我们知道 \(\ln\) 的导数 \(\ln'\) 为：

\[(\ln x)'=\frac{1}{x} \]

因为这里 \(f=\ln\)，于是有：

\[f'\Big(F(x)\Big)=\frac{1}{F(x)} \]

也就是：

\[G'(x)=\frac{F'(x)}{F(x)} \]

于是求导，求逆，卷积，积分，完了。

啊？

void lnp(ll *f,ll n){
	static ll g[2*N];
	for(int i=0;i<n;i++) g[i]=f[i];
	qiudao(f,n);
	invp(g,n);
	mul(f,g,n,n);
	jifen(f,n-1);
	ll m=get2p(n);
	for(int i=0;i<(m<<1ll);i++) g[i]=0;
}

多项式 \(\exp\)

\[\color{red}\text{我不理解它的实际含义。} \]

P4726 【模板】多项式指数函数（多项式 exp）

首先明确一下 \(\exp\) 的含义：

\[\exp(x)=e^x \]

有一个多项式 \(\exp\) 的麦克劳林级数定义：

\[\exp F(x)=\sum\limits_{i=0}^{\infty}\frac{F(x)^i}{i!} \]

这啥？

给定多项式 \(F(x)\)，求一个多项式 \(G(x)\)，满足：

\[G(x)\equiv \exp F(x) \pmod {x^n} \]

我们先全当做不是同余式，虽然我不理解此时的意义是什么，但我们还是求：

\[G(x)=\exp F(x) \]

或者说我们写成更和蔼的形式：

\[G(x)=e^{F(x)} \]

两边同时取 \(\ln\)，有：

\[\ln G(x)=F(x) \]

于是有 \(\ln G(x)-F(x)=0\)，有 \(0\) 了，考虑牛迭！

我们设 \(f\Big(G(x)\Big)=\ln G(x)-F(x)\)，对 \(G(x)\) 求导，有（这里 \(F(x)\) 视作常数）：

\[f'\Big(G(x)\Big)=\frac{1}{G(x)} \]

搞牛迭。我们记 \(G_0(x)\) 为模 \(x^\frac{n}{2}\) 意义下的答案，即：

\[G_0(x)\equiv e^{F(x)}\pmod {x^\frac{n}{2}} \]

我们带入牛迭式子，于是有：

\[\begin{aligned}G(x)&\equiv G_0(x)-\frac{\ln G_0(x)-F(x)}{\frac{1}{G_0(x)}}\\&\equiv G_0(x)-G_0(x)\ln G_0(x)+G_0(x)F(x)\\&\equiv G_0(x)\Big(1-\ln G_0(x)+F(x)\Big) \pmod{x^n}\end{aligned} \]

套模板迭代即可。

void expp(ll *f,ll n){
	static ll m=get2p(n),g[2*N],tmp[2*N];
	g[0]=1;
	for(ll len=2;len<=m;len<<=1){
		for(int i=0;i<(len>>1ll);i++) tmp[i] = g[i];
		lnp(g,len);
		for(int i=0;i<len;i++) g[i] = (f[i]-g[i]+mod)%mod; 
		g[0] = (g[0]+1) %mod;
		mul(g,tmp,len,len);
	}
	for(int i=0;i<n;i++) f[i] = g[i];
	for(int i=0;i<(m<<1ll);i++) tmp[i]=g[i]=0;
}

多项式快速幂

P5245 【模板】多项式快速幂

给定 \(F(x)\)，求多项式 \(G(x)\)，满足：

\[F(x)^k\equiv G(x)\pmod {x^n} \]

我们先不管同余式，当做等式计算。两边同时取 \(\ln\)，有：

\[\ln F(x)^k=\ln G(x) \]

根据 \(\ln\) 运算法则，\(\ln (a^b)=b\ln a\)，有：

\[k\ln F(x)=\ln G(x) \]

两边同时 \(\exp\) 回来，有：

\[e^{k\ln F(x)}= G(x) \]

模一下：

\[e^{k\ln F(x)}\equiv G(x)\pmod {x^n} \]

套板子就行了。

注意两点：

1. 这么做的前提是 \(f_0=1\)，不然需要别的办法处理。
2. 我们要让指数 \(k\) 对 \(p\) 取模，这样做的依据是什么？

参照题解：Light_Poet

若 \(F(x)\) 为 \(n\) 次多项式，\(p\) 为质数，有：

\[F(x)^p\equiv F(x^p)\pmod {p} \]

当我们计算 \(F(x)^p\bmod x^n\) 时，由于大多数情况下都有 \(n<p\)，于是有：

\[F(x)^p\equiv a_0\equiv 1\pmod {x^n} \]

故我们可以先让 \(k\) 对 \(p\) 取模后再计算。

void powp(ll *f,ll n,ll k){
	lnp(f,n);
	for(int i=0;i<n;i++) f[i] = (f[i] * k)%mod;
	expp(f,n);
}

啊啊啊调了好久，鉴定为求逆不清空导致的。

多项式三角函数

P5264 多项式三角函数

其实没啥用，写着玩。

我们有欧拉公式：

\[\mathrm{e}^{ix}=\cos x+\mathrm{i} \sin \]

将 \(-x\) 带入，有：

\[\mathrm{e}^{-ix}=\cos x-\mathrm{i}\sin x \]

两式相减，有：

\[2\mathrm{i}\sin x=\mathrm{e}^{ix}-\mathrm{e}^{-ix} \]

即：

\[\sin x=\frac{\mathrm{e}^{ix}-\mathrm{e}^{-ix}}{2\mathrm{i}} \]

两式相加，有：

\[2\cos x=\mathrm{e}^{ix}+\mathrm{e}^{-ix} \]

即：

\[\cos x=\frac{\mathrm{e}^{ix}+\mathrm{e}^{-ix}}{2} \]

将 \(F(x)\) 带入两式，有：

\[\sin F(x)=\frac{\mathrm{e}^{iF(x)}-\mathrm{e}^{-iF(x)}}{2\mathrm{i}} \]

\[\cos F(x)=\frac{\mathrm{e}^{iF(x)}+\mathrm{e}^{-iF(x)}}{2} \]

求个 \(\exp\) 就行了。

注意，根据 \(i\) 的定义式，在模 \(p\) 意义下有：

\[\mathrm{i}^2\equiv -1\pmod p \]

我们解二次同余方程，易得有两个解：

\[i_1=86583718,i_2=911660635 \]

将两个中任意一个带入就能算。

const int I=86583718;
void sinp(ll *f,ll n){
	static ll g[N];
	for(int i=0;i<n;i++) g[i]=((mod-I) * f[i])%mod;
	for(int i=0;i<n;i++) f[i]=( I * f[i])%mod;
	expp(f,n);expp(g,n);
	for(int i=0;i<n;i++){
		f[i] = (f[i] - g[i] + mod)%mod * qpow(2*I%mod,mod-2) %mod;
	}
}
void cosp(ll *f,ll n){
	static ll g[N];
	for(int i=0;i<n;i++) g[i]=((mod-I) * f[i])%mod;
	for(int i=0;i<n;i++) f[i]=( I * f[i])%mod;
	expp(f,n);expp(g,n);
	for(int i=0;i<n;i++){
		f[i] = (f[i] + g[i])%mod * qpow(2,mod-2) %mod;
	}
}

多项式反三角函数

P5265 多项式反三角函数

贴一篇反三角函数求导的证明

有这么三个式子，不会证：

\[\arcsin x=\int\frac{x'}{\sqrt{1-x^2}}\mathrm{d}x \]

\[\arccos x=-\int\frac{x'}{\sqrt{1-x^2}}\mathrm{d}x \]

\[\arctan x=\int\frac{x'}{1+x^2}\mathrm{d}x \]

把 \(F(x)\) 带入，有：

\[\arcsin F(x)=\int\frac{F'(x)}{\sqrt{1-F(x)^2}}\mathrm{d}x \]

\[\arccos F(x)=-\int\frac{F'(x)}{\sqrt{1-F(x)^2}}\mathrm{d}x \]

\[\arctan F(x)=\int\frac{F'(x)}{1+F(x)^2}\mathrm{d}x \]

套板子即可。

这里需要多项式在常数项不为 \(1\) 时开根，套个奇波拉。

void arcsinp(ll *f,ll n){
	static ll g[N];
	for(int i=0;i<n;i++) g[i] = f[i];
	mul(g,f,n,2*n);
	for(int i=0;i<2*n;i++) g[i] = mod - g[i];
	g[0] = (1 + g[0]) %mod;
	sqrtp(g,2*n);
	for(int i=n;i<2*n;i++) g[i]=0;
	invp(g,n);
	qiudao(f,n);
	mul(f,g,n,n);
	jifen(f,n);
	ll m=get2p(n);
    for(int i=0;i<m;i++) g[i]=0;
}
void arccosp(ll *f,ll n){
	arcsinp(f,n);
	for(int i=0;i<n;i++) f[i] = (mod-f[i])%mod;
	ll m=get2p(n);
    for(int i=0;i<m;i++) g[i]=0;
}
void arctanp(ll *f,ll n){
	static ll g[N];
	for(int i=0;i<n;i++) g[i] = f[i];
	mul(g,f,n,2*n);
	g[0] = (1 + g[0]) %mod;
	invp(g,2*n);
	qiudao(f,2*n);
	mul(f,g,n*2,n*2);
	jifen(f,n*2);
	ll m=get2p(n);
    for(int i=0;i<m;i++) g[i]=0;
}

---

最终融合板子

我们放一下最终的板子：

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define fr(a) freopen(a,"r",stdin)
#define fw(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=8e5+10,mod=998244353;
ll qpow(ll a,ll n){
    ll res=1;
    while(n){
        if(n&1) (res *= a) %= mod;
        (a *= a) %= mod,n>>=1;
    }return res;
}
const int G=3,invG=qpow(G,mod-2);
ll rev[N];
ll get2p(ll n){
	ll tmp=1;while((1ll<<tmp) < n) tmp++;
	return 1ll<<tmp;
}
void ntt(ll *f,ll n,bool flag){
	for(int i=0;i<n;i++) rev[i] = (rev[i>>1]>>1) | (i&1?(n>>1):0);
	for(int i=0;i<n;i++) if(i < rev[i]) swap(f[i],f[rev[i]]);
	for(int len=2;len<=n;len<<=1){
		ll tmp=qpow((flag==0)?G:invG,(mod-1)/len);
		for(int k=0;k<n;k+=len){
			ll now=1;
			for(int i=0;i<(len>>1);i++){
				ll l=i+k,r=i+k+(len>>1),tt=f[r] * now %mod;
				f[r] = f[l] - tt;if(f[r] < 0) f[r] += mod;
				f[l] = f[l] + tt;if(f[l]>mod) f[l] -= mod;
				now = now * tmp %mod;
			}
		}
	}if(flag){
		ll invl=qpow(n,mod-2);
		for(int i=0;i<n;i++) f[i] = f[i] * invl %mod;
	}return;
}
void mul(ll *f,ll *g,ll n,ll lim){
	static ll tmp[N];ll m=get2p(2*n);
	for(int i=0;i<m;i++) tmp[i] = g[i];
	ntt(f,m,0),ntt(tmp,m,0);
	for(int i=0;i<m;i++) f[i] = f[i] * tmp[i] %mod;
	ntt(f,m,1);
	for(int i=lim;i<=(m<<1);i++) f[i] = 0;
	for(int i=0;i<=(m<<1);i++) tmp[i] = 0;
}
void invp(ll *f,ll n){
	static ll tmpf[N],tmpg[N],G0[N];
	ll m=get2p(n);G0[0] = qpow(f[0],mod-2);
	for(int len=2;len <= m;len<<=1){
		for(int i=0;i<len;i++) tmpf[i] = f[i];
		for(int i=0;i<(len>>1);i++) tmpg[i] = (G0[i]<<1) %mod;
		ntt(G0,len<<1,0),ntt(tmpf,len<<1,0);
		for(int i=0;i<(len<<1);i++) G0[i] = G0[i] * G0[i] %mod * tmpf[i] %mod;
		ntt(G0,len<<1,1);
		for(int i=len;i<(len<<1);i++) G0[i] = 0;
		for(int i=0;i<len;i++) G0[i] = (tmpg[i] - G0[i] + mod) %mod;
	}for(int i=0;i<n;i++) f[i] = G0[i];
	for(int i=0;i<(m<<1);i++) G0[i] = tmpf[i] = tmpg[i] = 0;
}
/*----------------奇波拉----------------*/ 
ll w;
struct cplx{
    ll x,y;
    cplx (ll xin=0,ll yin=0){x=xin,y=yin;}
    void MO(){
        x=(x%mod+mod)%mod;
        y=(y%mod+mod)%mod;
    }
    friend cplx operator*(cplx a,cplx b){
        return cplx(a.x*b.x%mod+a.y*b.y%mod*w%mod,a.x*b.y%mod+a.y*b.x%mod);
    }
};
cplx cqpow(cplx a,ll n){
    cplx res(1,0);
    while(n){
        if(n&1) res=res*a;res.MO();
        a=a*a,a.MO(),n>>=1;
    }return res;
}
ll Cipolla(ll n){
    srand(time(0));
    if(!n) return 0;
    ll a;
    while(1){
        a=rand()%mod;
        if(qpow((a*a-n)%mod,(mod-1)/2)==mod-1) break;
    }
    w=((a*a-n)%mod+mod)%mod;
    cplx tmp=cplx(a,1) * cplx(a,1);
    cplx ans=cqpow(cplx(a,1),(mod+1)/2);ans.MO();
    ll res1=ans.x,res2=mod-ans.x;
    return min(res1,res2);
}
/*----------------奇波拉----------------*/ 
void sqrtp(ll *f,ll n){
	static ll G0[N],invg[N];
	ll m=get2p(n);G0[0] = Cipolla(f[0]);
	for(int len=2;len <= m;len<<=1){
		for(int i=0;i<(len>>1);i++) invg[i] = (G0[i]<<1) %mod;
		ntt(G0,len,0);
		for(int i=0;i<len;i++) G0[i] = G0[i] * G0[i] %mod;
		ntt(G0,len,1);
		for(int i=0;i<len;i++) G0[i] = (G0[i] + f[i]) %mod;
		invp(invg,len);mul(G0,invg,len,len);
		for(int i=0;i<len;i++) invg[i] = 0;
	}for(int i=0;i<n;i++) f[i] = G0[i];
	for(int i=0;i<(m<<1);i++) G0[i] = invg[i] = 0;
}
void qiudao(ll *f,ll n){
	for(int i=1;i<n;i++){
		f[i-1] = f[i] * i %mod;
	}f[n-1] = 0;return;
}
ll inv[N];
void getinv(ll lim){
	inv[1] = 1;
	for(int i=2;i<=lim;i++){
		inv[i]=inv[mod%i] * (mod - mod / i) %mod;
	}return;
}
void jifen(ll *f,ll n){
	for(int i=n;i>=1;i--){
		f[i] = inv[i] * f[i-1] %mod;
	}f[0] = 0;return;
}
void lnp(ll *f,ll n){
	static ll tmpf[N];ll m=get2p(2*n);
	for(int i=0;i<n;i++) tmpf[i] = f[i];
	qiudao(tmpf,n);invp(f,n);mul(f,tmpf,n,n);jifen(f,n-1);
	for(int i=0;i<(m<<1);i++) tmpf[i] = 0;return;
}
void expp(ll *f,ll n){
	static ll G0[2*N],tmpg[2*N];
	ll m=get2p(n);G0[0] = 1;
	for(int len=2;len<=m;len<<=1){
		for(int i=0;i<(len>>1);i++) tmpg[i] = G0[i];
		lnp(tmpg,len);
		for(int i=0;i<len;i++) tmpg[i] = (f[i] - tmpg[i] + mod) %mod;
		tmpg[0] = (1+tmpg[0]) %mod;
		mul(G0,tmpg,len,len);
	}for(int i=0;i<n;i++) f[i] = G0[i];
	for(int i=0;i<(m<<1);i++) G0[i] = tmpg[i] = 0;
}
void powp(ll *f,ll n,ll k){
    lnp(f,n);
    for(int i=0;i<n;i++) f[i] = (f[i] * k)%mod;
    expp(f,n);
}
void divp(ll *f,ll *g,ll n,ll m){
	if(n<m){
		for(int i=0;i<m;i++) g[i] = 0;
		for(int i=0;i<n;i++) g[i] = f[i] , f[i] = 0;
		return;
	}static ll q[N],r[N];
	ll lim=n-m+1,tmp=get2p(2*n);
	reverse(f,f+n),reverse(g,g+m);
	for(int i=0;i<lim;i++) q[i] = f[i],r[i] = g[i];
	invp(r,lim),mul(q,r,lim,lim);reverse(q,q+lim);//算 q
	reverse(g,g+m),reverse(f,f+n);mul(g,q,n,n);
	for(int i=0;i<n;i++) g[i] = (f[i] - g[i] + mod) %mod;
	for(int i=0;i<lim;i++) f[i] = q[i];
	for(int i=lim;i<=tmp;i++) f[i] = 0;
	for(int i=m;i<=tmp;i++) g[i] = 0;
	for(int i=0;i<=tmp;i++) q[i] = r[i] = 0;
}//f=gq+r;f->q,g->r;
int main(){
    getinv(N-10);
    return 0;
}

---

牛逼工具

任意模数NTT/FTT

P4245 【模板】任意模数多项式乘法

有时候我们会碰见毒瘤模数。

我们有两种处理方式：三模NTT/拆系数FFT

三模NTT

代码是对 ANIG 的拙劣模仿。

我们把前面提到的 NTT 三个模数，用作这里的三个模数。

我们分别计算出模这三个模数意义下的答案，然后用中国剩余定理合并。

具体的，我们有：

\[\begin{aligned}x\equiv a\pmod {p_1}\\x\equiv b\pmod {p_2}\\x\equiv c\pmod {p_3}\end{aligned} \]

我们解得一个最小的正整数解 \(x\)，这个东西再模上题目给的模数 \(p\) 即可。我们解这个 \(x\)：

设 \(x\equiv k_1p_1+a\pmod {p_1}\)，有：

\[k_1p_1+a\equiv b\pmod {p_2} \]

\[k_1p_1\equiv b-a\pmod {p_2} \]

\[k_1\equiv (b-a)p_1^{-1}\pmod {p_2} \]

于是我们解出 \(k_1\)。此时我们令 \(x_1=k_1p_1+a\)，于是我们设 \(x=k_2p_1p_2+x_1\)

有：

\[k_2p_1p_2+x_1\equiv c\pmod {p_3} \]

\[k_2p_1p_2\equiv c-x_1\pmod {p_3} \]

\[k_2\equiv (c-x_1) (p_1p_2)^{-1}\pmod {p_3} \]

这样我们就能解出来 \(k_2\)。此时方程的最小整数解就是 \(x=k_2p_1p_2+x_1\)。

//对 ANIG 的拙劣模仿 
//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define frr(a) freopen(a,"r",stdin)
#define fww(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=2e6+10;
const int mod1=998244353,mod2=1004535809,mod3=469762049,G=3;
ll qpow(ll a,ll n,ll mod){
	n%=mod-1,n+=mod-1,n%=mod-1;
	a%=mod;
	ll res=1;
	while(n){
		if(n&1) (res *= a) %= mod;
		(a *= a) %= mod,n>>=1;
	}return res;
}
const int inv1=qpow(mod1,mod2-2,mod2),inv2=qpow(1ll * mod1 * mod2 ,mod3-2,mod3);
ll mod;
struct node{
	ll a,b,c,d;
	node (ll aa=0,ll bb=0,ll cc=0,ll dd=0){a=aa,b=bb,c=cc,d=dd;}
	friend node operator+(node x,node y){return node((x.a+y.a)%mod1,(x.b+y.b)%mod2,(x.c+y.c)%mod3,(x.d+y.d)%mod);}
	friend node operator-(node x,node y){return node((x.a-y.a)%mod1,(x.b-y.b)%mod2,(x.c-y.c)%mod3,(x.d-y.d)%mod);}
	friend node operator*(node x,node y){return node((x.a*y.a)%mod1,(x.b*y.b)%mod2,(x.c*y.c)%mod3,(x.d*y.d)%mod);}
	friend node operator*(node x,ll t){return node((x.a*t)%mod1,(x.b*t)%mod2,(x.c*t)%mod3,(x.d*t)%mod);}
	void CRT(){
		a=(a+mod1)%mod1,b=(b+mod2)%mod2,c=(c+mod3)%mod3;
		ll k1=(b-a) * inv1 %mod2 , x = k1 * mod1 + a , k2=(c - x%mod3 +mod3)%mod3*inv2%mod3;
		d=k2%mod * mod1%mod * mod2%mod + x%mod;
		d=(d%mod+mod)%mod;
	}
}f[N],g[N],res[N];
ll n,m,rev[N];
ll get2p(ll n){
	ll tmp=1;
	while((1ll<<tmp) < n) tmp++; tmp=1ll<<tmp;
	return tmp;
}
node pows(ll a,ll b,ll c,ll d,ll n){
	return node(qpow(a,-1,mod1),qpow(b,-1,mod2),qpow(c,-1,mod3),1145141919810ll);
}
void ntt(node *f,ll n,bool flag){
	for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1) | (i&1?(n>>1):0);
	for(int i=0;i<n;i++) if(i<rev[i]) swap(f[i],f[rev[i]]);
	ll sgn=(flag==0)?1:-1;
	for(int len=2;len<=n;len<<=1){
		node tmp;
		tmp.a=qpow(3,sgn*(mod1-1)/len,mod1);
		tmp.b=qpow(3,sgn*(mod2-1)/len,mod2);
		tmp.c=qpow(3,sgn*(mod3-1)/len,mod3);
		tmp.d=0;
		for(int k=0;k<n;k+=len){
			node now=node(1,1,1,1);
			for(int i=0;i<(len>>1);i++){
				node tt=now*f[k+i+(len>>1)];
				f[k+i+(len>>1)]=f[k+i] - tt; 
				f[k+i]=f[k+i] + tt;
				now = now*tmp;
			}
		}
	}return;
}
int main(){
	scanf("%lld%lld%lld",&n,&m,&mod);n++,m++;
	ll tmp;
	for(int i=0;i<n;i++) scanf("%lld",&tmp),f[i]=node(tmp,tmp,tmp,0);
	for(int i=0;i<m;i++) scanf("%lld",&tmp),g[i]=node(tmp,tmp,tmp,0);
	tmp=get2p(n+m);
	ntt(f,tmp,0),ntt(g,tmp,0);
	for(int i=0;i<tmp;i++) res[i] = f[i] * g[i];
	ntt(res,tmp,1);
	node inv=pows(tmp,tmp,tmp,tmp,-1);
	for(int i=0;i<tmp;i++) res[i] = res[i] * inv , res[i].CRT();
	for(int i=0;i<n+m-1;i++) printf("%lld ",(res[i].d%mod+mod)%mod);
	return 0;
}

拆系数 FFT（又称MTT）

融合了 九阳哥 的代码和 Prean 大佬的思路。

我们取 \(c=2^{15}\)。

我们把题面描述一下：

求给定两个多项式，求 \(F(x)\cdot G(x)\)，系数对 \(p\) 取模。

我们设 \(F(x)=A(x)+c\cdot B(x)\)，\(G(x)=C(x)+c\cdot D(x)\)，其中 \(B,D\) 项的系数均小于 \(c\)。

于是我们有：

\[\begin{aligned}&F(x)\cdot G(x)\\=&\Big(A(x)+c\cdot B(x)\Big)\Big(C(x)+c\cdot D(x)\Big)\\=&A(x)\cdot C(x)+c\Big(A(x)\cdot D(x) + B(x)\cdot C(x)\Big) + c^2\cdot B(x)\cdot D(x)\end{aligned} \]

发现我们一共需要四个东西：

\[\begin{aligned}A(x)\cdot C(x)\;,\;A(x)\cdot D(x)\\B(x)\cdot C(x)\;,\;B(x)\cdot D(x)\end{aligned} \]

于是我们设 \(T(x)=C(x)+iD(x)\)，求出 \(A(x)\cdot T(x)\) 与 \(B(x)\cdot T(x)\)，有：

\[\begin{aligned}A(x)\cdot T(x)=A(x)\cdot C(x)+iA(x)\cdot D(x)\\B(x)\cdot T(x)=B(x)\cdot C(x)+iB(x)\cdot D(x)\end{aligned} \]

我们要的四个东西都出来了！正好在 \(A\cdot T\) 和 \(B\cdot T\) 的实部和虚部上。

于是我们分别对 \(A,B,T\) 进行 DFT，对 \(A\cdot T,B\cdot T\) 进行 IDFT，一合并就行了。合并部分代码不会写，遂抄自九阳哥。

九阳哥，你真是我的好大哥啊啊啊啊啊啊啊呜呜呜

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define frr(a) freopen(a,"r",stdin)
#define fww(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
#define double long double
using namespace std;
typedef long long ll;
const int N=2e6+10;
const double Pi=acos(-1);
ll n,m,mod,tr[2*N];
ll f[N],g[N],res[N];
struct cplx{
	cplx (double xin=0,double yin=0){x=xin,y=yin;}
	double x,y;
	friend cplx operator+(cplx a,cplx b){return cplx(a.x+b.x,a.y+b.y);}
	friend cplx operator-(cplx a,cplx b){return cplx(a.x-b.x,a.y-b.y);}
	friend cplx operator*(cplx a,cplx b){return cplx(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}A[N],B[N],t[N],at[N],bt[N];
ll get2p(ll n){
	ll tmp=1;
	while((1ll<<tmp) < n) tmp++; tmp=1ll<<tmp;
	return tmp;
}
void fft(cplx *f,ll n,bool flag){//第三版·迭代实现 
	for(int i=0;i<n;i++) tr[i]=(tr[i>>1]>>1) | (i&1?(n>>1):0);
	for(int i=0;i<n;i++) if(i<tr[i]) swap(f[i],f[tr[i]]);
	for(ll len=2;len<=n;len<<=1){
		cplx tmp(cos(2*Pi/len),(flag==0?1.0:-1.0)*sin(2*Pi/len));
		for(int k=0;k<n;k+=len){
			cplx now(1,0);
			for(int i=0;i<(len>>1);i++){
				cplx tt=now*f[k+i+(len>>1)];
				f[k+i+(len>>1)]=f[k+i] - tt;
				f[k+i]=f[k+i] + tt;
				now = now*tmp;
			}
		}
	}if(flag){
		for(int i=0;i<n;i++){
			f[i] = cplx(f[i].x/(double)n,f[i].y/(double)n); 
		}
	}
	return;
}
void mtt(ll *f,ll *g,ll *res,ll n){
	for(int i=0;i<n;i++) {
		A[i].x=f[i]>>15;A[i].y=0;
		B[i].x=f[i]&32767;B[i].y=0;
		t[i].x=g[i]>>15;
		t[i].y=g[i]&32767;
	}
	fft(A,n,0),fft(B,n,0),fft(t,n,0);
	for(ll i=0;i<n;i++) at[i]=A[i] * t[i],bt[i]=B[i] * t[i];
	fft(at,n,1),fft(bt,n,1);
	for (ll i = 0;i < n;i++) {
		ll a = (long long)floor(at[i].x+0.5)%mod;
		ll b = (long long)floor(at[i].y+0.5)%mod;
		ll c = (long long)floor(bt[i].x+0.5)%mod;
		ll d = (long long)floor(bt[i].y+0.5)%mod; 
		A[i]=B[i]=t[i]=at[i]=bt[i]=cplx();
		res[i] = (1ll * (1ll<<30)%mod*(long long)a%mod+(1ll<<15)%mod*((long long)b+c)%mod+(long long)d)%mod; 
		res[i] = (res[i]+mod)%mod; 
	}
} // a + (b+c) * 1<<15 + d*1<<30
int main(){
	scanf("%lld%lld%lld",&n,&m,&mod);n++,m++;
	for(int i=0;i<n;i++) scanf("%lld",&f[i]);
	for(int i=0;i<m;i++) scanf("%lld",&g[i]);
	ll tmp=get2p(n+m);
	mtt(f,g,res,tmp);
	for(int i=0;i<n+m-1;i++) printf("%lld ",res[i]);
	return 0;
}

注：这个题要用 long double。

分治 FFT/NTT 或 半在线卷积

原版

参考博客：九阳哥，command_block

P4721 【模板】分治 FFT

给定数组 \(F\)，求数组 \(G\)，使得：

\[G_0=1 \]

\[G_i=\sum\limits_{j=0}^{i}G_jF_{i-j} \]

保证 \(F_0=0\)。

我们使用一个强制中序遍历的 cdq 分治。

考虑我们当前在处理区间 \([l,r)\)，我们令 \(mid=\dfrac{l+r}{2}\)。

我们先向左递归处理 \([l,mid)\)，处理后 \(G\) 的 \([l,mid)\) 位置的值就算好了！

我们考虑把 \(G\) 的 \([l,mid)\) 向 \([mid,r)\) 提供贡献。

设 \(len=r-l\)，即区间长度。

我们将 \(F\) 的 \([0,len)\) 和 \(G\) 的 \([l,mid)\) 进行卷积（\(G\) 的不够位置补 \(0\)），将卷出来的东西的 \([\dfrac{len}{2},len)\) 部分累加到 \(G\) 的 \([mid,r)\) 位置，就算好了贡献。

区间长度为 \(1\) 时，算出来的就是 \(G\) 的真实值。

这样跑 cdq 就行。

求多项式 \(\exp\)

P4726 【模板】多项式指数函数（多项式 exp）

给定多项式 \(F(x)\)，求 \(G(x)\)，使得：

\[G(x)\equiv e^{F(x)}\pmod {x^n} \]

不考虑取模，我们推式子。

首先，把 \(e^{F(x)}\) 写成 \(\exp()\) 函数的形式，有：

\[G(x)=\exp\Big(F(x)\Big) \]

两边对 \(x\) 求导，右边根据复合函数求导法则，有：

\[G'(x)=\exp'\Big(F(x)\Big)\cdot F'(x) \]

由于 \(\exp'=\exp\)，有：

\[G'(x)=G(x)\cdot F'(x) \]

我们提取第 \(n\) 项的系数，有：

\[(n+1)g_{n+1}=\sum\limits_{i=0}^{n}f_{i+1}(i+1)g_{n-i} \]

右边更多的是 \(i+1\)，于是我们枚举 \(i+1\)，有：

\[(n+1)g_{n+1}=\sum\limits_{i=1}^{n+1}f_iig_{n+1-i} \]

我们令 \(n\) 替换 \(n+1\)，有：

\[ng_n=\sum\limits_{i=1}^{n}f_iig_{n-i} \]

将 \(n\) 除过去，有：

\[g_n=\frac{1}{n}\sum\limits_{i=1}^{n}f_iig_{n-i} \]

我们让另一个多项式系数为 \(f_ii\)，这样右边就变成了我们刚讨论过的问题，\(\dfrac{1}{n}\) 可以在分支过程中区间长度为 \(1\) 时乘上去。

void cdq(ll l,ll r,ll *f,ll *g){
	if(r-l==1) return (g[l] *= qpow(max(l,1ll),mod-2)) %=mod,void(0);
	ll mid=(r+l)>>1,len=r-l;
	cdq(l,mid,f,g);//[l,mid) 的 g 值已经算出来了。
	static ll tf[N],tg[N];
	for(int i=0;i<len;i++) tf[i] = f[i] , tg[i] = 0;
	for(int i=l;i<mid;i++) tg[i-l] = g[i];
	mul(tf,tg,len,len);
	for(int i=mid;i<r;i++) g[i] = (g[i] + tf[(len>>1) + (i-mid)])%mod;
	for(int i=0;i<len;i++) tf[i] = tg[i] = 0;
	cdq(mid,r,f,g);
}
void expp(ll *f,ll n){
	for(int i=0;i<n;i++) f[i] = f[i] * i %mod;
	static ll g[N];
	ll m=get2p(n);g[0]=1;
	cdq(0,m,f,g);
	for(int i=0;i<n;i++) f[i] = g[i];
	for(int i=0;i<(m<<1);i++) g[i] = 0;
}

FWT 快速沃尔什变换

P4717 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)

还得是 command_block 会讲！（按位拆分好抽象啊……）

位运算卷积

给出 \(A,B\) 两个下标范围为 \([0,n-1]\) 的数组（我们默认 \(n\) 是 \(2\) 的整次幂），求数组 \(C\)：

\[C_k=\sum\limits_{i\oplus j=k}A_i\cdot B_j \]

其中 \(\oplus\) 为某种位运算。

模仿

FFT 中，我们利用单位根的性质把两个系数表达式转化成了点值表达式，可以做到 \(O(n)\) 乘起来。这里我们类似地构造。

记 \({\rm FWT}(A)\) 为 \(A\) 经过 FWT 后的某种构造，一共有 \(n\) 项。

我们希望有：

\[\forall i\in[0,n),{\rm FWT}(C)_i={\rm FWT}(A)_i\cdot {\rm FWT}(B)_i\\ 简记作 {\rm FWT}(C)={\rm FWT}(A)\cdot {\rm FWT}(B) \]

相当于一个转点值的操作。这之后再经过一个 IFWT 就能把 \({\rm FWT}(C)\) 转回去了。

我们考虑构造一个 \(\rm{FWT}\) 变换。设 \(c(i,j)\) 为变换系数，即 \(A_j\) 对 \({\rm FWT}(A)_j\) 的贡献。我们有：

\[{\rm FWT}(A)_i=\sum\limits_{j=0}^{n-1}c(i,j)A_j \]

当然，\(c\) 里面扔一个自然数有些丑和冗余，我们根据位运算的独立性，分位考虑。

假设我们构造出 \(c(0或1,0或1)\) 的四种取值，我们就有：

\[c(i,j)=\prod_{k=0}c(i_k,j_k) \]

其中 \(i_k\) 为 \(i\) 二进制下的从高往低的第 \(k\) 位。这很棒。我们待会就要讨论 \(c\) 的四种取值的构造。

接下来我们考虑如何模仿 FFT，用分治和 \(c\) 来计算 FWT！

\[{\rm FWT}(A)_i=\sum\limits_{j=0}^{n-1}c(i,j)A_j \]

把这个式子对半分开。

\[{\rm FWT}(A)_i=\sum\limits_{j=0}^{(n/2)-1}c(i,j)A_j + \sum\limits_{j=(n/2)}^{n-1}c(i,j)A_j \]

设 \(i'\) 为 \(i\) 去掉首位后剩余的数（如果 \(n=8\)，则 \(i\) 的二进制一律按照三位来算），我们把首位拆开。

\[{\rm FWT}(A)_i=\sum\limits_{j=0}^{(n/2)-1}c(i_0,j_0)c(i',j')A_j + \sum\limits_{j=(n/2)}^{n-1}c(i_0,j_0)c(i',j')A_j \]

由于 \(n\) 为二的整次幂，我们会发现前面一个求和中，\(j\) 的最高位（也就是 \(j_0\)）一直为 \(0\)！而后面一个求和中最高位一直为 \(1\)！这样我们能提取公因式：

\[\begin{aligned} {\rm FWT}(A)_i&=c(i_0,0)\sum\limits_{j=0}^{(n/2)-1}c(i',j')A_j +c(i_0,1)\sum\limits_{j=(n/2)}^{n-1}c(i',j')A_j\\ &=c(i_0,0){\rm FWT}(A_0)_i +c(i_0,1){\rm FWT}(A_1)_i \end{aligned}\]

其中 \(A_0\) 是 \(A\) 的 \(0\sim \dfrac{n}{2}-1\) 部分组成的子序列，也就是 \(A\) 的前半段！\(A_1\) 为后半段。

很棒！我们已经分治出子问题了！接下来我们讨论 \(i_0\) 的取值，把整体的分支过程写出来。

当 \(i_0=0\) 时，对应着 \({\rm FWT}(A)_{0\sim(n/2)-1}\)：

\[\forall i\in [0,\frac{n}{2}),{\rm FWT}(A)_i=c(0,0){\rm FWT}(A_0)_i+c(0,1){\rm FWT}(A_1)_i \]

当 \(i_0=1\) 时，对应着 \({\rm FWT}(A)_{(n/2)\sim n-1}\)：

\[\forall i\in [0,\frac{n}{2}),{\rm FWT}(A)_{(n/2)+i}=c(1,0){\rm FWT}(A_0)_i+c(1,1){\rm FWT}(A_1)_i \]

诶，这不是我们 FFT 嘛！只不过 FFT 是乘了个 \(\omega_{n}^{i}\)，这里是乘了个 \(c\) 的某一位！

而且这个甚至不用蝴蝶变化，因为它只是把序列的前半段后半段分开，而不是像 FFT 一样把奇数的放一起，偶数的放一起！

就很棒！

于是乎，我们现在的目的就是要构造出 \(c(i,j)\) 在 \(i,j\in [0,1]\) 时的四种取值，扔进去就行了！

我们之后就把 \(c=\begin{bmatrix}c(0,0) &c(0,1)\\c(1,0) & c(1,1)\end{bmatrix}\) 这个矩阵搞出来就好了！这个矩阵被称之为位矩阵。

另外，求 FWT 的逆变换 IFWT 时，只需要把 \(c\) 矩阵的逆扔进去就行了。为什么之后解释？

位矩阵

我们借助这么一个推论来推理：

通过 \({\rm FWT}(C)={\rm FWT}(A)\cdot {\rm FWT}(B)\)，把 \(\rm{FWT}\) 的定义式带入，可推出：

\[\sum\limits_{j=0}^{n-1}\sum\limits_{k=0}^{n-1}c(i,j)c(i,k)A_jB_k=\sum\limits_{p=0}^{n-1}c(i,p)C_p \]

我们又有：

\[C_k=\sum\limits_{i\oplus j=k}A_i\cdot B_j \]

带入上式一通化简，我们就有：

\[\sum\limits_{j=0}^{n-1}\sum\limits_{k=0}^{n-1}c(i,j)c(i,k)A_jB_k=\sum\limits_{j=0}^{n-1}\sum\limits_{k=0}^{n-1}c(i,j\oplus k)A_jB_k \]

我们就有了：

\[c(i,j)c(i,k)=c(i,j\oplus k) \]

剩余的推导过程自己手玩去吧！command 大哥的博客写的很详细了！这里直接摆结论！

• 或运算

\[c=\begin{bmatrix}1 & 0\\1 & 1\end{bmatrix},c^{-1}=\begin{bmatrix}1 & 0\\-1 & 1\end{bmatrix} \]

• 与运算

\[c=\begin{bmatrix}1 & 1\\0 & 1\end{bmatrix},c^{-1}=\begin{bmatrix}1 & -1\\0 & 1\end{bmatrix} \]

• 异或运算

\[c=\begin{bmatrix}1 & 1\\1 & -1\end{bmatrix},c^{-1}=\begin{bmatrix}0.5 &0.5\\0.5 & -0.5\end{bmatrix} \]

提交记录

子集卷积

P6097 【模板】子集卷积

给定长度为 \(2^n\) 的数组 \(a,b\)，求：

\[c_k=\sum\limits_{i|j=k,i\&j=0}a_ib_j \]

转换条件。记 \(|a|\) 为 \(a\) 二进制下 \(1\) 的个数。

我们有 \(i\&j=0\Leftrightarrow |i|+|j|=|k|\)

我们考虑扩展一维。记 \(A_{len,i}=\begin{cases}a_i & len=|i|\\0&\text{otherwise}\end{cases}\)，\(B,C\) 同理。

我们有：

\[C_p=\sum\limits_{len=0}^{p}A_{len}*B_{p-len} \]

其中 \(*\) 为或的位运算卷积。

什么意思呢？差不多就是把 \(p\) 个 \(1\) 中 \(len\) 个分给 \(A\)，剩余的分给 \(B\)，然后只把 \(A\) 内部满足 \(1\) 的个数为 \(len\) 的项，和 \(B\) 内部满足 \(1\) 的个数为 \(p-len\) 的项卷起来。

答案就是 \(C_{|i|,i}\)。

提交记录

模板题目-CF914G Sum the Fibonacci

---

不知道有什么用

多项式多点求值

P5050 【模板】多项式多点求值

oiwiki

给定一个多项式 \(f(x)\) 和 \(n\) 个点 \(x_1,x_2\dots x_n\)，求：

\[f(x_1),f(x_2)\dots f(x_n) \]

首先，在这个题中，我们可以把多项式项数和点数补到 \(2\) 的整次幂，方便计算。

我们考虑分治解决问题。假设当前在求 \(F(x_1\dots x_n)\) 的值。

构造一个多项式：

\[G(x)=\prod\limits_{i=1}^{\frac{n}{2}}(x-x_i) \]

这个多项式满足 \(G(x_1),G(x_2)\dots G(x_{\frac{n}{2}})\) 都为 \(0\)。

我们求出一个 \(F(x)\bmod G(x)=F_0(x)\)，它有什么性质呢？

我们设 \(F(x)=Q(x)G(x) + F_0(x)\)

将 \(x=x_1,x=x_2\dots x=x_{\frac{n}{2}}\) 带入其中，会发现总是有：

\[\forall i\in [1,\frac{n}{2}],F(x_i)=F_0(x_i) \]

于是乎，我们就能用 \(F_0(x)\) 来代替 \(F(x)\)，求出 \(F_0(x_1\dots x_\frac{n}{2})\) 的值与 \(F(x_1\dots x_\frac{n}{2})\) 相同！

而 \(F_0(x)\) 的规模是和 \(F(x)\) 的一半！求出 \(F_0(x_1\dots x_\frac{n}{2})\) 的过程就是一个子问题。

递归边界：当 \(l=r\) 时，可以验证 \(F(x)\) 是个常数，则当 \(x=x_l\) 时，我们索求的答案正是这个常数。

这样就可以递归算了。

顺便提一嘴，求 \(G(x)\) 可以空间复杂度 \(O(n\log n)\) 预处理。

求 \(G(x)\) 和递归计算答案的时间复杂度均为：\(T(n)=2T(\dfrac{n}{2})+O(n\log n)=O(n\log ^2 n)\)

代码细节多的要命。譬如 \(G(x)\) 的存储问题啊，还有求出来的 \(F_0(x)\) 可以直接放到原本 \(F(x)\) 的位置上以减小空间之类的啊，很不好写。

提交记录

void modp(ll *f,ll *g,ll n,ll m){
	while(!f[n-1]) n--;
	while(!g[m-1]) m--;
	if(n<m){
		for(int i=0;i<m;i++) g[i] = 0;
		for(int i=0;i<n;i++) g[i] = f[i] , f[i] = 0;
		return;
	}static ll q[N],r[N];
	ll lim=n-m+1,tmp=get2p(2*n);
	reverse(f,f+n),reverse(g,g+m);
	for(int i=0;i<lim;i++) q[i] = f[i],r[i] = g[i];
	invp(r,lim),mul(q,r,lim,lim);reverse(q,q+lim);//算 q
	reverse(g,g+m),reverse(f,f+n);mul(g,q,n,n);
	for(int i=0;i<n;i++) g[i] = (f[i] - g[i] + mod) %mod;
//	for(int i=0;i<lim;i++) f[i] = q[i];
//	for(int i=lim;i<=tmp;i++) f[i] = 0;
	for(int i=m;i<=tmp;i++) g[i] = 0;
	for(int i=0;i<=tmp;i++) q[i] = r[i] = 0;
}//f=gq+r;f->q,g->r;
ll n,m,f[N],_x[N];
#define L (2*l-2)
#define R (2*r-1)
#define mL (2*mid-1)
#define mR (2*mid)
ll g[2*N][20];
//下标 l~r 对应空间的 (2l-2)~(2r-1) 
void getg(ll l,ll r,ll dep){//l~r 为下标，L~R 为空间位置 
	if(r==l){
		g[L][dep] = (mod-_x[l]) %mod;
		g[R][dep] = 1;
		return;
	}ll mid=(l+r)>>1;
	getg(l,mid,dep+1),getg(mid+1,r,dep+1);
	static ll tmpf[N],tmpg[N];
	ll len=mL-L+1;//左区间长度>=右区间 
	for(int i=L;i<=mL;i++) tmpf[i-L] = g[i][dep+1];
	for(int i=mR;i<=R;i++) tmpg[i-mR] = g[i][dep+1];
	mul(tmpf,tmpg,len,R-L);
	for(int i=L;i<=R;i++) g[i][dep] = tmpf[i-L];
	ll tmp=get2p(2*len);
	for(int i=0;i<=tmp;i++) tmpf[i] = tmpg[i] = 0;
	return;
}
ll ans[N],tmpf[N];
void solve(ll l,ll r,ll dep,ll *f,ll n){//传进去的是项数！！ 
	if(r==l){
		ans[l] = f[0];
		return;
	}ll mid=(l+r)>>1,len=n>>1;
	static ll tmpg[N]={0};
	for(int i=L;i<=mL;i++) tmpg[i-L] = g[i][dep+1];//取出来！ 
	modp(f,tmpg,n,n);//传进去项数！！！
	for(int i=0;i<len;i++) tmpf[i] = tmpg[i];
	
	for(int i=mR;i<=R;i++) tmpg[i-mR] = g[i][dep+1];
	modp(f,tmpg,n,n);
	for(int i=0;i<len;i++) tmpf[i+len] = tmpg[i];
	
	for(int i=0;i<n;i++) f[i] = tmpf[i];
	solve(l,mid,dep+1,f,len);
	solve(mid+1,r,dep+1,f+len,len);
	return;
}
int main(){
	scanf("%lld%lld",&n,&m);n++;
	for(int i=0;i<n;i++) scanf("%lld",&f[i]);
	for(int i=1;i<=m;i++) scanf("%lld",&_x[i]);
	ll tmp=get2p(max(n,m));
	getg(1,tmp,0);
	solve(1,tmp,0,f,tmp);
	for(int i=1;i<=m;i++){
		printf("%lld\n",ans[i]);
	}
	return 0;
}

多项式复合

P5373 【模板】多项式复合函数

给定多项式 \(F(x),G(x)\)，求 \(G\Big(F(x)\Big)\)。

分块。设 \(L=\sqrt n\)。

题解讲的很详细。我有空再补。

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

提交记录

多项式复合逆

给定 \(F(x),G(x)\)，若 \(G\Big(F(x)\Big)=x\)，则 \(G(x),F(x)\) 互为复合逆。这意味着：

\[G\Big(F(x)\Big)=x\Leftrightarrow F\Big(G(x)\Big)=x \]

得看你索求的只是某一项的系数还是全部系数了。

P5809 【模板】多项式复合逆 这 b 题要求全部系数。不会。感觉没啥用。

而如果已知 \(F(x)\)，只要求 \(G(x)\) 的某一项系数的值，有拉格朗日反演：

\[[x^n]G(x)=\frac{1}{n}[x^{-1}]F(x)^{-n} \]

且 \(F(x),G(x)\) 应常数项为零，一次项系数非零。

很抽象。需要扩域扩到分式域。

当然有个更方便的形式：

\[[x^{-1}]F(x)^{-n}=[x^{dn-1}]\left(\frac{F(x)}{x^d}\right)^{-n} \]

其中 \(d\) 为 \(F(x)\) 系数里前缀零的个数，换句话说，\(F(x)\) 除以 \(x^d\) 后就能求逆了。总而言之：

\[[x^n]G(x)=\frac{1}{n}[x^{dn-1}]\left(\frac{F(x)}{x^d}\right)^{-n} \]

有的题需要求复合逆某一项的值，比如说：

• BZOJ3684. 大朋友和多叉树 推荐这篇博客：Autoint's Blog ，讲了解法以及拉格朗日反演的东西，很清晰。

Chirp Z 变换

P6800 【模板】Chirp Z-Transform

又称 Bluestein 算法。

给定多项式 \(A(x)\) 和 \(c\)，求 \(A(1,c,c^2,c^3\dots c^{m-1})\)。

对于 \(A(c^j)=\sum\limits_{i=0}^{n-1}a_ic^{ij}\)，我们有：\(ij=\dbinom{i+j}{2}-\dbinom{i}{2}-\dbinom{j}{2}\)，于是：

\[\begin{aligned} A(c^j)&=\sum\limits_{i=0}^{n-1}a_ic^{\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}}\\ &=c^{-\binom{j}{2}}\sum\limits_{i=0}^{n-1}a_ic^{\binom{i+j}{2}}c^{-\binom{i}{2}} \end{aligned}\]

设 \(f_i=a_ic^{-\binom{i}{2}}\)，\(g_{n-1+m-1-i}=c^{\binom{i}{2}}\)，我们有：

\[\begin{aligned} A(c^j)&=c^{-\binom{j}{2}}\sum\limits_{i=0}^{n-1}f_ig_{n-1-m-1-i}\\ &=c^{-\binom{j}{2}}\sum\limits_{i=0}^{n-1+m-1-j}f_ig_{n-1-m-1-j-i} \end{aligned}\]

就可以卷积了。

提交记录

\[\begin{aligned}\\\\\\\\\\\\\\\\\\\\\\\\\end{aligned} \]