FFT/NTT复习笔记&多项式&生成函数学习笔记Ⅲ

第三波，走起~~

FFT/NTT复习笔记&多项式&生成函数学习笔记Ⅰ

FFT/NTT复习笔记&多项式&生成函数学习笔记Ⅱ

单位根反演

今天打多校时 1002 被卡科技了……赛场上看出来是个单位根反演但不会，所以只好现学这东西了（

首先你得知道单位根是什么东西，对于 \(n\) 次方程 \(x^n-1=0(x\in\mathbb{C})\)，在复数域上有 \(n\) 个根，其对应到复平面上就是单位圆的 \(n\) 等分点，我们将这些单位根从 \(x\) 轴正半轴开始顺时针依次标号 \(0,1,2,\cdots,n-1\)，记作 \(\omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}\)，显然单位根满足以下性质：

• \(\omega_n^x=\cos\dfrac{2\pi x}{n}+\sin\dfrac{2\pi x}{n}\text{i}\)
• \((\omega_n^x)^n=1\)
• \(\omega_n^x=(\omega_n^1)^x\)

相信学过 FFT 的同学都不难理解。

单位根反演说的大概是这样一件事，对于任意整数 \(n\) 和 \(x\)，有 \([n\mid x]=\dfrac{1}{n}\sum\limits_{i=0}^{n-1}(\omega_n^{x})^i\)。也就是说如果 \(x\) 是 \(n\) 的倍数那么 \(\dfrac{1}{n}\sum\limits_{i=0}^{n-1}(\omega_n^{x})^i=1\)，否则 \(\dfrac{1}{n}\sum\limits_{i=0}^{n-1}(\omega_n^{x})^i=0\)。

证明：

• 如果 \(n\mid x\)，那么显然 \(\omega_n^x=1\)，\(\sum\limits_{i=0}^{n-1}(\omega_n^x)^i=n\)
• 如果 \(n\nmid x\)，那么根据等比数列求和公式 \(\sum\limits_{i=0}^{n-1}(\omega_n^{x})^i=\dfrac{(\omega_n^x)^n-1}{\omega_n^x-1}\)，显然分子等于 \(0\)，而由于 \(n\nmid x\)，\(\omega_n^x-1\ne 0\)，因此分式值为 \(0\)。

那么什么样的题能用单位根反演呢？首先如果题目给你一个带组合数的式子，这个组合数中带有取模符号 \(\bmod\)，而往往这个组合数上底数非常大（高达 \(10^9\) 甚至 \(10^{18}\)），而下底数中含枚举变量，那么此时这道题大概率就是单位根反演，运用二项式定理将组合数转化为幂的形式求解。注意，这里的 \(\bmod y\) 后面的 \(y\) 应当在输入中确定了，如果 \(y\) 作为枚举变量则无法二项式定理，还有一般单位根反演的题都涉及到取模，如果 \(y\nmid MOD-1\) 那么也无法二项式反演，因为我们要通过 \(g^{(MOD-1)/y}\) 求出 \(\omega_y^1\bmod MOD\)（譬如前天晚上的 D1C 如果模数 \(=10^9+9\) 就可以单位根反演了？）

说了这么多，还是具体题目具体分析吧。

60. LOJ #6485. LJJ 学二项式定理

首先这题组合数的 \(n\) 非常大，直接算的话第一步就爆了，因此我们需要想着用二项式定理将这里的 \(n\) 放到指数上去，怎么办呢？注意到这里涉及到 \(\bmod 4\)，而刚好有 \(4\mid 998244352\)，因此可以考虑单位根反演，上来推一波式子：

\[\begin{aligned} ans&=\sum\limits_{i=0}^n\dbinom{n}{i}s^i\sum\limits_{j=0}^3a_j[j\equiv i\pmod 4]\\ &=\sum\limits_{j=0}^3a_j\sum\limits_{i=0}^n\dbinom{n}{i}s^i[4\mid(i-j)]\\ &=\sum\limits_{j=0}^3a_j\sum\limits_{i=0}^n\dbinom{n}{i}s^i\sum\limits_{k=0}^3\omega_4^{(i-j)k}·\dfrac{1}{4}\\ &=\dfrac{1}{4}\sum\limits_{j=0}^3a_j\sum\limits_{k=0}^3\dfrac{1}{\omega_4^{jk}}\sum\limits_{i=0}^n\dbinom{n}{i}s^i\omega_4^{ik}\\ &=\dfrac{1}{4}\sum\limits_{j=0}^3a_j\sum\limits_{k=0}^3\dfrac{1}{\omega_4^{jk}}(s\omega_4^k+1)^n \end{aligned} \]

