学习笔记 - 单位根反演 | Lucky_Glass's Blog
0%

学习笔记 - 单位根反演

听课的时候有一道题要单位根反演;
看起来是个轻量级的算法(说起来像真的一样),然后就来学一下


# 引理

定义 $\omega_n$ 表示 $n$ 次单位复根,即 $\omega_n^n=1$,对于任意正整数 $k$,有 $$n\cdot[n\mid k]=\sum_{i=0}^{n-1}\omega_n^{ik}$$

证明只需用到单位复根的性质。

如果 $n\mid k$,则 $\omega_n^{ik}=1$,

此时等比数列的公比不为 $\mathbf1$($n\mid k$ 的时候公比为 $1$,所以不能用等比数列求和),直接用等比数列求和得到

又因为 $\omega_n^{nk}=(\omega_n^n)^k=1$,所以


# 直接推导

求 $n$ 次多项式 $f(x)$ 的所有 $x^{ik}$ 项的系数之和($k\le n$)。形式化的,求

对式子化简:

就把 $k$ 个单位复根代入式子计算。当然一般来说式子的次数非常高,但是如果式子有比较快速的计算一次的方法,单位根反演就比较有用了。

# DFT推导

另外,还可以用 DFT 理解单位根反演,设 $a_i$ 为 $[x^i]f(x)$:

因为单位根 $k$ 次一循环,所以可以找到下图所示的规律:
png1

于是可以将矩阵“压缩”成 $k\times k$ 的:
png2

设 $b_t=\sum a_{ik+t}$,则压缩后的式子就是

不难发现上式就相当于是对 $b_i$ 的矩阵求了一个 DFT,那么我们只需要 IDFT 就可以得到 $b_i$ 了。

而直接推导求出的 $\sum a_{ik}$ 的公式则是 IDFT 的一个特例。


# 例题 - 生成树计数加强版(LOJ)

一道不这么板子题的例题

看到不进位加法,而且还是三进制的……只能想到一种做法,拆位。那么只用考虑边权只有 $[0,2]$ 的情况。

对于生成树计数,我们有一个很熟悉的算法叫做矩阵树定理。矩阵树定理最基础的应用就是求生成树的总方案数,但是还有一个比较常用的技巧:把矩阵 $M$ 中的整数换成其他类型——多项式。

如果存在边 $(u,v)$ 的权值为 $w$,则:

1
2
M[u][u]+=x^w,M[v][v]+=x^w;
M[u][v]-=x^w,M[v][u]-=x^w;

最后做完矩阵树定理,我们会得到一个次数很高的多项式 $F(x)$。$F(x)$ 第 $i$ 位的系数 $f_i$ 表示生成树的边权之和为 $i$ 的方案数,而我们要求的答案是

那么只需要分别求出 $F(x)$ 的模 $3$ 余 $0,1,2$ 的项的系数即可,直接套用单位根反演——做三次行列式,算出 $F(\omega_k^0),F(\omega_k^1),F(\omega_k^2)$,然后再 IDFT 就可以了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
/*Lucky_Glass*/
#include<cstdio>
#include<cassert>
#include<cstring>
#include<algorithm>
using namespace std;

const int MOD=1e9+7,VARA=500000003,VARB=541031193,N=105,M=N*N;
#define cint const int &
//(a+bi)^3=1 (mod MOD)
inline int Rint(int &r){
int b=1,c=getchar();r=0;
while(c<'0' || '9'<c) b=c=='-'? -1:b,c=getchar();
while('0'<=c && c<='9') r=(r<<1)+(r<<3)+(c^'0'),c=getchar();
return r*=b;
}
inline int Mul(cint a,cint b){return 1ll*a*b%MOD;}
inline int Pow(cint a,cint b){return b? Mul(Pow(Mul(a,a),b>>1),(b&1)? a:1):1;}
inline int Add(cint a,cint b){return a+b>=MOD? a+b-MOD:a+b;}
inline int Sub(cint a,cint b){return a-b<0? a-b+MOD:a-b;}
#define square(a) Mul(a,a)
const int INV3=Pow(3,MOD-2);

