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

stargazer.

于 2020-01-19 19:07:51 发布

阅读量174

点赞数

分类专栏：莫比乌斯反演整除分块

本文链接：https://blog.csdn.net/qq_42555009/article/details/104044056

版权

莫比乌斯反演同时被 2 个专栏收录

20 篇文章 0 订阅

订阅专栏

整除分块

5 篇文章 0 订阅

订阅专栏

传送门

首先简单莫反后可以得到

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

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

考虑 $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_i\le 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<<'\n';
    }
}