[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\)，式子就变成：

\[\begin{align} s(n,m)&=y\cdot\sum_{i=1}^{m}\varphi(wi)\\ &=y\cdot \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\)化一下，即：

\[\begin{align} s(n,m)&=y\cdot \sum_{i=1}^{m}\varphi(\frac{w}{(i,w)})\cdot\varphi(i)\cdot\sum_{d|i\&d|w}\varphi(d) \\ &=y\cdot \sum_{i=1}^{m}\varphi(i)\cdot\sum_{d|i\&d|w}\varphi( \frac{(i,w)}{d}) \varphi(\frac{w}{(i,w)})\\ &=y\cdot \sum_{i=1}^{m}\varphi(i)\cdot\sum_{d|i\&d|w}\varphi( \frac{w}{d})\\ &=y\cdot \sum_{d|w}\varphi(\frac{w}{d})\sum_{i=1}^{\lfloor\frac{m}{d}\rfloor}\varphi(id)\\ &=y\cdot \sum_{d|w}\varphi(\frac{w}{d})s(d,\lfloor\frac{m}{d}\rfloor) \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]='\n';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;
}