然后直接计算就行了，众所周知 \(998244353\) 的一个原根 \(g=3\)，那么我们可以通过 \(g^{998244352/4}\) 求出 \(\omega_4^1\)

const int MOD=998244353;
const int INV4=(3ll*MOD+1)>>2;
int w[4];
int qpow(int x,int e){
	int ret=1;
	for(;e;e>>=1,x=1ll*x*x%MOD) if(e&1) ret=1ll*ret*x%MOD;
	return ret;
}
ll n;int s,a[4];
void solve(){
	int res=0;scanf("%lld%d%d%d%d%d",&n,&s,&a[0],&a[1],&a[2],&a[3]);n%=(MOD-1);
	for(int i=0;i<4;i++) for(int j=0;j<4;j++)
		res=(res+1ll*a[i]*w[(16-j*i)&3]%MOD*qpow((1ll*w[j]*s+1)%MOD,n))%MOD;
	printf("%d\n",1ll*res*INV4%MOD);
}
int main(){
	w[0]=1;w[1]=qpow(3,MOD-1>>2);w[2]=1ll*w[1]*w[1]%MOD;
	w[3]=1ll*w[1]*w[1]%MOD*w[1]%MOD;
	int qu;scanf("%d",&qu);while(qu--) solve();
	return 0;
}

61. P5591 小猪佩奇学数学

还是注意到此题 \(n\) 很大，而 \(k\) 是 \(2\) 的整数次幂，众所周知 \(998244353-1=119·2^{23}\) 刚好是 \(k\) 的倍数，因此题目已经疯狂暗示了此题的算法，然后就考虑一阵推式子：

\[\begin{aligned} ans&=\sum\limits_{i=0}^n\dbinom{n}{i}p^i·\dfrac{i-i\bmod k}{k}\\ &=\dfrac{1}{k}\sum\limits_{i=0}^n\dbinom{n}{i}p^i·i-\dfrac{1}{k}\sum\limits_{i=0}^n\dbinom{n}{i}p^i·i\bmod k\\ &=\dfrac{1}{k}\sum\limits_{i=0}^n\dbinom{n-1}{i-1}p^i·n-\dfrac{1}{k}\sum\limits_{i=0}^n\dbinom{n}{i}p^i·\sum\limits_{r=0}^{k-1}r[k\mid i-r]\\ &=\dfrac{np}{k}\sum\limits_{i=0}^n\dbinom{n-1}{i-1}p^{i-1}-\dfrac{1}{k^2}\sum\limits_{i=0}^n\dbinom{n}{i}p^i·\sum\limits_{r=0}^{k-1}r\sum\limits_{j=0}^{k-1}\omega_k^{j(i-r)}\\ &=\dfrac{np}{k}(p+1)^{n-1}-\dfrac{1}{k^2}\sum\limits_{i=0}^n\dbinom{n}{i}p^i·\sum\limits_{r=0}^{k-1}r\sum\limits_{j=0}^{k-1}\omega_k^{j(i-r)}\\ &=\dfrac{np}{k}(p+1)^{n-1}-\dfrac{1}{k^2}\sum\limits_{j=0}^{k-1}(\sum\limits_{r=0}^{k-1}r·\omega_k^{-jr})·(\sum\limits_{i=0}^n\dbinom{n}{i}p^i\omega_k^{ji})\\ &=\dfrac{np}{k}(p+1)^{n-1}-\dfrac{1}{k^2}\sum\limits_{j=0}^{k-1}(\sum\limits_{r=0}^{k-1}r·\omega_k^{-jr})·(p^i\omega_k^{i}+1)^n \end{aligned} \]

式子推到这里，聪明的读者们不难发现我们复杂度瓶颈在于如何计算 \(\sum\limits_{r=0}^{k-1}r·\omega_k^{-jr}\)，显然这东西可以被抽象为 \(\sum\limits_{i=0}^{n-1}i·a^i\)，因此接下来我们考虑如何高效求出这个式子的值（其实《具体数学》上有这方面的推导内容，大概在 P27），我们设 \(S=\sum\limits_{i=0}^{n}ia^i\)，那么 \(aS=\sum\limits_{i=0}^{n}ia^{i+1}=\sum\limits_{i=1}^{n+1}(i-1)a^i\)，两式相减可以得到 \((a-1)S=(n+1)a^{n+1}-\sum\limits_{i=1}^{n+1}a^i\)，然后一个等比数列求和公式带走即可，注意特判 \(a=1\)。

复杂度 \(k\log k\)。

62. UOJ #450. 【集训队作业2018】复读机

