【BZOJ3512】DZY Loves Math IV（杜教筛）

点此看题面

大致题意： 求\(\sum_{i=1}^n\sum_{j=1}^m\phi(ij)\)。

推式子

好妙的一道题。。。

考虑\(n\)和\(m\)的大小相差很多，因此我们可以考虑枚举\(n\)，即令：

\[S_n(m)=\sum_{i=1}^m\phi(ni) \]

这样有什么好处呢？最明显的一个就是，我们可以把原本无比棘手的含有两个变量的\(\phi\)给拆开。

假设\(n=\prod P_i^{k_i}\)，则\(\phi(n)=\phi(\prod P_i)\times \prod P_i^{k_i-1}\)。

不妨令\(p=\prod P_i^{k_i-1},q=\prod P_i\)，则：

\[S_n(m)=p\times \sum_{i=1}^m\phi(qi)=p\times \sum_{i=1}^m\phi(\frac q{gcd(q,i)})\phi(i)\times gcd(q,i) \]

此时式子中虽然有一个\(gcd\)，但并不能无脑上反演（说不定可以，反正我不会）。

看到式子中有这么多\(\phi\)，且\(\phi(\frac q{gcd(q,i)})\)的分母中也刚好有\(gcd(q,i)\)，因此我们可以利用\(\phi*I=Id\)，得到：

\[S_n(m)=p\times \sum_{i=1}^m\phi(\frac q{gcd(q,i)})\phi(i)\times\sum_{d|gcd(q,i)}\phi(d) \]

把第一个\(\phi\)和最后一个\(\phi\)合并，也就得到：

\[S_n(m)=p\times\sum_{i=1}^m\phi(i)\sum_{d|gcd(q,i)}\phi(\frac qd) \]

套路地枚举\(d\)得到：

\[S_n(m)=p\times \sum_{d|q}\phi(\frac qd)\sum_{i=1}^{\lfloor\frac md\rfloor}\phi(di)=p\times \sum_{d|q}\phi(\frac qd)S_d(\lfloor\frac md\rfloor) \]

这个式子可以递归求解，于是就可以上杜教筛了。

代码

#pragma GCC optimize(2)
#pragma GCC optimize("inline")
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 100000
#define X 1000000007
#define Inc(x,y) ((x+=(y))>=X&&(x-=X))
using namespace std;
int n,m,sm;
class DuSieve//杜教筛
{
	private:
		#define SZ 100000
		#define V(x) ((x)<=sm?f1[x]:f2[m/(x)])
		int Pt,P[SZ+5],phi[SZ+5],sphi[N+5],f1[SZ+5],f2[SZ+5];map<int,int> s[SZ+5];
		I int SPhi(CI x)//求Phi
		{
			if(x<=SZ) return sphi[x];if(V(x)) return V(x);RI l,r,t=1LL*x*(x+1)/2%X;
			for(l=2;l<=x;l=r+1) r=x/(x/l),t=(t-1LL*(r-l+1)*SPhi(x/l)%X+X)%X;return V(x)=t;
		}
	public:
		I DuSieve()//预处理
		{
			phi[1]=1;for(RI i=2,j,t;i<=SZ;++i)
				for(!P[i]&&(phi[P[++Pt]=i]=i-1),j=1;j<=Pt&&(t=i*P[j])<=SZ;++j)
					if(P[t]=1,i%P[j]) phi[t]=phi[i]*(P[j]-1);else {phi[t]=phi[i]*P[j];break;}
			for(RI i=1;i<=SZ;++i) sphi[i]=(sphi[i-1]+phi[i])%X;
		}
		I int S(RI n,CI m)//求S
		{
			if(!m) return 0;if(n==1) return SPhi(m);if(m==1) return phi[n];//边界
			if(s[n][m]) return s[n][m];RI i,k=n,p=1,q=1,t=0;//记忆化
			for(i=1;P[i]*P[i]<=k;++i) if(!(k%P[i])) for(k/=P[i],q*=P[i];!(k%P[i]);k/=P[i],p*=P[i]);//质因数分解
			for(k^1&&(q*=k),i=1;i*i<=q;++i) !(q%i)&&//暴力枚举q的因数
				(t=(1LL*phi[q/i]*S(i,m/i)+t)%X,i^(q/i)&&(t=(1LL*phi[i]*S(q/i,m/(q/i))+t)%X));
			return s[n][m]=1LL*p*t%X;//注意最后乘上p
		}
}D;
int main()
{
	RI i,j,t=0;for(scanf("%d%d",&n,&m),sm=sqrt(m),i=1;i<=n;++i) Inc(t,D.S(i,m));//暴力枚举n
	return printf("%d\n",t),0;
}