[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

Solution

神题。。

考虑(n)比较小，试着枚举(n)，处理(m)。

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

那么考虑拆解(n)，设(n=prod_{i=1}^{k}a_i^{p_i})，对于(n)的每一个质因子，提一个出来，剩下的可以直接乘上去，即：

设(y=prod_{i=1}^k a_i^{p_i-1})，(w=prod_{i=1}^ka_i)，式子就变成：

[egin{align} s(n,m)&=ycdotsum_{i=1}^{m}varphi(wi)\ &=ycdot sum_{i=1}^{m}varphi(frac{w}{(i,w)})cdot varphi(i)cdot (i,w) end{align} ]

其中((i,j))表示(gcd(i,j))。

然后由于：

[sum_{d|n}varphi(d)=n ]

可以反过来把后面的(gcd)化一下，即：

[egin{align} s(n,m)&=ycdot sum_{i=1}^{m}varphi(frac{w}{(i,w)})cdotvarphi(i)cdotsum_{d|i&d|w}varphi(d) \ &=ycdot sum_{i=1}^{m}varphi(i)cdotsum_{d|i&d|w}varphi( frac{(i,w)}{d}) varphi(frac{w}{(i,w)})\ &=ycdot sum_{i=1}^{m}varphi(i)cdotsum_{d|i&d|w}varphi( frac{w}{d})\ &=ycdot sum_{d|w}varphi(frac{w}{d})sum_{i=1}^{lfloorfrac{m}{d} floor}varphi(id)\ &=ycdot sum_{d|w}varphi(frac{w}{d})s(d,lfloorfrac{m}{d} floor) end{align} ]

然后这个东西变成了一个递归的形式，直接照着式子爆算一下就好了，记得记忆化。

当(n=1)时这个东西就变成了(sum_{i=1}^{m}varphi(i))，杜教筛一下就好了。

#include<bits/stdc++.h>
//#include<tr1/unordered_map>
using namespace std;
//using namespace std :: tr1;

#define ll long long 
//#define ONLINE_JUDGE

#ifdef ONLINE_JUDGE
#define getchar() ((p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin)),p1==p2)?EOF:*p1++)
#endif

namespace fast_IO {
	char buf[1<<21],*p1=buf,*p2=buf;

	template <typename T> inline void read(T &x) {
		x=0;T f=1;char ch=getchar();
		for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
		for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
	}
	template <typename T,typename... Args> inline void read(T& x,Args& ...args) {
		read(x),read(args...);
	}

	char buf2[1<<21],a[80];int p,p3=-1;

	inline void flush() {fwrite(buf2,1,p3+1,stdout),p3=-1;}
	template <typename T> inline void write(T x) {
		if(p3>(1<<20)) flush();
		if(x<0) buf2[++p3]='-',x=-x;
		do {a[++p]=x%10+48;} while(x/=10);
		do {buf2[++p3]=a[p];} while(--p);
		buf2[++p3]='
';flush();
	}
	template <typename T,typename... Args> inline void write(T x,Args ...args) {
		write(x),write(args...);
	}
}

using fast_IO :: read;
using fast_IO :: write;
using fast_IO :: flush;

const int maxn = 3e6+10;
const int mod = 1e9+7;
const int N = 1e5+10;

int phi[maxn],pri[maxn],vis[maxn],tot,w[maxn],varphi[maxn];

void sieve() {
	phi[1]=1;w[1]=1;
	for(int i=2;i<maxn;i++) {
		if(!vis[i]) pri[++tot]=i,phi[i]=i-1,w[i]=i;
		for(int j=1;j<=tot&&i*pri[j]<maxn;j++) {
			vis[i*pri[j]]=1;
			if(i%pri[j]==0) {phi[i*pri[j]]=phi[i]*pri[j],w[i*pri[j]]=w[i];break;}
			phi[i*pri[j]]=phi[i]*(pri[j]-1);
			w[i*pri[j]]=w[i]*pri[j];
		}
	}
	for(int i=1;i<maxn;i++) varphi[i]=phi[i],phi[i]=(phi[i]+phi[i-1])%mod;
}

map<int,int> Phi,S[N];

int sum_phi(int x) {
	if(x<maxn) return phi[x];
	if(Phi[x]) return Phi[x];
	int res=1ll*x*(x+1)/2%mod,i=2;
	while(i<=x) {
		int pre=i;i=x/(x/i);
		res=(res-1ll*(i-pre+1)*sum_phi(x/i)%mod)%mod;i++;
	}return Phi[x]=(res+mod)%mod;
}

vector <int > Div;

int s(int n,int m) {
	if(!n) return 0;
	if(!m) return 0;
	if(n==1) return sum_phi(m);
	if(S[n][m]) return S[n][m];
	int y=n/w[n],res=0;
	for(int i=1;i*i<=w[n];i++) {
		if(w[n]%i) continue;
		int j=w[n]/i;
		res=(res+1ll*varphi[j]*s(i,m/i)%mod)%mod;
		if(i!=j) res=(res+1ll*varphi[i]*s(j,m/j)%mod)%mod;
	}
	return S[n][m]=1ll*res*y%mod;
}

int main() {
	sieve();
	int n,m;read(n,m);int ans=0;
	for(int i=1;i<=n;i++) ans=(ans+s(i,m))%mod;
	write(ans);
	flush();
	return 0;
}