首先当 \(d=1\) 时答案显然是 \(m^n\)

当 \(k=2\) 时每个复读机的 EGF 显然是 \(\sum\limits_{2\mid i}\dfrac{x^i}{i!}=\cosh x=\dfrac{e^x+e^{-x}}{2}\)，因此答案即为 \(n^k\)，二项式定理一通暴算即可，复杂度 \(\mathcal O(k)\)。

当 \(k=3\) 时类似地有每个复读机的 EGF 为 \(\sum\limits_{3\mid i}\dfrac{x^i}{i!}=\sum\limits_{i}\dfrac{x^i}{i!}[3\mid i]=\dfrac{1}{3}\sum\limits_{j=0}^2\sum\limits_{i}\dfrac{x^i·\omega_3^{ij}}{i!}\)，化简一下可以得到 \(\dfrac{1}{3}(\sum\limits_{j=0}^2e^{\omega_3^jx})\)，然后三项式定理（大雾）枚举一下即可，复杂度 \(\mathcal O(k^2)\)

63. HDU 7013 String Mod

好了，既然你已经学会了这么多，就来试一下昨天的多校 1002 吧（

首先对于一对 \((i,j)\)，答案显然可以写成 \(\sum\limits_{p}\sum\limits_{q}\dbinom{L}{p,q,L-p-q}·(k-2)^{L-p-q}[p\equiv i\pmod{n}][q\equiv j\pmod{n}]\)，看到 \(\bmod\) 以及 \(n\mid P-1\)，可以考虑单位根反演，即

\[\begin{aligned} ans&=\sum\limits_{p}\sum\limits_{q}\dbinom{L}{p,q,L-p-q}·(k-2)^{L-p-q}[p\equiv i\pmod{n}][q\equiv j\pmod{n}]\\ &=\dfrac{1}{n^2}\sum\limits_{p}\sum\limits_{q}\dbinom{L}{p,q,L-p-q}·(k-2)^{L-p-q}\sum\limits_{x=0}^{n-1}\sum\limits_{y=0}^{n-1}\omega_n^{(p-i)x}\omega_n^{(q-j)y}\\ &=\dfrac{1}{n^2}\sum\limits_{x=0}^{n-1}\sum\limits_{y=0}^{n-1}\dfrac{1}{\omega^{ix}}\dfrac{1}{\omega^{jy}}\sum\limits_{p}\sum\limits_{q}\dbinom{L}{p,q,L-p-q}·(k-2)^{L-p-q}(\omega_n^{x})^p(\omega_n^{y})^q\\ &=\dfrac{1}{n^2}\sum\limits_{x=0}^{n-1}\sum\limits_{y=0}^{n-1}\dfrac{1}{\omega^{ix}}\dfrac{1}{\omega^{jy}}(\omega_n^x+\omega_n^y+k-2)^L \end{aligned} \]

直接算是 \(n^2\) 的，再加上前面 \(i,j\) 的复杂度，总复杂度 \(n^4\)，不过发现这个式子中每一项要么与 \(i,x\) 有关，要么与 \(j,y\) 有关，要么与 \(x,y\) 有关，并没有哪一项与 \(i,j\) 直接相关，因此考虑矩阵乘法，我们记矩阵 \(A_{i,j}=\dfrac{1}{\omega_n^{ij}}\)，矩阵 \(B_{i,j}=(\omega_n^i+\omega_n^j+k-2)^L\)，那么答案矩阵即为 \(A\times B\times A\) 乘 \(\dfrac{1}{n^2}\)。

总复杂度 \(\mathcal O(n^3)\)​

FFT 三次变两次的小技巧

一个非常 trivial 的小 trick。8 个月前学 FFT 时没学懂现在来补了。

首先注意到 \(\text{DFT}\) 是一个线性变换，因此对于两个序列 \(a,b\)，记 \(va\) 为 \(a\) 每一项都乘以 \(v\) 后的结果，那么有 \(\text{DFT}(a+b\sqrt{-1})=\text{DFT}(a)+\text{DFT}(b)\sqrt{-1}\)，也就是说，假设 \(P=\text{DFT}(a),Q=\text{DFT}(b)\)，那么 \(\text{DFT}(a+b\sqrt{-1})=P+\sqrt{-1}Q\)。

