Description
给定n,m，求
模10^9+7的值。
Solution
设 (S(n,m)) 表示 (sum_{i=1}^{m}phi(n*i))
(Ans=sum_{i=1}^{n}S(i,m))
(S(n,m)=sum_{i=1}^{m}phi(n*i))
如果 (mu(n)!=0)
则有 (sum_{i=1}^{m}phi(frac{n}{gcd(i,n)})*phi(i)*gcd(i,n)) (因为要保证除完(gcd)之后,两数不含相同的质因子,所以 (mu(n)!=0))
(sum_{i=1}^{m}phi(frac{n}{gcd(i,n)})*phi(i)*sum_{d|i,d|n}phi(d))
因为第一项和第三项是互质的 , 所以可以合并.
(sum_{i=1}^{m}phi(i)*sum_{d|n,d|i}phi(frac{n}{d}))
(sum_{d=1}^{n}phi(frac{n}{d})*sum_{i=1}^{frac{m}{d}}phi(d*i))
(sum_{d=1}^{n}phi(frac{n}{d})*S(d,lfloorfrac{m}{d}
floor))
递归处理即可
如果 (mu(n)=0)
我们直接提出 (n) 的多出的质因子之积 (a),使得 (mu(frac{n}{a})!=0)
那么 (S(n,m)=sum_{i=1}^{m}phi(n*i)) 中也可以提出 (a) 了,因为相同的质因子只会被算一次
根据定义式 (phi(n)=n*Pi p_i),所以 (a) 唯一的贡献就是使前面的 (n) 乘了个 (a)
(a*S(n,m)=sum_{i=1}^{m}phi(frac{n}{a}*i)=a*S(frac{n}{a},m))
递归处理即可
边界条件 (m=1) 时,结果为 (phi(n)), (n=1) 时,跑一个杜教筛就行了
#include<bits/stdc++.h>
using namespace std;
const int N=1e5+10,mod=1e9+7;
int n,prime[N],num=0,m,phi[N],mu[N],pre[N],la[N],s[N];bool vis[N];
void priwork(){
phi[1]=s[1]=1;
for(int to,i=2;i<N;i++){
if(!vis[i])prime[++num]=i,phi[i]=i-1,mu[i]=-1,pre[i]=i;
for(int j=1;j<=num && i*prime[j]<N;j++){
vis[to=i*prime[j]]=1;pre[to]=prime[j];
if(i%prime[j])phi[to]=phi[i]*(prime[j]-1),mu[to]=-mu[i];
else {phi[to]=phi[i]*prime[j];break;}
}
s[i]=(s[i-1]+phi[i])%mod;
}
for(int i=2;i<=n;i++){
if(mu[i])continue;
int last=0,x=i;la[i]=1;
while(x>1){
if(pre[x]==last)la[i]*=pre[x];
last=pre[x];x/=pre[x];
}
}
}
map<int,int>S[N],T;
inline int calc(int n){
if(n<N)return s[n];
if(T.find(n)!=T.end())return T[n];
int ret=(1ll*n*(n+1)>>1)%mod;
for(int i=2,r;i<=n;i=r+1){
r=n/(n/i);
ret=(ret-1ll*calc(n/i)*(r-i+1))%mod;
}
if(ret<0)ret+=mod;
return T[n]=ret;
}
inline int solve(int n,int m){
if(m==1)return phi[n];
if(n==1)return calc(m);
if(S[n].find(m)!=S[n].end())return S[n][m];
if(!mu[n])return 1ll*la[n]*solve(n/la[n],m)%mod;
int ret=0,lim=min(m,(int)sqrt(n));
for(int i=1;i<=lim;i++){
if(n%i==0){
if(i*i!=n)ret=(ret+1ll*phi[n/i]*solve(i,m/i)+1ll*phi[i]*solve(n/i,m/(n/i)))%mod;
else ret=(ret+1ll*phi[n/i]*solve(i,m/i))%mod;
}
}
return S[n][m]=ret;
}
int main(){
freopen("pp.in","r",stdin);
freopen("pp.out","w",stdout);
cin>>n>>m;
priwork();
int ans=0;
for(int i=1;i<=n;i++)ans=(ans+solve(i,m))%mod;
printf("%d
",ans);
return 0;
}