【bzoj3512】DZY Loves Math IV 杜教筛+记忆化搜索+欧拉函数

Description

给定n,m，求\(\sum_{i=1}^{n}\sum_{j=1}^{m}\varphi(ij)\)模10^9+7的值。

Input

仅一行，两个整数n,m。

Output

仅一行答案。

Sample Input

100000 1000000000

Sample Output

857275582

数据规模

1<=n<=10^5，1<=m<=109。

sol

%%%ranwen！！！

前置技能：

1. \(n=\sum_{d|n}\varphi(d)\)
2. \(\varphi(ij)=\varphi(\frac{i}{d})*\varphi(j)*d\quad[d=(i,j)]\)
3. \([\mu(x)=1||-1]\quad\varphi(\frac{x}{d})*\varphi(\frac{d}{e})=\varphi(\frac{x}{e})\quad[e|d]\)

证明：1不多说，2的意思就是让ij互质然后直接分解成两个phi，接着把gcd产生的倍数贡献乘回去，3因为x没有平方因子。除以d之后就不会有d 的因子了，所以与\(\frac{d}{e}\)互质，满足积性函数性质，乘起来即可。

解法：

首先观察数据范围可知，n的范围较小，可以进行枚举，而m的范围极大，已经超过了线性筛的范围。

我们考虑枚举i，然后推一波式子：

设\(s(n,m)=\sum_{i=1}^{m}\varphi(ni)\)

设w是n所有质因子一次方的乘积，v=n/w，则：

\(s(n,m)=v*\sum_{i=1}^{m}\varphi(iw)\)

设\(d=(i,w)\)，然后用公式2得：

\(s(n,m)=v*\sum_{i=1}^{m}\varphi(i)*\varphi(\frac{w}{d})*d\)

用公式1，得：

\(s(n,m)=v*\sum_{i=1}^{m}\varphi(i)*\varphi(\frac{w}{d})*\sum_{e|d}\varphi(\frac{d}{e})\)

用公式3，得：

\(s(n,m)=v*\sum_{i=1}^{m}\varphi(i)*\sum_{e|i,e|w}\varphi(\frac{w}{e})\)

因为 \(d=(i,w)\)，所以：

\(s(n,m)=v*\sum_{e|w}\varphi(\frac{w}{e})*\sum_{i=1}^{\lfloor\frac{m}{e}\rfloor}\varphi(ie)\)

根据\(s(n,m)\)的定义得：

\(s(n,m)=v*\sum_{e|w}\varphi(\frac{w}{e})*s(e,\frac{m}{e})\)

当n等于1的时候，可以直接使用杜教筛计算。

直接记忆化搜索即可。

因为每一步的m都是\(\lfloor\frac{x}{y}\rfloor\)的形式，所以可以使用下底分块法来计算。

时间复杂度\(O(n\sqrt{m}+n^{\frac{2}{3}})\)。

代码

#include <bits/stdc++.h>
using namespace std;
int n,m,tot,ans,phi[1000005],sum[1000005],pri[1000005],low[1000005],vis[1000005],P=1e9+7;
map<pair<int,int>,int>a;map<int,int>b;
int djs(int x)
{
	if(x<=1e6) return sum[x];
	if(b.count(x)) return b[x];
	int tp=1ll*x*(x+1)/2%P,last;
	for(int i=2;i<=x;i=last+1) last=x/(x/i),tp=(tp-1ll*(last-i+1)*djs(x/i)%P+P)%P;
	b.insert(make_pair(x,tp));return tp;
}
int solve(int x,int y)
{
	if(!x||!y) return 0;
	if(x==1) return djs(y);
	if(y==1) return phi[x];
	if(a.count(make_pair(x,y))) return a[make_pair(x,y)];
	int w=low[x],v=x/w,lim=floor(sqrt(w)+0.1),tp=0;
	for(int i=1;i<=lim;i++) if(w%i==0)
	{
		tp=(tp+1ll*phi[w/i]*solve(i,y/i)%P)%P;
		if(i!=w/i) i=w/i,tp=(tp+1ll*phi[w/i]*solve(i,y/i)%P)%P,i=w/i;
	}
	tp=1ll*tp*v%P;a.insert(make_pair(make_pair(x,y),tp));
	return tp;
}
int main()
{
	phi[1]=low[1]=sum[1]=1;
	for(int i=2;i<=1000000;sum[i]=(sum[i-1]+phi[i])%P,i++)
	{
		if(!vis[i]){pri[++tot]=i;phi[i]=i-1;low[i]=i;}
		for(int j=1;j<=tot&&i*pri[j]<=1000000;j++)
		{
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];low[i*pri[j]]=low[i];break;}
			phi[i*pri[j]]=phi[i]*(pri[j]-1),low[i*pri[j]]=low[i]*pri[j];
		}
	}
	scanf("%d%d",&n,&m);
	if(n>m) swap(n,m);
	for(int i=1;i<=n;i++) ans=(ans+solve(i,m))%P;
	printf("%d\n",ans);
}