这样还是无法确定出 \(P,Q\)​​，不过注意到有个东西叫共轭复数，具体来说，对于 \(z=a+b\sqrt{-1}\)​​，其共轭复数 \(\overline{z}=a-b\sqrt{-1}\)​​，因此对于一个序列 \(a\)​​，我们考虑定义其共轭序列（瞎起名字 ing）\(\overline{a}\)​​ 满足 \(\overline{a}_i=\overline{a_i}\)​​，那么 \(\text{DFT}(\overline{a+b\sqrt{-1}})=\text{DFT}(a-b\sqrt{-1})=P-\sqrt{-1}Q\)​​，因此如果我们能高效求出 \(\text{DFT}(\overline{a+b\sqrt{-1}})\)​，即 \(\text{DFT}(a-b\sqrt{-1})\)​，那么对于每一项 \(i\)​，如果我们设 \(X=\text{DFT}(a+b\sqrt{-1})_i,Y=\text{DFT}(a-b\sqrt{-1})_i\)​，有 \(P_i+\sqrt{-1}Q_i=X,P_i-\sqrt{-1}Q_i=Y\)​​​，简单二元一次方程组即可求出 \(P_i,Q_i\)​​。具体来说，假设 \(x_1=P_i\) 的实部，\(y_1=P_i\) 的虚部，\(x_2,y_2\) 也同理，那么：

\[\begin{cases} x_1-y_2=\text{Re}(X)\\ y_1+x_2=\text{Im}(X)\\ x_1+y_2=\text{Re}(Y)\\ y_1-x_2=\text{Im}(Y) \end{cases} \]

解得

\[\begin{cases} x_1=\dfrac{\text{Re}(X)+\text{Re}(Y)}{2}\\ y_1=\dfrac{\text{Im}(X)+\text{Im}(Y)}{2}\\ x_2=\dfrac{\text{Im}(X)-\text{Im}(Y)}{2}\\ y_2=\dfrac{\text{Re}(Y)-\text{Re}(Y)}{2}\\ \end{cases} \]

那么怎么求出 \(\text{DFT}(\overline{a+b\sqrt{-1}})\)​ 呢？注意到 DFT 的过程实际上是将每个 \(a_i\)​ 变为 \(f(\omega_n^i)\)​​​​，因此 \(\text{DFT}(\overline{a})_i=\sum\limits_{j=0}^n\overline{a}_j(\omega_n^i)^j\)​，又有：

• \(\overline{ab}=\overline{a}·\overline{b}\)
• \(\overline{\omega_{n}^i}=\omega_n^{n-i}\)

因此 \(\text{DFT}(\overline{a})_i=\sum\limits_{j=0}^n\overline{a}_j\overline{(\omega_n^{n-i})^j}=\overline{\sum\limits_{j=0}^na_j(\omega_n^{n-i})^j}\)，不难发现对于 \(i=0\) 而言，由于 \(\omega_{n}^n=1\)，后面那东西就是 \(\overline{\text{DFT}(a)_0}\)，而对于 \(i\ne 0\)，后面那东西就是 \(\overline{\text{DFT}(a)_{n-i}}\)，因此我们求出 \(\text{DFT}(a+b\sqrt{-1})\) 后，把这个序列第 \(1\) 项至第 \(n-1\) 项 reverse 一下并把每一项的虚部变为其相反数即可得到 \(\text{DFT}(\overline{a+b\sqrt{-1}})\)​。这样即可一次 DFT+一次 IDFT 实现 FFT（

const int MAXP=1<<21;
const double Pi=acos(-1);
int n,m,LEN=1;
struct comp{
	double x,y;
	comp(double _x=0,double _y=0):x(_x),y(_y){}
	comp operator +(const comp &rhs){return comp(x+rhs.x,y+rhs.y);}
	comp operator -(const comp &rhs){return comp(x-rhs.x,y-rhs.y);}
	comp operator *(const comp &rhs){return comp(x*rhs.x-y*rhs.y,x*rhs.y+y*rhs.x);}
	comp operator /(const double &rhs){return comp(x/rhs,y/rhs);}
} a[MAXP+5],b[MAXP+5];
int rev[MAXP+5];
void FFT(comp *a,int len,int type){
	int lg=31-__builtin_clz(len);
	for(int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<lg-1);
	for(int i=0;i<len;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=2;i<=len;i<<=1){
		comp W=comp(cos(2*Pi/i),type*sin(2*Pi/i));
		for(int j=0;j<len;j+=i){
			comp w=comp(1,0);
			for(int k=0;k<(i>>1);k++,w=w*W){
				comp X=a[j+k],Y=w*a[(i>>1)+j+k];
				a[j+k]=X+Y;a[(i>>1)+j+k]=X-Y;
			}
		}
	} if(!~type) for(int i=0;i<len;i++) a[i].x=(int)(a[i].x/len+0.5);
}
int main(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%lf",&a[i].x);
	for(int i=0;i<=m;i++) scanf("%lf",&a[i].y);
	while(LEN<=n+m) LEN<<=1;FFT(a,LEN,1);
	for(int i=0;i<LEN;i++){
		comp x=a[i],y=(!i)?a[i]:a[LEN-i];y.y=-y.y;
		comp A=(x+y)/2.0,B=(x-y)/2.0;swap(B.x,B.y);B.y=-B.y;
		b[i]=A*B;
	} FFT(b,LEN,-1);
	for(int i=0;i<=n+m;i++) printf("%d%c",(int)b[i].x," \n"[i==n+m]);
	return 0;
}

