莫比乌斯反演

前言

反演实际上是利用莫比乌斯函数实现 和函数与原函数之间的转换

我们需要求原函数，而实际上和函数是容易求得的

常用反演规律

表示..为n的方案数

表示..为n的倍数的方案数

反演

例题

公约数的和 可以当作模版题 反演 这一题还可以对式子进行变换，对于j=i+1这种和式超级难受

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void solve(){
    // de("\n\n");
    cin>>n;
    getMu();
    int ans=0;
    //先算出F[i]的值，公倍数为i的倍数的方案数
    rep(i,1,n){
        F[i]=0;
        //枚举倍数，公倍数为j 乘法原理计算方案
        for(int j=i;j<=n;j+=i){
            F[i]+=(n/i-j/i);
        }
    }
    rep(i,1,n){
        f[i]=0;
        //枚举倍数d 根据反演式计算F[i]
        for(int d=i;d<=n;d+=i){
            f[i]+=F[d]*mu[d/i];
        }
        ans+=f[i]*i;
    }
    cout<<ans<<endl;
}

式子变换之后的写法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
void solve(){
    // de("\n\n");
    cin>>n;
    init();
    int ans=0;
    rep(i,1,n)  F[i]=(n/i)*(n/i);
    rep(i,1,n){
        f[i]=0;
        for(int d=i;d<=n;d+=i){
            f[i]+=F[d]*mu[d/i];
        }
        ans+=f[i]*i;
    }
    ans-=(1+n)*n/2;
    ans/=2;
    cout<<ans<<endl;
}

能量采集

序列

f(n)表示gcd为n的方案数，F(n)表示gcd为n的倍数的方案数，所求答案即为f(1)

跳蚤

相同的方法反演，使用乘法原理就可以算出gcd为k的倍数的方案数

公约数之和加强版

LCMS

和式可以变换为

上面的式子是比较好算的，交叉相乘为和的乘积，利用反演将原函数用上述和函数表示

Unknown environment 'alig' \begin{aligned} f(n)=\sum_{n|d}F(d)\mu(\frac{d}{n})\\ ans=\sum_{d=1}^n\frac{1}{d}f(d) \end{alig} 最小公倍数之和

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
void solve(){
    // de("\n\n");
    cin>>n;
    init(N-10);
    int sum=0,mx=0;
    rep(i,1,n){
        cin>>a[i];
        cnt[a[i]]++;
        mx=max(mx,a[i]);
        sum+=a[i];
        sum%=p;
    }
    n=mx;
    rep(i,1,n){
        //两个因数都是i的倍数
        //枚举i的倍数为j
        for(int j=i;j<=n;j+=i){
            F[i]+=cnt[j]*j;
            F[i]%=p;
        }
        F[i]*=F[i];
        F[i]%=p;
    }
    int ans=0;
    rep(i,1,n){
        //公约数为i的方案数，计算反演式
        for(int d=i;d<=n;d+=i){
            //mu[d/i]可能为负
            f[i]=(f[i]+F[d]*mu[d/i]%p+p)%p;
        }
        ans=(ans+f[i]*inv[i]%p+p)%p;
    }
    ans=(ans-sum+p)%p;
    ans=ans*inv[2]%p;
    cout<<ans<<endl;
}

分块优化

ZAP Queries

显然是非常好计算的，但是在多case的情况下，如果不进行优化很容易超时

优化前

1
2
3
4
for(int i=d;i<=min(a,b);i+=d){
    F[i]=(a/i)*(b/i);
    ans+=F[i]*mu[i/d];
}

发现F[i]其实可以分块来做，在一定区间内，F[i]的值是不变的，统计莫比乌斯函数区间和即可

t时枚举量，实际上除数是a/d,b/d

1
2
3
4
for(int l=1,r;l<=min(a/d,b/d);l=r+1){
    r=min(a/d/(a/d/l),b/d/(b/d/l));
    ans+=(a/(l*d))*(b/(l*d))*(mu[r]-mu[l-1]);
}
1
2
3
4
5
a/=d,b/=d;
for(int l=1,r;l<=min(a,b);l=r+1){
    int r=min(a/(a/l),b/(b/l));
    ans+=(a/l)*(b/l)*(sum[r]-sum[l-1]);
}

数表

不需要确保枚举量是d的倍数，使用树状数组求前面一部分前缀和，只有d的倍数才会被添加到树状数组中

离线查找答案

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
//记录1e5以下的数据和约数和
struct node{
    int sum,val;
    bool operator<(const node& o)const{
        return sum<o.sum;
    }
}b[N];
struct QUERY{
    int n,m,a,id;
    bool operator<(const QUERY& o)const{
        return a<o.a;
    }
}Q[N];
int solve2(int n,int m){
    int t=min(n,m),ans=0;
    for(int l=1,r;l<=t;l=r+1){
        r=min(n/(n/l),m/(m/l));
        //query只会查到[l,r]区间内d的倍数的位置
        ans+=(query(r)-query(l-1))*(n/l)*(m/l);
    }
    return ans;
}
int ans[N];
void solve(){
    // de("\n\n");
    int q;
    cin>>q;
    rep(i,1,q){
        cin>>Q[i].n>>Q[i].m>>Q[i].a;
        Q[i].id=i;
    }
    sort(b+1,b+100001);
    sort(Q+1,Q+1+q);
    for(int i=1,j=1;i<=q;i++){
        while(b[j].sum<=Q[i].a and j<=1e5){
            //将满足条件的约数和的倍数塞入到树状数组中
            for(int d=b[j].val;d<=1e5;d+=b[j].val){
                add(d,b[j].sum*mu[d/b[j].val]);  
            }
            j++;
        }
        //分块查找
        ans[Q[i].id]=solve2(Q[i].n,Q[i].m);
    }
    rep(i,1,q){
        if(ans[i]<0) ans[i]+=pow(2,31);
        cout<<ans[i]<<endl;
    }
}