【BZOJ3512】DZY Loves Math IV(杜教筛)

点此看题面

大致题意:(sum_{i=1}^nsum_{j=1}^mphi(ij))

推式子

好妙的一道题。。。

考虑(n)(m)的大小相差很多,因此我们可以考虑枚举(n),即令:

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

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

假设(n=prod P_i^{k_i}),则(phi(n)=phi(prod P_i) imes prod P_i^{k_i-1})

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

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

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

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

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

把第一个(phi)和最后一个(phi)合并,也就得到:

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

套路地枚举(d)得到:

[S_n(m)=p imes sum_{d|q}phi(frac qd)sum_{i=1}^{lfloorfrac md floor}phi(di)=p imes sum_{d|q}phi(frac qd)S_d(lfloorfrac md floor) ]

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

代码

#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
",t),0;
}
原文地址:https://www.cnblogs.com/chenxiaoran666/p/BZOJ3512.html