各种多项式的 \(n^2\) 版本

qwq 大概是为了为下面子集卷积的高级操作做铺垫（？）

博客戳这儿

子集卷积的一些高级操作

首先在开始正题之前我们先得知道什么是集合幂级数，集合幂级数大概就是一个形如 \(f(x)=\sum\limits_{S\subseteq U}a_Sx^S\) 的幂级数 \(f\)，其中 \(U\) 是一个给定的有限集合。与形式幂级数不同的是，\(x^S\) 上方的指数在这种定义下全部是全集 \(U\) 的一个个子集 instead of 一个个非负整数（当然，在 OI 中，由于用二进制表示集合这种方法的存在，我们也可以将 \(x\) 上方的指数视作一个个在 \([0,2^{|U|})\) 中的二进制数，这样全集可视为 \(2^{|U|}-1\)）。与形式幂级数的加法卷积类似，我们也可以定义 or 卷积、and 卷积、xor 卷积等，分别表示两个集合幂级数 \(F,G\) 相乘，\([x^S]F(x)\times[x^T]G(x)\) 的系数贡献到 \(x^{S\cup T},x^{S\cap T},x^{S\oplus T}\) 的情况，对于这些卷积你都可以在我的快速沃尔什变换&快速莫比乌斯变换小记中找到，当然由于这里我们将着重探讨子集卷积，所以对于其他的卷积咱也大可不必花心思赘述了（

子集卷积，大概也是 \([x^S]F(x)\times[x^T]G(x)\)​ 的系数贡献到 \(x^{S\cup T}\)​ 的情况，不过有一个前提就是 \(S\cap T=\varnothing\)​ 时才能产生贡献，否则贡献为 \(0\)​。那么怎么快速求出两个集合幂级数 \(F,G\)​ 的子集卷积呢？注意到关于集合的并，有一个性质 \(|S|+|T|\ge|S\cup T|\)​，因此我们考虑将集合幂级数写成一个二元函数的形式（形式而已，并不是真正在集合幂级数后面加一维），即将每个 \(a_Sx^S\)​ 扩展成 \(a_Sx^Sy^{|S|}\)​​ 的形式，这样每个集合幂级数可以写作 \(\sum\limits_{n\ge 0}\sum\limits_{m\ge 0}a_{n,m}x^ny^m\)​​，对于两个集合幂级数 \(F,G\)，假设它们的系数分别为 \(f_{n,m}\) 和 \(g_{n,m}\)，那么我们重定义它们的乘法为 \(F\times G=\sum\limits_{n\ge 0}\sum\limits_{m\ge 0}(\sum\limits_{p|q=n}\sum\limits_{r+s=m}f_{p,r}g_{q,s})x^ny^m\)，也就是 \(f_{p,r}g_{q,s}\) 贡献到 \(x^{p|q}y^{r+s}\) 的位置，这样我们取它们乘积 \(x^Sy^{|S|}\) 前的系数作为它们子集卷积 \(x^S\) 前的系数即可，不难发现这样一来，两个集合 \(S,T\) 如果满足 \(S\cap T\ne\varnothing\)，那么其就一定不满足 \(|S\cap T|=|S|+|T|\)，也就一定不会贡献到 \(x^{S\cup T}y^{|S\cup T|}\) 处，也就证明了上述做法的正确性。

那么这东西该怎么实现呢，我们考虑将集合幂级数 \(F(x)\) 扩充之后的样子写成 \(\sum\limits_{S}x^SG_S(y)\)​ 的形式，其中 \(G_S(y)\) 是一个关于 \(y\)​​ 的形式幂级数，这样不难发现我们的乘法运算可以视作：将每个集合幂级数前的系数视作一个形式幂级数，系数的乘法视作多项式的乘法，然后对它们做一遍 or 卷积。因此考虑按照套路 FWTor，然后对每个 \(x^S\)，将两个集合幂级数 \(x^S\) 前的系数进行多项式乘法，然后 IFWTor 回去即可。复杂度 \(2^nn^2\)。

