第三波，走起~~
单位根反演
今天打多校时 1002 被卡科技了……赛场上看出来是个单位根反演但不会，所以只好现学这东西了（
首先你得知道单位根是什么东西，对于 (n) 次方程 (x^n-1=0(xinmathbb{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=cosdfrac{2pi x}{n}+sindfrac{2pi x}{n} ext{i})
- ((omega_n^x)^n=1)
- (omega_n^x=(omega_n^1)^x)
相信学过 FFT 的同学都不难理解。
单位根反演说的大概是这样一件事，对于任意整数 (n) 和 (x)，有 ([nmid x]=dfrac{1}{n}sumlimits_{i=0}^{n-1}(omega_n^{x})^i)。也就是说如果 (x) 是 (n) 的倍数那么 (dfrac{1}{n}sumlimits_{i=0}^{n-1}(omega_n^{x})^i=1)，否则 (dfrac{1}{n}sumlimits_{i=0}^{n-1}(omega_n^{x})^i=0)。
证明：
- 如果 (nmid x)，那么显然 (omega_n^x=1)，(sumlimits_{i=0}^{n-1}(omega_n^x)^i=n)
- 如果 (n mid x)，那么根据等比数列求和公式 (sumlimits_{i=0}^{n-1}(omega_n^{x})^i=dfrac{(omega_n^x)^n-1}{omega_n^x-1})，显然分子等于 (0)，而由于 (n mid x)，(omega_n^x-1 e 0)，因此分式值为 (0)。
那么什么样的题能用单位根反演呢？首先如果题目给你一个带组合数的式子，这个组合数中带有取模符号 (mod)，而往往这个组合数上底数非常大（高达 (10^9) 甚至 (10^{18})），而下底数中含枚举变量，那么此时这道题大概率就是单位根反演，运用二项式定理将组合数转化为幂的形式求解。注意，这里的 (mod y) 后面的 (y) 应当在输入中确定了，如果 (y) 作为枚举变量则无法二项式定理，还有一般单位根反演的题都涉及到取模，如果 (y mid MOD-1) 那么也无法二项式反演，因为我们要通过 (g^{(MOD-1)/y}) 求出 (omega_y^1mod MOD)（譬如前天晚上的 D1C 如果模数 (=10^9+9) 就可以单位根反演了？）
说了这么多，还是具体题目具体分析吧。
60. LOJ #6485. LJJ 学二项式定理
首先这题组合数的 (n) 非常大，直接算的话第一步就爆了，因此我们需要想着用二项式定理将这里的 (n) 放到指数上去，怎么办呢？注意到这里涉及到 (mod 4)，而刚好有 (4mid 998244352)，因此可以考虑单位根反演，上来推一波式子：
然后直接计算就行了，众所周知 (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
",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) 的倍数，因此题目已经疯狂暗示了此题的算法，然后就考虑一阵推式子：
式子推到这里，聪明的读者们不难发现我们复杂度瓶颈在于如何计算 (sumlimits_{r=0}^{k-1}r·omega_k^{-jr})，显然这东西可以被抽象为 (sumlimits_{i=0}^{n-1}i·a^i)，因此接下来我们考虑如何高效求出这个式子的值（其实《具体数学》上有这方面的推导内容，大概在 P27），我们设 (S=sumlimits_{i=0}^{n}ia^i)，那么 (aS=sumlimits_{i=0}^{n}ia^{i+1}=sumlimits_{i=1}^{n+1}(i-1)a^i)，两式相减可以得到 ((a-1)S=(n+1)a^{n+1}-sumlimits_{i=1}^{n+1}a^i)，然后一个等比数列求和公式带走即可，注意特判 (a=1)。
复杂度 (klog k)。
62. UOJ #450. 【集训队作业2018】复读机
首先当 (d=1) 时答案显然是 (m^n)
当 (k=2) 时每个复读机的 EGF 显然是 (sumlimits_{2mid i}dfrac{x^i}{i!}=cosh x=dfrac{e^x+e^{-x}}{2})，因此答案即为 (n^k)，二项式定理一通暴算即可，复杂度 (mathcal O(k))。
当 (k=3) 时类似地有每个复读机的 EGF 为 (sumlimits_{3mid i}dfrac{x^i}{i!}=sumlimits_{i}dfrac{x^i}{i!}[3mid i]=dfrac{1}{3}sumlimits_{j=0}^2sumlimits_{i}dfrac{x^i·omega_3^{ij}}{i!})，化简一下可以得到 (dfrac{1}{3}(sumlimits_{j=0}^2e^{omega_3^jx}))，然后三项式定理（大雾）枚举一下即可，复杂度 (mathcal O(k^2))
63. HDU 7013 String Mod
好了，既然你已经学会了这么多，就来试一下昨天的多校 1002 吧（
首先对于一对 ((i,j))，答案显然可以写成 (sumlimits_{p}sumlimits_{q}dbinom{L}{p,q,L-p-q}·(k-2)^{L-p-q}[pequiv ipmod{n}][qequiv jpmod{n}])，看到 (mod) 以及 (nmid P-1)，可以考虑单位根反演，即
直接算是 (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 imes B imes A) 乘 (dfrac{1}{n^2})。
总复杂度 (mathcal O(n^3))
FFT 三次变两次的小技巧
一个非常 trivial 的小 trick。8 个月前学 FFT 时没学懂现在来补了。
首先注意到 ( ext{DFT}) 是一个线性变换，因此对于两个序列 (a,b)，记 (va) 为 (a) 每一项都乘以 (v) 后的结果，那么有 ( ext{DFT}(a+bsqrt{-1})= ext{DFT}(a)+ ext{DFT}(b)sqrt{-1})，也就是说，假设 (P= ext{DFT}(a),Q= ext{DFT}(b))，那么 ( ext{DFT}(a+bsqrt{-1})=P+sqrt{-1}Q)。
这样还是无法确定出 (P,Q)，不过注意到有个东西叫共轭复数，具体来说，对于 (z=a+bsqrt{-1})，其共轭复数 (overline{z}=a-bsqrt{-1})，因此对于一个序列 (a)，我们考虑定义其共轭序列（瞎起名字 ing）(overline{a}) 满足 (overline{a}_i=overline{a_i})，那么 ( ext{DFT}(overline{a+bsqrt{-1}})= ext{DFT}(a-bsqrt{-1})=P-sqrt{-1}Q)，因此如果我们能高效求出 ( ext{DFT}(overline{a+bsqrt{-1}}))，即 ( ext{DFT}(a-bsqrt{-1}))，那么对于每一项 (i)，如果我们设 (X= ext{DFT}(a+bsqrt{-1})_i,Y= ext{DFT}(a-bsqrt{-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) 也同理，那么：
解得
那么怎么求出 ( ext{DFT}(overline{a+bsqrt{-1}})) 呢？注意到 DFT 的过程实际上是将每个 (a_i) 变为 (f(omega_n^i))，因此 ( ext{DFT}(overline{a})_i=sumlimits_{j=0}^noverline{a}_j(omega_n^i)^j)，又有：
- (overline{ab}=overline{a}·overline{b})
- (overline{omega_{n}^i}=omega_n^{n-i})
因此 ( ext{DFT}(overline{a})_i=sumlimits_{j=0}^noverline{a}_joverline{(omega_n^{n-i})^j}=overline{sumlimits_{j=0}^na_j(omega_n^{n-i})^j})，不难发现对于 (i=0) 而言，由于 (omega_{n}^n=1)，后面那东西就是 (overline{ ext{DFT}(a)_0})，而对于 (i e 0)，后面那东西就是 (overline{ ext{DFT}(a)_{n-i}})，因此我们求出 ( ext{DFT}(a+bsqrt{-1})) 后，把这个序列第 (1) 项至第 (n-1) 项 reverse 一下并把每一项的虚部变为其相反数即可得到 ( ext{DFT}(overline{a+bsqrt{-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,"
"[i==n+m]);
return 0;
}
各种多项式的 (n^2) 版本
qwq 大概是为了为下面子集卷积的高级操作做铺垫（？）
子集卷积的一些高级操作
首先在开始正题之前我们先得知道什么是集合幂级数，集合幂级数大概就是一个形如 (f(x)=sumlimits_{Ssubseteq 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) imes[x^T]G(x)) 的系数贡献到 (x^{Scup T},x^{Scap T},x^{Soplus T}) 的情况，对于这些卷积你都可以在我的快速沃尔什变换&快速莫比乌斯变换小记中找到，当然由于这里我们将着重探讨子集卷积，所以对于其他的卷积咱也大可不必花心思赘述了（
子集卷积，大概也是 ([x^S]F(x) imes[x^T]G(x)) 的系数贡献到 (x^{Scup T}) 的情况，不过有一个前提就是 (Scap T=varnothing) 时才能产生贡献，否则贡献为 (0)。那么怎么快速求出两个集合幂级数 (F,G) 的子集卷积呢？注意到关于集合的并，有一个性质 (|S|+|T|ge|Scup T|)，因此我们考虑将集合幂级数写成一个二元函数的形式（形式而已，并不是真正在集合幂级数后面加一维），即将每个 (a_Sx^S) 扩展成 (a_Sx^Sy^{|S|}) 的形式，这样每个集合幂级数可以写作 (sumlimits_{nge 0}sumlimits_{mge 0}a_{n,m}x^ny^m)，对于两个集合幂级数 (F,G)，假设它们的系数分别为 (f_{n,m}) 和 (g_{n,m})，那么我们重定义它们的乘法为 (F imes G=sumlimits_{nge 0}sumlimits_{mge 0}(sumlimits_{p|q=n}sumlimits_{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) 如果满足 (Scap T evarnothing)，那么其就一定不满足 (|Scap T|=|S|+|T|)，也就一定不会贡献到 (x^{Scup T}y^{|Scup T|}) 处，也就证明了上述做法的正确性。
那么这东西该怎么实现呢，我们考虑将集合幂级数 (F(x)) 扩充之后的样子写成 (sumlimits_{S}x^SG_S(y)) 的形式，其中 (G_S(y)) 是一个关于 (y) 的形式幂级数，这样不难发现我们的乘法运算可以视作：将每个集合幂级数前的系数视作一个形式幂级数，系数的乘法视作多项式的乘法，然后对它们做一遍 or 卷积。因此考虑按照套路 FWTor，然后对每个 (x^S)，将两个集合幂级数 (x^S) 前的系数进行多项式乘法，然后 IFWTor 回去即可。复杂度 (2^nn^2)。
注意到这里涉及形式幂级数的操作，因此我们也可以按照形式幂级数定义它们的逆、(ln)、(exp)。具体来说，对于集合幂级数 (F)，我们定义：
- 逆：定义 (F) 的逆为满足 (F imes G=epsilon) 的集合幂级数 (G)，其中 ( imes) 在这部分默认为子集卷积，(epsilon) 为满足 ([x^S]epsilon=[S=varnothing]) 的集合幂级数。
- (exp)：定义 (exp F) 为满足 (G=sumlimits_{nge 0}dfrac{F^n}{n!}) 的集合幂级数 (G)
- (ln)：定义 (ln F) 为满足 (exp G=F) 的集合幂级数 (G)
- (exp_{le k})：定义 (exp_{le k}) 为满足 (G=sumlimits_{n=0}^kdfrac{F^n}{n!}) 的集合幂级数
由于我们子集卷积可以看作是将集合幂级数都一遍 FWTor，然后每一项按照形式幂级数的套路卷起来，再 IFWTor，因此上述操作也可视作，将集合幂级数 FWTor 一遍，然后每一个 (x^S) 前的系数分别 ( ext{inv}/ln/exp/exp_{le k})。注意，由于集合幂级数的题 (n) 一般都很小，因此在这里就大可不必 NTT 实现多项式 ( ext{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，但这样会算重，因此考虑容斥，即
由于 Latex 上用中文太鸡肋了就换用了英文
注意到这东西可以写成子集卷积的形式，具体来说，即 (F) 为 (f) 的“生成集合幂级数”，(G(x)) 为满足 ([x^T]G(x)=[T ext{ is an independent set}]·(-1)^{|T|+1})，那么 (FG=F-1)，即 (F=dfrac{1}{1-G})，子集 inv 一波即可。
66. LOJ #154 集合划分计数
搞不懂模板题卡常是啥心态
记 (F) 为满足 ([x^S]F(x)=[xinmathcal F])，那么我们对 (F) 作 (exp_{le k}) 即可。
只不过此题过于卡常，所以交了 114514191981019260817998244353 发
卡常技巧：
- 求 ( ext{exp}_k) 时候，如果该多项式最低非零项次数 (>dfrac{n}{k}) 就不用快速幂了，直接将快速幂得到的数组每一项设为 (0) 即可。
- 在暴力多项式乘法时，可以找到两个多项式最低和最高的非零项 (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=sumlimits_{nge 0}G^n)，也就是 (exp) 里面去掉分母上的 (n!)，也就相当于对每个 (n) 个连通块的方案乘了个 (n!)，注意到这东西显然等于 (dfrac{1}{1-G})，因此再一波子集 (ln) 求出 (H)，最终答案就是 (H_U)。