【BZOJ3512】DZY Loves Math IV

题目大意

给出$n,m$，求$\sum _{i=1}^n\sum _{j=1}^m\varphi (ij)$

$1\le n\le 10^5,1\le m\le 10^9$，答案模$10^9+7$。

题解

前置知识：杜教筛

因为$n$比较小，所以我们可以直接枚举$n$。

设$f(n,m)=\sum _{i=1}^m\varphi (ni)$

将$n$分为若干个质数的乘积，$n=p_1^{a_1}{p}_2^{a_2}\cdots p_k^{a_k}$，令$p=p_1^{a_1-1}{p}_2^{a_2-1}\cdots p_k^{a_k-1}$，$q=n/p$

$$\sum _{i=1}^m\varphi (ni)=\sum _{i=1}^{m}p\varphi (qi)=\sum _{i=1}^{m}p\varphi (\frac{q}{\gcd (q,i)}\times i)\times \gcd (q,i)$$

因为$\varphi (i)$为积性函数，所以$\varphi (\frac{q}{\gcd (q,i)}\times i)=\varphi (\frac{q}{\gcd (q,i)})\times \varphi (i)$，原式变为

$$p\sum _{i=1}^m\varphi (\frac{q}{\gcd (q,i)})\times \varphi (i)\times \gcd (q,i)$$

由欧拉函数的性质，$n=\sum _{d|n}\varphi (\frac{n}{d})$，可得

$$p\sum _{i=1}^m\varphi (i)\varphi (\frac{q}{\gcd (q,i)})\sum _{j|i,j|q}\varphi (\frac{\gcd (q,i)}{j})$$

因为$\frac{q}{\gcd (q,i)}$与$\gcd (q,i)$互质（$q=p_1{p}_2\cdots p_k$，同一个质数要不不出现，要不出现在$\frac{q}{\gcd (q,i)}$或$\gcd (q,i)$，不可能两个都出现）而$\varphi$为积性函数，所以$\varphi (\frac{q}{\gcd (q,i)})\times \varphi (\frac{\gcd (q,i)}{j})=\varphi (\frac{q}{j})$，可得

$$p\sum _{i=1}^m\varphi (i)\sum _{j|i,j|q}\varphi (\frac{q}{j})$$

枚举$j$得

$$p\sum _{j|q}\varphi (\frac{q}{j})\sum _{i=1}^{\lfloor \frac{m}{j}\rfloor }\varphi (ij)$$

因为$f(n,m)=\sum _{i=1}^m\varphi (ni)$，所以

$$p\sum _{j|q}\varphi (\frac{q}{j})\times f(j,\lfloor \frac{m}{j}\rfloor )$$

于是我们就得到了递推式

$$f(n,m)=p\sum _{j|q}\varphi (\frac{q}{j})\times f(j,\lfloor \frac{m}{j}\rfloor )$$

预处理出前$n^{\frac{2}{3}}$个$\varphi (i)$并求前缀和，用杜教筛$\varphi (i)$的前缀和，利用$\varphi$的前缀和作差得出每个$\varphi (\frac{q}{j})$，然后就可以用杜教筛求$f(n,m)$了。

注意$f(n,m)$的一些特殊值。

• 当$n=1$时，$f(n,m)=\sum _{i=1}^m\varphi (i)$
• 当$m=1$时，$f(n,m)=\varphi (n)$

最终的答案为$\sum _{i=1}^{n}f(i,m)$。

code

#include<bits/stdc++.h>
using namespace std;
const int N=2000000;
const long long mod=1000000007;
int z[N+5],p[N+5],mn[N+5];
long long ans,ph[N+5],s[N+5];
map<long long,long long>sph,sum[100005];
void init(){
	ph[1]=1;
	for(int i=2;i<=N;i++){
		if(!z[i]){
			p[++p[0]]=i;
			ph[i]=i-1;mn[i]=i;
		}
		for(int j=1;j<=p[0]&&i*p[j]<=N;j++){
			z[i*p[j]]=1;
			mn[i*p[j]]=p[j];
			if(i%p[j]==0){
				ph[i*p[j]]=ph[i]*p[j]%mod;
				break;
			}
			ph[i*p[j]]=ph[i]*(p[j]-1)%mod;
		}
	}
	for(int i=1;i<=N;i++){
		s[i]=(s[i-1]+ph[i])%mod;
	}
}
long long gt(int n){
	if(n<=N) return s[n];
	if(sph.count(n)) return sph[n];
	long long re=0;
	for(int l=2,r;l<=n;l=r+1){
		r=min(n,n/(n/l));
		re=(re+1ll*(r-l+1)*gt(n/l)%mod)%mod;
	}
	re=(1ll*n*(n+1)/2-re+mod)%mod;
	return sph[n]=re;
}
long long djs(int n,int m){
	if(!m) return 0;
	if(n==1) return gt(m);
	if(m==1) return ph[n];
	if(sum[n].count(m)) return sum[n][m];
	long long re=0;
	vector<int>v;
	int p=1,q=1,t=n,l=0;
	while(t>1){
		int x=mn[t];
		q*=x;t/=x;
		v.push_back(x);++l;
		while(t%x==0){
			p*=x;t/=x;
		}
	}
	for(int i=0;i<1<<l;i++){
		int vt=1;
		for(int j=0;j<l;j++)
		if(i&(1<<j)) vt*=v[j];
		re=(re+ph[q/vt]*djs(vt,m/vt)%mod)%mod;
	}
	return sum[n][m]=re*p%mod;
}
int main()
{
	int n,m;
	init();
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n;i++){
		ans=(ans+djs(i,m))%mod;
	}
	printf("%lld",ans);
	return 0;
}