zoukankan      html  css  js  c++  java
  • 浅(kou)谈(hu)杜教筛

    引入

    让我们从一道例题开始
    这道题目需要我们求出\(\varphi(n)\)\(\mu(n)\)这两个函数的前缀和
    显然,这可以通过线性筛\(\mathcal O(n)\)完成,但这道题目的\(n\)的范围是\(\le 2^{31}\),线性筛无法完成,我们需要一种更高效的算法————杜教筛
    test


    test

    前置知识

    1.线性筛
    2.积性函数:如果函数\(f(n)\)满足当\(gcd(n,m)=1\)\(f(nm)=f(n)f(m)\),那么我们称\(f(n)\)为积性函数。
    3.完全积性函数:就是积性函数去掉那个互质的条件。
    4.艾弗森约定:对于布尔型变量\(x\)\([x]\)表示\(x\)为真时取\(1\),为假时取\(0\)
    5.一些常见积性函数

    • 元函数\(\epsilon(n)=[n=1]\),是完全积性的。
    • 幂函数\(id_k(n)=n^k\),是完全积性的。
    • 单位函数\(I(n)=id_0(n)\)
    • 约数幂函数\(\sigma_k(n)=\sum_{d|n}d^k\),是积性的。
    • 欧拉函数\(\varphi(n)=\sum_{d≤n}[gcd(d,n)=1]\),是积性的.
    • 莫比乌斯函数\(μ(n)\)满足\(\sum_{d|n}μ(d)=[n=1]=\epsilon(n)\)是积性的

    6.狄利克雷卷积
    \(f,g\)是两个数论函数,它们的狄利克雷卷积是:\((f \otimes g)(n)=\sum_{d|n} f(d) g(\frac nd)\)
    性质:满足交换律,结合律


    概述

    杜教筛,因由圈内著名的杜老师(杜瑜皓)引进而得名,可以以\(\mathcal O(n^{\frac 23})\)的时间复杂度求解一类积性函数的前缀和
    假设\(f(n)\)是我们需要求出的积性函数,令

    \[S(n)=\sum_{i=1}^{n} f(i) \]

    我们再寻找一个积性函数\(g(n)\),令\(h(n)=(f \otimes g)(n)\),那么

    \[\sum_{i=1}^{n} h(i) =\sum_{i=1}^{n} \sum_{d|i} f(d)g(\frac id) \]

    \[=\sum_{d=1}^{n} \sum_{i=1}^{\lceil \frac nd \rceil} g(d)f(i) \]

    \[=\sum_{d=1}^{n} g(d) \sum_{i=1}^{\lceil \frac nd \rceil} f(i) \]

    \[=\sum_{d=1}^{n} g(d)S(\lceil \frac nd \rceil) \]

    移项后得到

    \[g(1)S(n)=\sum_{i=1}^{n} h(i)-\sum_{d=2}^{n}g(d)S(\lceil \frac nd \rceil) \]

    那么,如果我们能够较快地求出\(g(n)\)\(h(n)\)的前缀和,那么减号右边的柿子可以整除分块后递归处理,我们就能够在低于线性的复杂度下完成对\(S(n)\)的求解
    可以证明直接递归并使用\(unordered\_map\)或哈希表记忆化后该算法的复杂度是\(\mathcal O(n^{\frac 34})\)
    如果在此基础上,先线性筛求出前\(n^{\frac 23}\)\(f\)的前缀和,那么我们就可以进一步优化到\(\mathcal O(n^{\frac 23})\)
    这就是杜教筛的核心思想


    例题

    洛谷P4213【模板】杜教筛(Sum)
    就是引入中的那道题。。。
    sol.
    根据经典卷积柿子\(\mu \otimes I=\epsilon,\varphi \otimes I=id_1\)
    \(\epsilon,I,id_1\)的前缀和都可以\(\mathcal O(1)\)求出
    于是直接上模板就行了

    #include<bits/stdc++.h>
    using namespace std;
    const int N=5000000;
    typedef long long ll;
    typedef unsigned int uint;
    int T;
    uint n;
    ll phi[N+10],mu[N+10];
    unordered_map<uint,ll> ans1,ans2;
    inline ll sumphi(uint m){ 
    	if(m<=N) return phi[m];
    	if(ans1.count(m)) return ans1[m];
    	ll ret=1ll*m*(m+1)/2;uint r;
    	for(uint l=2;l<=m;l=r+1){
    		r=m/(m/l);
    		ret-=1ll*(r-l+1)*sumphi(m/l);
    	}
    	return ans1[m]=ret;
    }
    inline ll summu(uint m){ 
    	if(m<=N) return mu[m];
    	if(ans2.count(m)) return ans2[m];
    	ll ret=1;uint r;
    	for(uint l=2;l<=m;l=r+1){
    		r=m/(m/l);
    		ret-=1ll*(r-l+1)*summu(m/l);
    	}
    	return ans2[m]=ret;
    }
    int p[N+10],f[N+10],pcnt;
    inline void init(){
    	mu[1]=phi[1]=1;
    	for(int i=2;i<=N;++i){
    		if(!f[i]) p[++pcnt]=i,phi[i]=i-1,mu[i]=-1;
    		for(int j=1;j<=pcnt&&i*p[j]<=N;++j){
    			int d=i*p[j];f[d]=1;
    			if(i%p[j]) mu[d]=-mu[i],phi[d]=phi[i]*phi[p[j]];
    			else{
    				mu[d]=0;phi[d]=phi[i]*p[j];
    				break;
    			}
    		}
    	}
    	for(int i=2;i<=N;++i) phi[i]+=phi[i-1],mu[i]+=mu[i-1];
    } 
    int main(){
    	scanf("%d",&T);
    	init();
    	while(T--){
    		scanf("%u",&n);
    		printf("%lld %lld\n",sumphi(n),summu(n));
    	}
    	return 0;
    }
    

    洛谷P3768 简单的数学题
    题目要求

    \[\sum_{i=1}^{n} \sum_{j=1}^{n}ij gcd(i,j) \]

    sol.
    因为

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

    \[原式=\sum_{i=1}^{n} \sum_{j=1}^{n} ij \sum_{d|n} \varphi(d) \]

    \[=\sum_{d=1}^{n} \varphi(d) \sum_{d|i} i\sum_{d|j} j \]

    \[=\sum_{d=1}^{n} \varphi(d)d^2 \sum_{i=1}^{\lceil \frac nd \rceil} i \sum_{j=1}^{\lceil \frac nd \rceil} j \]

    \[=\sum_{d=1}^{n} \varphi(d)d^2 (\sum_{i=1}^{\lceil \frac nd \rceil} i)^2 \]

    其中\((\sum_{i=1}^{\lceil \frac nd \rceil} i)^2\)可以数论分块后\(\mathcal O(\sqrt{n})\)求出
    于是我们现在就只需要求\(f(x)=\varphi(x)x^2\)的前缀和了
    因为

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

    又因为\(d*\frac nd=n\),于是

    \[\sum_{d|n}\varphi(d)d^2 (\frac nd)^2 =n^3 \]

    于是令\(g=n^2,h=n^3\)即可
    显然\(g\)的前缀和\(=\frac {n(n+1)(2*n+1)}{6}\),\(h\)的前缀和\(=(\frac{n(n+1)}{2})^2\),可以\(\mathcal O(1)\)求出
    于是上杜教筛即可

    #include<bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    const int N=5e6;
    int phi[N+10],f[N+10],p[N+10],pcnt,mod;
    ll n;
    inline int add(int x,int y){return (x+y>=mod)?x+y-mod:x+y;}
    inline int dec(int x,int y){return (x-y<0)?x-y+mod:x-y;}
    inline void init(){
    	phi[1]=1;
    	for(int i=2;i<=N;++i){
    		if(!f[i]) p[++pcnt]=i,phi[i]=i-1;
    		for(int j=1;j<=pcnt&&i*p[j]<=N;++j){
    			f[i*p[j]]=1;
    			if(i%p[j]) phi[i*p[j]]=1ll*phi[i]*phi[p[j]]%mod;
    			else{
    				phi[i*p[j]]=1ll*phi[i]*p[j]%mod;
    				break;
    			}
    		}
    	}
    	for(int i=1;i<=N;++i) f[i]=add(1ll*phi[i]*i%mod*i%mod,f[i-1]);
    }
    unordered_map<ll,int> ans;
    inline int ksm(int x,int y){
    	int ret=1;
    	for(;y;x=1ll*x*x%mod,y>>=1) if(y&1) ret=1ll*ret*x%mod;
    	return ret;
    }
    int iv2,iv6;
    inline int sum(ll x){x%=mod;return 1ll*x*(x+1)%mod*iv2%mod;}
    inline int sumpow(ll x){x%=mod;return 1ll*x*(x+1)%mod*(2*x+1)%mod*iv6%mod;}
    inline int sumf(ll x){
    	if(x<=N) return f[x];
    	if(ans.count(x)) return ans[x];
    	ll ret=1ll*sum(x)*sum(x)%mod,r;
    	for(ll l=2;l<=x;l=r+1){
    		r=x/(x/l);
    		ret=dec(ret,1ll*dec(sumpow(r),sumpow(l-1))*sumf(x/l)%mod);
    	}
    	ans[x]=ret;
    	return ret;
    }
    int main(){
    //	freopen("lx.out","w",stdout);
    	scanf("%d%lld",&mod,&n);
    	init();
    	int ret=0;iv2=ksm(2,mod-2);iv6=ksm(6,mod-2);ll r;
    	for(ll l=1;l<=n;l=r+1){
    		r=n/(n/l);
    		ret=add(ret,1ll*dec(sumf(r),sumf(l-1))*sum(n/l)%mod*sum(n/l)%mod);
    	}
    	printf("%d\n",ret);
    	return 0;
    }
    

    想必你也注意到了,我们的杜教筛往往筛的是一些与\(\varphi,\mu\)等积性函数相关的函数,于是我们这里给出一些常见的卷积,其他的柿子大部分都能由这些导出

    \[\mu \otimes I=\epsilon \]

    \[\phi \otimes I=id_1 \]

    \[id_k \otimes I=\sigma_k \]

    \[f \otimes \epsilon=f \]

  • 相关阅读:
    51nod 1138 【数学-等差数列】
    hdoj3665【简单DFS】
    hdoj3664【DP】
    51nod1270 【dp】
    51nod 1069【思维】
    关于一些数学符号和概率的阐述;
    51nod 1428【贪心】
    51nod 1133【贪心】
    51nod1127【尺取】
    51nod1126【矩阵快速幂】
  • 原文地址:https://www.cnblogs.com/tqxboomzero/p/14039589.html
Copyright © 2011-2022 走看看