单位根反演

Posted wangziqiang123

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了单位根反演相关的知识，希望对你有一定的参考价值。

先扔个题引入吧(好像大家学这玩意第一题都写的是这个)：

给定(n,s,a_0,a_1,a_2,a_3)，求： [ sum_{i=0}^{n}binom{n}{i}s^ia_{ibmod 4} ] 其中(nleqslant 1e18)。

题目链接。

简单变换一下，我们需要算这个： [ sum_{r=0}^3a_rsum_{i=0}^{n}binom{n}{i}s^i[ibmod 4=r] ] 设生成函数(f(x)=sum_{i=0}^{n}binom{n}{i}s^ix^i=(sx+1)^n)。

那么我们现在要分别求(f)的(bmod 4=i)的项之和。

我们引入单位根反演，就是说在(rm ntt)结束的时候证明一般是用到了这个式子： [ frac{1}{n}sum_{i=0}^{n-1}omega_{n}^{ik}=[n|k] ] 当(n|k)的时候显然，每一项都是(1)，否则可以利用等比数列求和得到。

我们把单位根带进前面的生成函数，假设生成函数系数为(a_i)： [ sum_{r=0}^{m-1}f(omega_m^r)=sum_{i=0}^{n}a_isum_{t=0}^{m-1}omega_m^{it}=msum_{i=0}^{n}a_i[ibmod m=0] ] 那么我们就得到了(f)的(bmod 4=0)的项之和，至于(bmod 4)为其他值的和我们可以通过把(f)右移若干位得到。

那么最终答案就是： [ frac{1}{4}sum_{r=0}^3a_{(-r)bmod 4}sum_{t=0}^3f(omega_4^t)omega_4^{rt} ] 代码对着式子抄一遍就行了，我就不放了

在看看这个题。

给定(n,k)，求： [ sum_{i=0}^{n}binom{n}{i}F_i[k|i] ] 其中(F_i)为斐波那契数列，(nleqslant 1e18,kleqslant 2e4)。

我们把这玩意搞成矩阵形式： [ sum_{i=0}^{n}binom{n}{i}A^i[k|i] ] 其中(A)是斐波那契的转移矩阵。

那么生成函数： [ f(x)=sum_{i=0}^{n}binom{n}{i}A^ix^i=(Ax+I)^n ] (I)是单位矩阵。

所以套一下上面的式子就写完了： [ frac{1}{k}sum_{i=0}^{k-1}f(w_k^i) ]

code:

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
114
115
116
117
118
119
120
121
122
123

#include<bits/stdc++.h>
using namespace std;

void (int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}

void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('n');}

#define lf double
#define ll long long 

#define pii pair<int,int >
#define vec vector<int >

#define pb push_back
#define mp make_pair
#define fr first
#define sc second

#define FOR(i,l,r) for(int i=l,i##_r=r;i<=i##_r;i++)

const int maxn = 1e6+10;
const int inf = 1e9;
const lf eps = 1e-8;

int k,p,g,w[maxn];ll n;
vector<int > r;

int qpow(int a,int x) {
    int res=1;
    for(;x;x>>=1,a=1ll*a*a%p) if(x&1) res=1ll*res*a%p;
    return res;
}

void gen() {
    int e=p-1;r.clear();
    for(int i=2;i*i<=e;i++) {
        if(e%i) continue;
        r.pb(i);while(!(e%i)) e/=i;
    }if(e!=1) r.pb(e);
    for(g=2;;g++) {
        int bo=1;
        for(int x=0;x<r.size();x++)
            if(qpow(g,(p-1)/r[x])==1) {bo=0;break;}
        if(bo) break;
    }
}

int add(int x,int y) {x+=y;return x>=p?x-p:x;}

struct matrix {
    int a[2][2];

    void clear() {memset(a,0,sizeof a);}

    matrix operator * (const matrix &r) const {
        matrix res;res.clear();
        
        //     for(int j=0;j<=1;j++)
        //         for(int k=0;k<=1;k++)
        //             res.a[i][j]=add(res.a[i][j],1ll*a[i][k]*r.a[k][j]%p);
        res.a[0][0]=add(1ll*a[0][0]*r.a[0][0]%p,1ll*a[0][1]*r.a[1][0]%p);
        res.a[0][1]=add(1ll*a[0][0]*r.a[0][1]%p,1ll*a[0][1]*r.a[1][1]%p);
        res.a[1][0]=add(1ll*a[1][0]*r.a[0][0]%p,1ll*a[1][1]*r.a[1][0]%p);
        res.a[1][1]=add(1ll*a[1][0]*r.a[0][1]%p,1ll*a[1][1]*r.a[1][1]%p);
        return res;
    }

    matrix operator * (const int &r) const {
        matrix res;res.clear();
        
        //     for(int j=0;j<=1;j++)
        //         res.a[i][j]=1ll*a[i][j]*r%p;
        res.a[0][0]=1ll*a[0][0]*r%p;
        res.a[0][1]=1ll*a[0][1]*r%p;
        res.a[1][0]=1ll*a[1][0]*r%p;
        res.a[1][1]=1ll*a[1][1]*r%p;
        return res;
    }

    matrix operator + (const matrix &r) const {
        matrix res;res.clear();
        
        //     for(int j=0;j<=1;j++)
        //         res.a[i][j]=add(a[i][j],r.a[i][j]);
        res.a[0][0]=add(a[0][0],r.a[0][0]);
        res.a[0][1]=add(a[0][1],r.a[0][1]);
        res.a[1][0]=add(a[1][0],r.a[1][0]);
        res.a[1][1]=add(a[1][1],r.a[1][1]);
        return res;
    }
}I,A;

matrix qpow(matrix a,ll x) {
    matrix res=I;
    for(;x;x>>=1,a=a*a) if(x&1) res=res*a;
    return res;
}

void solve() {
    // int st=clock();
    scanf("%lld",&n),read(k),read(p);gen();
    I.a[0][0]=I.a[1][1]=1;w[0]=1,w[1]=qpow(g,(p-1)/k);
    A.a[0][0]=A.a[0][1]=A.a[1][0]=1;
    for(int i=2;i<k;i++) w[i]=1ll*w[i-1]*w[1]%p;
    matrix ans;ans.clear();
    for(int i=0;i<k;i++) ans=ans+qpow(A*w[i]+I,n);
    write(1ll*ans.a[0][0]*qpow(k,p-2)%p);
    // cerr<<(lf)(clock()-st)/1e3<<endl;
}

int main() {
    int t;read(t);while(t--) solve();
    return 0;
}

原文:大专栏  单位根反演


以上是关于单位根反演的主要内容，如果未能解决你的问题，请参考以下文章

单位根反演和循环卷积

loj#6485. LJJ 学二项式定理（单位根反演）

单纯看懂公式的单位根反演

单纯看懂公式的单位根反演

单位根反演

UOJ450 复读机 生成函数单位根反演