[BZOJ3512] DZY Loves Math IV

传送门
昨天花了好久好像是看懂了，那今天早上尝试自己推一遍柿子 顺便水了一篇博客

\(\bf {Description}\)

求$$\sum_{i=1}^n \sum_{i=1}^m\varphi(ij)$$

\(1 \leq n \leq 10^5 , 1 \leq m \leq 10^9\)，模 \(10^9+7\)

\(\bf Solution\)

先给个结论，当 \(|\mu(n)|=1\) 时，有

\[\varphi(ni)=\varphi(i) \sum_{d|(n,i)}\varphi(\frac{n}{d}) \]

简单证明一下（想了一节晚自习才想通呢QAQ）
（gcd的表示有点乱，不要在意）

\[\varphi(ni)=\varphi(n) \cdot \varphi(\frac{i}{gcd(n,i)}) \cdot gcd(n,i) \]

\[=\varphi(i) \cdot \varphi(\frac{n}{gcd(n,i)}) \cdot gcd(n,i) \]

\[=\varphi(i) \cdot \varphi(\frac{n}{gcd(n,i)}) \cdot \sum_{d|(n,i)}\varphi(d) \]

\[=\varphi(i) \cdot \varphi(\frac{n}{gcd(n,i)}) \cdot \sum_{d|(n,i)}\varphi(\frac{gcd(n,i)}{d}) \]

\[=\varphi(i) \sum_{d|(n,i)}\varphi(\frac{n}{d}) \]

推的时候要注意 \(n\) 的性质，还要熟悉 \(\varphi\)，不然就推不粗来QAQ

然后令 \(S(n,m)=\sum_{i=1}^m \varphi(ni)\)，令 \(P\) 为 \(n\) 所有质因子的乘积，\(Q=\dfrac{n}{P}\)，然后可得

\[S(n,m)=Q \cdot \sum_{i=1}^m \varphi(i) \sum_{d|(P,i)} \varphi(\frac{P}{d}) \]

先枚举因数，得到

\[S(n,m)=Q \cdot \sum_{d|P}\varphi(\frac{P}{d}) \sum_{i=1}^{\lfloor \frac{m}{d} \rfloor} \varphi(di) \]

\[=Q \cdot \sum_{d|P} \varphi(\frac{P}{d}) \ S(d,\left \lfloor \frac{m}{d} \right \rfloor) \]

然后就阔以记忆化啦？
当 \(n=1\) 的时候就是 \(\varphi\) 的前缀和，杜教筛板子啦？
复杂度玄学？不会证啦。
代码？还没写出来啦。
写出来再贴吧。

---

写出来啦，十分暴力QAQ

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
#include<bits/stdc++.h>
#include<tr1/unordered_map>
#define LL long long
#define re register
#define fr(i,x,y) for(int i=(x);i<=(y);i++)
#define rf(i,x,y) for(int i=(x);i>=(y);i--)
#define frl(i,x,y) for(int i=(x);i<(y);i++)
using namespace std;
using namespace tr1;
const int N=2000002;
const int INF=2147483647;
const int p=1e9+7;
int n,m;

void read(int &x){ scanf("%d",&x); }

void Add(int &x,int y){
	x+=y;
	while(x<0) x+=p;
	while(x>=p) x-=p;
}

int phi[N];
int pri[N/10],b[N],L;
void init(){
	phi[1]=1;
	frl(i,2,N){
		if (!b[i]) pri[++L]=i,phi[i]=i-1;
		for(int j=1;j<=L&&i*pri[j]<N;j++){
			b[i*pri[j]]=1;
			if (i%pri[j]==0){
				phi[i*pri[j]]=phi[i]*pri[j];
				break;
			}
			phi[i*pri[j]]=phi[i]*(pri[j]-1);
		}
	}
	frl(i,2,N) Add(phi[i],phi[i-1]);
}

map<int,int> s,mp[N];
int S(int n,int m){
	//cout<<n<<' '<<m<<endl;
	if (n==0||m==0) return 0;
	if (n==1){
		if (m<N) return phi[m];
		if (s.count(m)) return s[m];
		int ans=(1LL*m*(m+1)/2)%p;
		for(re int L=2,r=2;L<=m;L=r+1)
		 r=m/(m/L),Add(ans,-1LL*S(n,m/L)*(r-L+1)%p);
		return s[m]=ans;
	}
	if (mp[n].count(m)) return mp[n][m];
	int ans=0;
	for(re int i=1;i*i<=n;i++)
	 if (n%i==0){
	 	Add(ans,1LL*(phi[n/i]-phi[n/i-1])*S(i,m/i)%p);
	 	if (n/i!=i) Add(ans,1LL*(phi[i]-phi[i-1])*S(n/i,m/(n/i))%p);
	 }
	return mp[n][m]=ans;
}

int main(){
	//freopen("1.in","r",stdin);
	init();
	read(n);read(m);
	int ans=0;
	fr(xx,1,n){
		int s=1,x=xx;
		for(int i=2;i*i<=x;i++)
		 if (x%i==0){
		 	s*=i;
		 	while(x%i==0) x/=i;
		 }
		s*=x;
		Add(ans,1LL*S(s,m)*(xx/s)%p);
		//cout<<S(s,m)*(xx/s)<<endl;
	}
	cout<<ans<<endl;
	//cout<<S(s,m)<<endl;
	//printf("%lld\n",1LL*S(s,m)*(n/s)%p);
	return 0;
}