【BZOJ 3309】DZY Loves Math（莫比乌斯反演）

传送门

首先简单莫反后可以得到

$ans=sum_{T}frac a Tfrac b Tsum_{d|T}mu(d)f(frac{T}{d})$
$=sum_{T}frac a Tfrac b T g(T)$

考虑怎么求$g$
一个显然的方法是$nlogn$枚举，但是过不了

考虑$g(t)=sum_{d|t}f(d)mu(t/d)$
$t=prod_{i=1}^kp_i^{b_i},t/d=prod_{i}p_i^{a_i}$
由于$mu$的性质，所以所有有贡献的$t/d$中$a_ile 1$

那么此时$f(d)=f(t)$或$f(t)-1$

考虑当$t$分解出的$b_i$不是都相等时
那么就一定存在某个$a_i$的取值不影响$f$
但是由于对于$mu$的贡献是$+1/-1$
所以会恰好抵消

而当$b_i$都相等时

除了当$a_i$都为1时$f$的值要少1
如果不看这个$f$少1的话和应该为0的
所以减去这个差量也就是贡献为$-(-1)^k=(-1)^{k+1}$

所以只有$t=(prod_{i=1}^kp_i)^r$时
$g(t)=(-1)^{k+1}$
否则$g(t)=0$

这个东西可以线筛了
利用线筛是用最小质因子筛，维护最小质因子的次数以及除去最小质因子后的值就可以了
具体实现可以看代码

#include<bits/stdc++.h>
using namespace std;
#define cs const
#define re register
#define pb push_back
#define pii pair<int,int>
#define fi first
#define ll long long
#define se second
#define bg begin
cs int RLEN=(1<<20)+1;
inline char gc(){
    static char ibuf[RLEN],*ib,*ob;
    (ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
    return (ib==ob)?EOF:*ib++;
}
inline int read(){
    char ch=gc();
    int res=0;bool f=1;
    while(!isdigit(ch))f^=ch=='-',ch=gc();
    while(isdigit(ch))res=res*10+(ch^48),ch=gc();
    return f?res:-res;
}
template<class tp>inline void chemx(tp&a,tp b){a<b?a=b:0;}
template<class tp>inline void chemn(tp&a,tp b){a>b?a=b:0;}
cs int N=10000007;
int g[N],pw[N],ecp[N],pr[N/10+1];
int tot;
bitset<N> vis;
inline void init(cs int len=N-7){
    for(int i=2;i<=len;i++){
        if(!vis[i])pr[++tot]=i,g[i]=pw[i]=ecp[i]=1;
        for(int j=1;j<=tot&&i*pr[j]<=len;j++){
            int p=i*pr[j];
            vis[p]=1;
            if(i%pr[j]==0){
                ecp[p]=ecp[i];
                pw[p]=pw[i]+1;
                if(ecp[p]==1)g[p]=1;
                else g[p]=(pw[ecp[p]]==pw[p])?-g[ecp[p]]:0;
                break;
            }
            ecp[p]=i,pw[p]=1,g[p]=(pw[i]==1?-g[i]:0);
        }
    }
    for(int i=2;i<=len;i++)g[i]+=g[i-1];
}
int main(){
    #ifdef Stargazer
    freopen("lx.in","r",stdin);
    #endif
    init();
    int T=read();
    while(T--){
        int a=read(),b=read();
        ll ret=0;
        for(int l=min(a,b),i=1,nxt;i<=l;i=nxt+1){
            nxt=min(a/(a/i),b/(b/i));
            ret+=1ll*(a/i)*(b/i)*(g[nxt]-g[i-1]);
        }
        cout<<ret<<'
';
    }
}