题意
给出n,m,求
∑ni=1∑nj=1φ(ij)
。答案模1e9+7.
n<=105,m<=109
分析
看到n辣么小,考虑枚举n。
设
S(n,m)=∑i=1mφ(ni)
那么我们要求的就是
∑i=1nS(i,m)
首先有一个结论就是说若 |μ(n)|=1 那么有
φ(nk)=∑d|nφ(nd)φ(k)
证明:设 d=gcd(n,k) ,那么有
φ(nk)=d∗φ(ndk)=d∗φ(nd)φ(k)
把 ∑i|dφ(i)=d 带进去得
φ(nk)=φ(k)∑i|dφ(i)φ(nd)=φ(k)∑d|nφ(nd)
证毕。
把结论代进 S(n,m) 的式子里可以得到
S(n,m)=∑i=1mφ(ni)
=∑i=1m∑d|n,d|iφ(i)φ(nd)
=∑d|nφ(nd)∑i=1⌊md⌋φ(di)
=∑d|nφ(nd)S(d,⌊md⌋)
当 μ(n)=0 的时候,我们可以找到最大的 k 满足
当 n>1 时可以根据上式求 S(n,m) ,当 n=1 时就可以直接杜教筛了。
记得用map记忆化一下。
代码
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<map>
using namespace std;
typedef long long LL;
const int N=1000005;
const int MOD=1000000007;
int n,m,phi[N],low[N],tot,prime[N],s[N];
bool not_prime[N];
map<LL,int> ma;
void get_prime(int n)
{
phi[1]=low[1]=1;
for (int i=2;i<=n;i++)
{
if (!not_prime[i]) prime[++tot]=i,phi[i]=i-1,low[i]=i;
for (int j=1;j<=tot&&i*prime[j]<=n;j++)
{
not_prime[i*prime[j]]=1;
if (i%prime[j]==0)
{
phi[i*prime[j]]=phi[i]*prime[j];
low[i*prime[j]]=low[i];
break;
}
phi[i*prime[j]]=phi[i]*(prime[j]-1);
low[i*prime[j]]=low[i]*prime[j];
}
}
for (int i=1;i<=n;i++) s[i]=s[i-1]+phi[i],s[i]-=s[i]>=MOD?MOD:0;
}
LL point(int x,int y)
{
return (LL)(x-1)*m+(LL)y;
}
int solve(int n,int m)
{
if (m==0) return 0;
if (m==1) return phi[n];
if (n==1)
{
if (m<=1000000) return s[m];
if (ma[point(n,m)]) return ma[point(n,m)];
int ans=(LL)m*(m+1)/2%MOD;
for (int i=2,last;i<=m;i=last+1)
{
last=m/(m/i);
ans+=MOD-(LL)(last-i+1)*solve(1,m/i)%MOD;
ans-=ans>=MOD?MOD:0;
}
return ma[point(n,m)]=ans;
}
else
{
int k=low[n],ans=0;
if (ma[point(k,m)]) return (LL)ma[point(k,m)]*(n/k)%MOD;
for (int i=1;i*i<=k;i++)
if (k%i==0)
{
ans+=(LL)phi[k/i]*solve(i,m/i)%MOD;
ans-=ans>=MOD?MOD:0;
if (i*i==k) break;
ans+=(LL)phi[i]*solve(k/i,m/(k/i))%MOD;
ans-=ans>=MOD?MOD:0;
}
ma[point(k,m)]=ans;
ans=(LL)ans*(n/k)%MOD;
return ans;
}
}
int main()
{
get_prime(1000000);
scanf("%d%d",&n,&m);
if (n>m) swap(n,m);
int ans=0;
for (int i=1;i<=n;i++) ans+=solve(i,m),ans-=ans>=MOD?MOD:0;
printf("%d",ans);
return 0;
}