链接: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
输出描述:
一行一个整数表示输出结果一行一个整数表示输出结果
小A最近开始研究数论题了,这一次他随手写出来一个式子,∑ni=1∑mj=1gcd(i,j)2∑i=1n∑j=1mgcd(i,j)2,但是他发现他并不太会计算这个式子,你可以告诉他这个结果吗,答案可能会比较大,请模上1000000007。
输入描述:
一行两个正整数n,m一行两个正整数n,m
输出描述:
一行一个整数表示输出结果一行一个整数表示输出结果
输入:
2 2
输出:
7
2 2
输出:
7
1=<n,m<=1e6
解题思路:这题应该算是一题莫比乌斯反演的套路题了,感觉莫比乌斯真的好难啊,学了好久感觉也没懂,这题算是它的一个简单应用。
具体可以参考博客:https://blog.sengxian.com/algorithms/mobius-inversion-formula
具体可以参考博客:https://blog.sengxian.com/algorithms/mobius-inversion-formula
莫比乌斯反演经典套路:
现在有个积性函数f(n),设n<m,则:
![](https://img2018.cnblogs.com/blog/1396791/201904/1396791-20190423231455986-1115876032.png)
于是原来的式子就变成了求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; }