struct NCOMPLEX{
#define cnc const NCOMPLEX &
#define oper(typ) inline friend typ operator
int a,b;
NCOMPLEX(){}
NCOMPLEX(cint A,cint B):a(A),b(B){}
oper(NCOMPLEX) *(cnc A,cnc B){return NCOMPLEX(Sub(Mul(A.a,B.a),Mul(A.b,B.b)),Add(Mul(A.a,B.b),Mul(B.a,A.b)));}
oper(NCOMPLEX) *(cnc A,cint B){return NCOMPLEX(Mul(A.a,B),Mul(A.b,B));}
oper(NCOMPLEX) +(cnc A,cnc B){return NCOMPLEX(Add(A.a,B.a),Add(A.b,B.b));}
oper(NCOMPLEX) -(cnc A,cnc B){return NCOMPLEX(Sub(A.a,B.a),Sub(A.b,B.b));}
oper(NCOMPLEX) -(cnc A){return NCOMPLEX(Sub(0,A.a),Sub(0,A.b));}
inline void operator +=(cnc A){(*this)=(*this)+A;}
inline void operator -=(cnc A){(*this)=(*this)-A;}
inline void operator *=(cnc A){(*this)=(*this)*A;}
inline friend NCOMPLEX inv(cnc A){
return NCOMPLEX(A.a,Sub(0,A.b))
*Pow(Add(square(A.a),square(A.b)),MOD-2);
}
#undef cnc
#undef oper
};

const NCOMPLEX con[3]={NCOMPLEX(1,0),NCOMPLEX(VARA,VARB),NCOMPLEX(VARA,VARB)*NCOMPLEX(VARA,VARB)};

int n,m;
int edg[M][3];
NCOMPLEX mat[N][N];

NCOMPLEX Det(){
NCOMPLEX ans(1,0);
for(int i=1;i<n;i++){
for(int j=i;j<n;j++)
if(mat[j][i].a || mat[j][i].b){
if(i==j) break;
swap(mat[i],mat[j]),ans=-ans;
break;
}
ans=ans*mat[i][i];
NCOMPLEX tmp=inv(mat[i][i]);
for(int j=i+1;j<n;j++){
NCOMPLEX now=mat[j][i]*tmp;
for(int k=i;k<n;k++) mat[j][k]-=now*mat[i][k];
}
}
return ans;
}
inline void IDFT(NCOMPLEX *p){
int i;
NCOMPLEX tmp[3];
tmp[0]=p[0],tmp[1]=p[1],tmp[2]=p[2];
p[0]=tmp[0]+tmp[1]+tmp[2];
p[1]=tmp[0]+tmp[1]*inv(con[1])+tmp[2]*inv(con[2]);
p[2]=tmp[0]+tmp[1]*inv(con[2])+tmp[2]*inv(con[1]);
for(i=0;i<3;++i) p[i]=p[i]*INV3;
}
int Calc(){
NCOMPLEX res[3];
for(int i=0;i<3;i++){
memset(mat,0,sizeof mat);
NCOMPLEX w[3]={con[0],con[i],con[i*2%3]};
for(int j=1;j<=m;j++){
int typ=edg[j][2]%3;
mat[edg[j][0]][edg[j][0]]+=w[typ];
mat[edg[j][1]][edg[j][1]]+=w[typ];
mat[edg[j][0]][edg[j][1]]-=w[typ];
mat[edg[j][1]][edg[j][0]]-=w[typ];
}
res[i]=Det();
}
IDFT(res);
return Add(res[1].a,Add(res[2].a,res[2].a));
}
int main(){
freopen("sum.in","r",stdin);
freopen("sum.out","w",stdout);
Rint(n),Rint(m);
for(int i=1;i<=m;i++)
Rint(edg[i][0]),Rint(edg[i][1]),Rint(edg[i][2]);
int ans=0;
for(int tim=1;;tim=Mul(tim,3)){
ans=Add(ans,Mul(tim,Calc()));
bool exi=false;
for(int i=1;i<=m;i++)
if(edg[i][2]/=3)
exi=true;
if(!exi) break;
}
printf("%d\n",ans);
return 0;
}

THE END

Thanks for reading!

> Linked 星星在唱歌-网易云