注意到这里涉及形式幂级数的操作，因此我们也可以按照形式幂级数定义它们的逆、\(\ln\)、\(\exp\)。具体来说，对于集合幂级数 \(F\)，我们定义：

• 逆：定义 \(F\) 的逆为满足 \(F\times G=\epsilon\) 的集合幂级数 \(G\)，其中 \(\times\) 在这部分默认为子集卷积，\(\epsilon\) 为满足 \([x^S]\epsilon=[S=\varnothing]\) 的集合幂级数。
• \(\exp\)：定义 \(\exp F\) 为满足 \(G=\sum\limits_{n\ge 0}\dfrac{F^n}{n!}\) 的集合幂级数 \(G\)
• \(\ln\)：定义 \(\ln F\) 为满足 \(\exp G=F\) 的集合幂级数 \(G\)
• \(\exp_{\le k}\)：定义 \(\exp_{\le k}\)​ 为满足 \(G=\sum\limits_{n=0}^k\dfrac{F^n}{n!}\) 的集合幂级数

由于我们子集卷积可以看作是将集合幂级数都一遍 FWTor，然后每一项按照形式幂级数的套路卷起来，再 IFWTor，因此上述操作也可视作，将集合幂级数 FWTor 一遍，然后每一个 \(x^S\) 前的系数分别 \(\text{inv}/\ln/\exp/\exp_{\le k}\)。注意，由于集合幂级数的题 \(n\) 一般都很小，因此在这里就大可不必 NTT 实现多项式 \(\text{inv}/\ln/\exp/\exp_{\le k}\) 了，可以平方地递推，具体递推方法可见上面那篇 blog。

