BZOJ 3529 （ 莫比乌斯反演 + 树状数组维护狄利克雷卷积 ）

参考博客：https://www.lydsy.com/JudgeOnline/problem.php?id=3529

#include<bits/stdc++.h>

using namespace std;

const int maxn=1e5+10;

int mu[maxn],ans[maxn],c[maxn],vis[maxn],p[maxn],t[maxn],g[maxn];

struct F{int d,num;}f[maxn];  //f[i].d表示f[i].num的约数和

struct Q{int n,m,a,id;}q[maxn]; //q表示询问。

bool cmp1(Q a,Q b){
    return a.a<b.a;
}

bool cmp2(F a,F b){
    return a.d<b.d;
}

int quick_power(int a,int b){
    int res=1;
    while(b){
        if(b&1) res*=a;
        a*=a;
        b>>=1;
    }
    return res;
}

void init(){
    mu[1]=1;
    f[1].d=f[1].num=1;
    for(int i=2;i<maxn;i++){
        f[i].num=i;
        if(!vis[i]) mu[i]=-1,f[i].d=t[i]=i+1,g[i]=1,p[++p[0]]=i;
        //g  数组表示指数， t 数组表示 ∑p^k
        for(int j=1;j<=p[0]&&i*p[j]<maxn;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                g[i*p[j]]=g[i]+1;
                t[i*p[j]]=t[i]+quick_power(p[j],g[i]+1);
                f[i*p[j]].d=f[i].d/t[i]*t[i*p[j]];
                break;
            }else{
                mu[i*p[j]]=-mu[i];
                f[i*p[j]].d=f[i].d*f[p[j]].d;
                g[i*p[j]]=1;
                t[i*p[j]]=p[j]+1;
            }
        }
        //在纸上画一下，就会发现t与f之间的关系，因为线性筛的缘故，所以只要判断你最小的那个素数就可以了。
    }
}

int lowbit(int x){
    return x&(-x);
}

void add(int x,int val){
    for(int i=x;i<maxn;i+=lowbit(i)) c[i]+=val;
}

int query(int x){
    int sum=0;
    for(int i=x;i>=1;i-=lowbit(i)){
        sum+=c[i];
    }
    return sum;
}

int main(){
    init();
    int T;
    scanf("%d",&T);
    for(int i=1;i<=T;i++){
        scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a);
        q[i].id=i;
    }
    sort(q+1,q+1+T,cmp1);
    sort(f+1,f+maxn,cmp2);
    for(int now=1,i=1;i<=T;i++){
        for(;now<maxn&&f[now].d<=q[i].a;now++){
            for(int j=1;j*f[now].num<maxn;j++){
                add(j*f[now].num,mu[j]*f[now].d);
            }
        }
        int n=q[i].n,m=q[i].m;
        if(n>m) swap(n,m);
        for(int l=1,r;l<=n;l=r+1){
            r=min(n/(n/l),m/(m/l));
            ans[q[i].id]+=(n/l)*(m/l)*(query(r)-query(l-1));
        }
        ans[q[i].id]&=0x7fffffff;
    }
    for(int i=1;i<=T;i++){
        printf("%d\n",ans[i]);
    }
    return 0;
}
/*
2
4 4 3
10 10 5
*/