【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！！！

前置技能：

(n=sum_{d|n}varphi(d))
(varphi(ij)=varphi(frac{i}{d})*varphi(j)*dquad[d=(i,j)])
([mu(x)=1||-1]quadvarphi(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}^{lfloorfrac{m}{e} floor}varphi(ie))

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

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

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

直接记忆化搜索即可。

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

时间复杂度(O(nsqrt{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
",ans);
}