牛客小白月赛13-J小A的数学题 （莫比乌斯反演）

链接：https://ac.nowcoder.com/acm/contest/549/J
来源：牛客网

题目描述
小A最近开始研究数论题了，这一次他随手写出来一个式子，∑ni=1∑mj=1gcd(i,j)2∑i=1n∑j=1mgcd(i,j)2，但是他发现他并不太会计算这个式子，你可以告诉他这个结果吗，答案可能会比较大，请模上1000000007。
输入描述:
一行两个正整数n，m一行两个正整数n，m
输出描述:
一行一个整数表示输出结果一行一个整数表示输出结果

 

输入：
2 2
输出：
7

1=<n,m<=1e6

解题思路：这题应该算是一题莫比乌斯反演的套路题了，感觉莫比乌斯真的好难啊，学了好久感觉也没懂，这题算是它的一个简单应用。
具体可以参考博客：https://blog.sengxian.com/algorithms/mobius-inversion-formula

莫比乌斯反演经典套路：

现在有个积性函数f(n)，设n<m，则：

于是原来的式子就变成了求f∗μ了，再用上整数分块就可以快速搞定了。

自己推演了一遍：

代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#include<string>
#include<set>
#include<cmath>
#include<list>
#include<deque>
#include<cstdlib>
#include<bitset>
#include<stack>
#include<map>
#include<cstdio>
#include<queue>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
#define lson l,mid,rt<<1
#define rson mid+1,r,rt<<1|1
const int INF=0x3f3f3f3f;
const double PI=acos(-1.0);
const double eps=1e-6;
const int maxn=1000005;
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll lcm(ll a,ll b){return a/gcd(a,b)*b;}
const int mod=1e9+7;
const int dir[4][2]={{1,0},{-1,0},{0,1},{0,-1}};
const int N=1e6+10;
ll n,m,prime[N],mu[N],tot;
void getMu(){
    for(int i=1;i<=1e6+5;i++) prime[i]=1;
    mu[1]=1;
    for(int i=2;i<=1e6+5;i++){
        if(prime[i]){
            prime[++tot]=i;
            mu[i]=-1;
        }
            for(int j=1;j<=tot&&prime[j]*i<=1e6+5;j++){
                prime[i*prime[j]]=0;
                if(i%prime[j]==0){
                    mu[i*prime[j]]=0;
                    break;
                }else mu[i*prime[j]]=-mu[i];
            }
    }
}
int main(){
    cin>>n>>m;
    getMu();
    ll ans=0;
    for(ll i=1;i<=min(n,m);i++){
        ll tmp=0;
        for(ll j=i;j<=min(n,m);j+=i){
            tmp=(tmp+mu[j/i]*(n/j)*(m/j))%mod;
        }
        ans=(ans+tmp*i*i%mod)%mod;
    }
    cout<<ans<<endl;
    return 0;
}

整除分块优化：

#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
const int maxn=1e6+7;
const int mod=1e9+7;
ll n,m,prime[maxn],mu[maxn],sum[maxn],tot,ans;
void getMobius(int N){
    for(int i=1;i<=N;i++)prime[i]=1;
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(prime[i]){
            prime[tot++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<tot&&i*prime[j]<=N;j++){
            prime[i*prime[j]]=0;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
}
ll solve(ll a,ll b){
    ll res=0;
    for(int l=1,r;l<=min(a,b);l=r+1){
        r=min(a/(a/l),b/(b/l));
        res=(res+(sum[r]-sum[l-1])%mod*(a/l)%mod*(b/l)%mod)%mod;
    }
    return res;
}
int main(){
    scanf("%lld%lld",&n,&m);
    if(n>m) swap(n,m);
    getMobius(1e6);
    sum[1]=1;
    for(int i=2;i<=1e6;i++) sum[i]=sum[i-1]+mu[i];
    for(ll l=1,r;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        ll dd=(r*(r+1)*(2*r+1)/6-(l-1)*l*(2*l-1)/6)%mod;
        ans=(ans+dd*solve(n/l,m/l)%mod)%mod;
    }
    printf("%lld
",ans);
    return 0;
}