64. P6570 [NOI Online #3 提高组] 优秀子序列

题解

65. P6846 [CEOI2019] Amusement Park

注意到对于一种翻转边集的状态，如果将所有边的状态反转，那么得到的图仍是一个 DAG，也就是说我们可以将所有翻转某个边集，使得得到的图是一个 DAG 的方案两两配对，并使每对中翻转的边集数量之和恰好等于 \(m\)。因此我们只用算出有多少种翻转边集之后是 DAG 的方案然后乘上 \(\dfrac{m}{2}\) 即可。

考虑怎么算翻转之后是 DAG 的方案，我们考虑设 \(f_S\) 表示集合 \(S\)​ 中形成 DAG 的方案数，那么考虑钦定一些点入度为 \(0\)，记这个集合为 \(T\)，那么显然 \(T\) 必须是一个独立集，否则总会存在 \(T\)​ 中的点入度非零，考虑枚举这个集合 \(T\)，那么显然 \(S-T\) 也应是一个 DAG，但这样会算重，因此考虑容斥，即

\[f_S=\sum\limits_{\varnothing\ne T\subseteq S}[T\text{ is an independent set}]·(-1)^{|T|+1}f_{S-T} \]

由于 Latex 上用中文太鸡肋了就换用了英文

注意到这东西可以写成子集卷积的形式，具体来说，即 \(F\) 为 \(f\) 的“生成集合幂级数”，\(G(x)\) 为满足 \([x^T]G(x)=[T\text{ is an independent set}]·(-1)^{|T|+1}\)，那么 \(FG=F-1\)，即 \(F=\dfrac{1}{1-G}\)，子集 inv 一波即可。

66. LOJ #154 集合划分计数

搞不懂模板题卡常是啥心态

记 \(F\)​ 为满足 \([x^S]F(x)=[x\in\mathcal F]\)，那么我们对 \(F\) 作 \(\exp_{\le k}\) 即可。

只不过此题过于卡常，所以交了 114514191981019260817998244353 发

卡常技巧：

1. 求 \(\text{exp}_k\) 时候，如果该多项式最低非零项次数 \(>\dfrac{n}{k}\) 就不用快速幂了，直接将快速幂得到的数组每一项设为 \(0\) 即可。
2. 在暴力多项式乘法时，可以找到两个多项式最低和最高的非零项 \(la,ra\) 和 \(lb,rb\)，然后只枚举次数在 \([la,ra]\) 和 \([lb,rb]\) 的范围内的数即可，这样效率可以高不少。

不过还是几乎喜提最劣解

人傻常数大，需要狠命卡

67. LOJ #6673 EntropyIncreaser 与山林

首先先不考虑连通这个条件，考虑怎样求有多少张子图 \(G’=(V’,E’)\)，满足 \(V’=S\)，且 \(G’\) 中每个连通块都有欧拉回路，记这个方案数为 \(f_S\)。显然，对于一张图而言，每个连通块都有欧拉回路的充要条件是每个点的度都是偶数，因此我们考虑对每条边新建一个布尔型变量表示其有没有选，这样我们可以列出 \(|S|\) 个方程组，高斯消元一下，那么 \(f_S\) 就是 \(2\) 的自由元个数次方。

接下来考虑带上连通这个条件后怎么处理，我们记 \(F(x)\) 为满足 \([x^S]F(x)=f_S\) 的集合幂级数，再设 \(G(x)\) 表示 \([x^S]G(x)=\) 以 \(S\) 为点集的、存在欧拉回路（不是每个连通块都有欧拉回路）的子图个数，那么根据 P4841 的套路有 \(\exp G=F\)。现在已知 \(F\)，子集 \(\ln\) 一波即可求出 \(G\)。

68. UOJ #94 【集训队互测2015】胡策的统计

考虑记 \(f_S\)​ 为以 \(S\)​ 为点集的子图个数，那么显然有 \(f_S=2^{cnt_S}\)​，其中 \(cnt_S\)​ 为 \(S\)​ 内部边的数量，再设 \(g_S\)​ 为以 \(S\)​ 为点集的连通子图个数，那么设 \(F,G\)​ 分别为 \(f,g\)​ 对应的集合幂级数，显然就有 \(F=\exp G\)​，\(\ln\)​ 一下可求出 \(G\)​。再设 \(h_S\)​ 为以集合 \(S\)​ 为点集的所有子图的连通值之和，也对应地设其生成集合幂级数为 \(H\)​，那么 \(H=\sum\limits_{n\ge 0}G^n\)，也就是 \(\exp\) 里面去掉分母上的 \(n!\)，也就相当于对每个 \(n\) 个连通块的方案乘了个 \(n!\)，注意到这东西显然等于 \(\dfrac{1}{1-G}\)，因此再一波子集 \(\ln\) 求出 \(H\)，最终答案就是 \(H_U\)。

69. NFLSOJ #1060. 【2021 六校联合训练 NOI #40】白玉楼今天的饭

题解

多项式多点求值

好，关于多点求值，它终于出现在我的 blog 中了（

所谓多点求值，就是给定多项式 \(f(x)\) 以及 \(a_1,a_2,\cdots,a_n\)，要快速求出其在 \(a_i\) 处的点值，即求出 \(f(a_1),f(a_2),\cdots,f(a_n)\)。直接做是平方的，考虑如何优化它。

首先有引理 \(f(v)=[x^0](f(x)\bmod (x-v))\)​​​，其中 mod 为多项式带余除法。证明就考虑归纳，对于 \(0\)​​​ 次多项式命题显然成立，而对于 \(n\)​​​ 次多项式 \(F(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n\)​​​，我们首先将 \(n\)​​​ 次项次数变为 \(0\)​​​，即减去 \(a_nx^{n}-a_n·v·x^{n-1}\)​​​，此时多项式变为 \(G(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-2}x^{n-2}+(a_{n-1}+a_nv)·x^{n-1}\)​​​，根据归纳假设有 \(F(x)\bmod (x-v)=G(x)\bmod(x-v)=\sum\limits_{i=0}^{n-2}a_i·v^i+(a_{n-1}+a_nv)·v^{n-1}=\sum\limits_{i=0}^{n}a_i·v^i\)​，如此归纳下去即可。

知道了这个结论之后考虑怎样求解答案，我们先对点值序列做一遍类似于线段树的东西，即考虑分治，递归到 \([l,r]\) 时，我们预处理出 \(\prod\limits_{i=l}^r(x-a_i)\)，这个可以分治+FFT 处理出来，然后我们将原多项式放在根结点处然后往下递归，每到一个结点，就令当前多项模上该节点上的 \(\prod\limits_{i=l}^r(x-a_i)\)，递归到叶子时多项式肯定只剩 \(1\) 项了，此时这一项的值就是该点的点值。

算下复杂度，在第 \(i\) 层中我们要进行 \(2^i\) 次 fmod 操作，而第 \(i\) 层的多项式长度为 \(\dfrac{n}{2^i}\)，进行 fmod 的时间复杂度为 \(\dfrac{n}{2^i}·\log n\)，因此总时间复杂度 \(n\log^2n\)，只不过由于常数巨大跑起来不比 3log 快多少罢了（

多项式 fmod 除了模板没有实际用处！——tzc 早期迷惑发言（

赋上常数巨大的代码：

const int MAXN = 64000;
const int MAXP = 1 << 18;
const int pr = 3;
const int ipr = 332748118;
const int MOD = 998244353;
int qpow(int x, ll e) {
	int ret = 1;
	for (; e; e >>= 1, x = 1ll * x * x % MOD)
		if (e & 1) ret = 1ll * ret * x % MOD;
	return ret;
}
int n, m, v[MAXN + 5], rev[MAXP + 5], res[MAXN + 5];
void FFT(vector<int> &a, int len, int type) {
	int lg = 31 - __builtin_clz(len);
	for (int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg - 1);
	for (int i = 0; i < len; i++) if (rev[i] < i) swap(a[i], a[rev[i]]);
	static int W[MAXP + 5];
	for (int i = 2; i <= len; i <<= 1) {
		int w = qpow((type < 0) ? ipr : pr, (MOD - 1) / i);
		for (int j = (W[0] = 1); j < (i >> 1); j++) W[j] = 1ll * W[j - 1] * w % MOD;
		for (int j = 0; j < len; j += i) {
			for (int k = 0; k < (i >> 1); k++) {
				int X = a[j + k], Y = 1ll * W[k] * a[(i >> 1) + j + k] % MOD;
				a[j + k] = (X + Y) % MOD; a[(i >> 1) + j + k] = (X - Y + MOD) % MOD;
			}
		}
	}
	if (!~type) {
		int ivn = qpow(len, MOD - 2);
		for (int i = 0; i < len; i++) a[i] = 1ll * a[i] * ivn % MOD;
	}
}
vector<int> conv(vector<int> a, vector<int> b) {
	int LEN = 1, lim = a.size() + b.size() - 1;
	while (LEN < a.size() + b.size()) LEN <<= 1;
	a.resize(LEN, 0); b.resize(LEN, 0); FFT(a, LEN, 1); FFT(b, LEN, 1);
	for (int i = 0; i < LEN; i++) a[i] = 1ll * a[i] * b[i] % MOD;
	FFT(a, LEN, -1); while (a.size() > lim) a.ppb(); return a;
}
vector<int> getinv(vector<int> a, int len) {
	vector<int> b(len, 0); b[0] = qpow(a[0], MOD - 2); assert(a[0]);
	for (int i = 2; i <= len; i <<= 1) {
		vector<int> c(a.begin(), a.begin() + i);
		vector<int> d(b.begin(), b.begin() + (i >> 1));
		d = conv(conv(d, d), c);
		for (int j = 0; j < i; j++) b[j] = (2ll * b[j] - d[j] + MOD) % MOD;
	}
	return b;
}
vector<int> getfmod(vector<int> a, vector<int> b) {
	if (a.size() < b.size()) return a;
	if (a.size() == b.size()) {
		int c = 1ll * a.back() * qpow(b.back(), MOD - 2) % MOD;
		for (int i = 0; i < a.size(); i++) a[i] = (a[i] - 1ll * c * b[i] % MOD + MOD) % MOD;
		a.ppb(); return a;
	}
	vector<int> aR = a, bR = b;
	reverse(aR.begin(), aR.end()); reverse(bR.begin(), bR.end());
	int LEN = 1; while (LEN < bR.size()) LEN <<= 1; bR.resize(LEN, 0);
	vector<int> ibR = getinv(bR, LEN), mul = conv(ibR, aR);
	vector<int> q(a.size() - b.size() + 1);
	for (int i = 0; i <= a.size() - b.size(); i++) q[i] = mul[i];
	reverse(q.begin(), q.end()); mul = conv(q, b);
	vector<int> res(b.size() - 1);
	for (int i = 0; i + 1 < b.size(); i++) res[i] = (a[i] - mul[i] + MOD) % MOD;
	return res;
}
vector<int> vec[MAXN * 4 + 5];
void init(int k, int l, int r) {
	if (l == r) {vec[k].clear(); vec[k].pb(MOD - v[l]); vec[k].pb(1); return;}
	int mid = l + r >> 1; init(k << 1, l, mid); init(k << 1 | 1, mid + 1, r);
	vec[k] = conv(vec[k << 1], vec[k << 1 | 1]);
}
void solve(vector<int> cur, int k, int l, int r) {
	if (cur.size() >= vec[k].size()) cur = getfmod(cur, vec[k]);
	if (l == r) return res[l] = cur[0], void();
	int mid = l + r >> 1;
	solve(cur, k << 1, l, mid); solve(cur, k << 1 | 1, mid + 1, r);
}
int main() {
	read(n); read(m); vector<int> A;
	for (int i = 0, x; i <= n; i++) read(x), A.pb(x);
	for (int i = 1; i <= m; i++) read(v[i]);
	init(1, 1, m); solve(A, 1, 1, m);
	for (int i = 1; i <= m; i++) print(res[i], '\n');
	print_final();
	return 0;
}