• 【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;
    }
    
  • 相关阅读:
    Linux数据链路层的包解析
    Nmap的活跃主机探测常见方法
    甲方安全建设推进思路
    重新学习python类
    python装饰器
    记录一次奇葩渗透中的点点滴滴
    安全情报总结
    机器学习基础
    tensorflow学习笔记(四十五):sess.run(tf.global_variables_initializer()) 做了什么?
    tensorflow学习笔记(二十五):ConfigProto&GPU
  • 原文地址:https://www.cnblogs.com/chenxiaoran666/p/BZOJ3512.html
Copyright © 2020-2023  润新知