zoukankan      html  css  js  c++  java
  • 杜教筛学习笔记

    杜教筛学习笔记

    接着 Dirichlet 卷积,继续学习杜教筛。

    本文中一些函数的定义见此博文

    应用

    通过杜教筛,我们可以快速求出某一数论函数 (f) 的前缀和,即,我们可以在低于线性的时间复杂度内求出 (S(n)=sumlimits_{i=1}^n f(i))

    方法

    杜教筛主要运用一个公式,通过这个公式我们建立了 (S(n))(S(frac{n}{i})) 的关系式,已知两个数论函数 (f,g)(S(n)=sumlimits_{i=1}^n f(i)),则有:

    [sumlimits_{i=1}^n (fast g)(i)=sumlimits_{i=1}^n g(i)cdot S(lfloorfrac{n}{i} floor) ]

    我们给个证明:

    [egin{aligned} sumlimits_{i=1}^n (fast g)(i) \=sumlimits_{i=1}^nsumlimits_{dmid i} f(d)cdot g(frac{i}{d})\=sumlimits_{d=1}^nsumlimits_{i=1}^{lfloor frac{n}{d} floor} g(d)cdot f(i)\=sumlimits_{d=1}^n g(d)cdotsumlimits_{i=1}^{lfloor frac{n}{d} floor} f(i)\=sumlimits_{d=1}^n g(d)cdot S(lfloor frac{n}{d} floor)\=sumlimits_{i=1}^n g(i)cdot S(lfloor frac{n}{i} floor)\&square end{aligned}]

    证明的关键就在于中间对于 (d,i) 枚举顺序的变化,这个证明值得细品。

    那么我们有了这个公式有什么用呢?我们再来导一下:

    [egin{aligned} sumlimits_{i=1}^n (fast g)(i)=sumlimits_{i=1}^n g(i)cdot S(lfloorfrac{n}{i} floor)\sumlimits_{i=1}^n (fast g)(i)=g(1)cdot S(n)+sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor)\g(1)cdot S(n)=sumlimits_{i=1}^n (fast g)(i)-sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor) end{aligned}]

    我们可以轻松地求出 (S(n)) 的值。而 (S(lfloorfrac{n}{i} floor)) 可以通过数论分块来解决。由此看来,杜教筛的时间复杂度为低于线性复杂度。

    例子

    看了理论,可能还有点晕,我们来看几个例子。

    求莫比乌斯函数前缀和

    (S(n)=sumlimits_{i=1}^n mu(i) (n<2^{31}))

    我们根据公式来求,我们令 (f=mu,g=1),由 Dirichlet 卷积知识,我们知道 (muast 1=varepsilon)。那么:

    [egin{aligned} g(1)cdot S(n)=sumlimits_{i=1}^n (fast g)(i)-sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor)\1cdot S(n)=sumlimits_{i=1}^n varepsilon(i)-sumlimits_{i=2}^n 1cdot S(lfloorfrac{n}{i} floor)\S(n)=1-sumlimits_{i=2}^n S(lfloorfrac{n}{i} floor) end{aligned}]

    根据数论分块的知识,(lfloorfrac{n}{i} floor) 的取值有 (O(sqrt{n})) 个。那么时间复杂度为 (O(sumlimits_{i=1}2^{frac{1}{2^i}}) approx O(n^{frac{3}{4}}))

    这明显会超时,我们需要进一步优化。其实可以这样想,一部分前缀和我们先预处理出来,另一部分再通过杜教筛来算。我们设预处理前 (k) 个,后 (n-k) 个用杜教筛做,那么此时时间复杂度是 (O(k+(n-k)^{frac{3}{4}})),由基本不等式得,(k=(n-k)^{frac{3}{4}}) 原式取最小值。我们算一下发现 (k=n^{frac{2}{3}}) 是不错的选择,此时复杂度大约为 (O(n^{frac{2}{3}}))

    //Don't act like a loser.
    //You can only use the code for studying or finding mistakes
    //Or,you'll be punished by Sakyamuni!!!
    #include<bits/stdc++.h>
    #define int long long
    using namespace std;
    
    int read() {
    	char ch=getchar();
    	int f=1,x=0;
    	while(ch<'0'||ch>'9') {
    		if(ch=='-')
    			f=-1;
    		ch=getchar();
    	}
    	while(ch>='0'&&ch<='9') {
    		x=x*10+ch-'0';
    		ch=getchar();
    	}
    	return f*x;
    }
    
    const int maxn=4e7+10;
    
    int n,mu[maxn],p[maxn/10],cnt,s[maxn];
    bool is[maxn];
    map<long,long> mp;//记忆化搜索
    
    void sieve_mu(int n) {
    	fill(is,is+n+1,1);
    	mu[1]=1;
    	
    	for(int i=2;i<=n;i++) {
    		if(is[i]) {
    			p[++cnt]=i;
    			mu[i]=-1;
    		}
    		for(int j=1;j<=cnt;j++) {
    			if(p[j]*i>n) {
    				break;
    			}
    			
    			is[p[j]*i]=0;
    			if(i%p[j]==0) {
    				mu[i*p[j]]=0;
    				break;
    			}
    			mu[i*p[j]]=-mu[i];
    		}
    	}
    	for(int i=1;i<=n;i++){
    		s[i]=s[i-1]+mu[i];
    	}
    }
    
    int pre(int n) {
    	if(n<=4e7) {
    		return s[n];
    	}
    	if(mp[n]) {
    		return mp[n];
    	}
    	int j=2,ret=1;
    	for(int i=2;i<=n;i=j+1) {
    		j=(n/(n/i));//数论分块
    		
    		ret-=(j-i+1)*pre(n/i);
    	}
    	return mp[n]=ret;
    }
    
    signed main() {
    	sieve_mu(4e7);
    	int t=read();
    	
    	while(t--) {
    		printf("%lld
    ",pre(read()));
    	}
    	return 0;
    }
    

    求欧拉函数前缀和

    (S(n)=sumlimits_{i=1}^n varphi(i) (n<2^{31}))

    我们根据公式来求,我们令 (f=varphi,g=1),由 Dirichlet 卷积知识,我们知道 (varphiast 1=id)。那么:

    [egin{aligned} g(1)cdot S(n)=sumlimits_{i=1}^n (fast g)(i)-sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor)\1cdot S(n)=sumlimits_{i=1}^n id(i)-sumlimits_{i=2}^n 1cdot S(lfloorfrac{n}{i} floor)\S(n)=frac{n^2+n}{2}-sumlimits_{i=2}^n S(lfloorfrac{n}{i} floor) end{aligned}]

    和之前一样,我们还是采取一部分直接算,一部分杜教筛的想法,时间复杂度为 (O(n^{frac{2}{3}}))

    //Don't act like a loser.
    //You can only use the code for studying or finding mistakes
    //Or,you'll be punished by Sakyamuni!!!>
    #define int long long
    using namespace std;
    
    int read() {
    	char ch=getchar();
    	int f=1,x=0;
    	while(ch<'0'||ch>'9') {
    		if(ch=='-')
    			f=-1;
    		ch=getchar();
    	}
    	while(ch>='0'&&ch<='9') {
    		x=x*10+ch-'0';
    		ch=getchar();
    	}
    	return f*x;
    }
    
    const int maxn=1e7+10;
    
    int n,phi[maxn],p[maxn/10],cnt,s[maxn];
    bool is[maxn];
    map<long,long> mp;
    
    void sieve_euler(int n) {
    	fill(is,is+n+1,1);
    	phi[1]=1;
    	
    	for(int i=2;i<=n;i++) {
    		if(is[i]) {
    			p[++cnt]=i;
    			phi[i]=i-1;
    		}
    		for(int j=1;j<=cnt;j++) {
    			if(p[j]*i>n) {
    				break;
    			}
    			
    			is[p[j]*i]=0;
    			if(i%p[j]==0) {
    				phi[i*p[j]]=p[j]*phi[i];
    				break;
    			}
    			phi[i*p[j]]=(p[j]-1)*phi[i];
    		}
    	}
    	for(int i=1;i<=n;i++){
    		s[i]=s[i-1]+phi[i];
    	}
    }
    
    int pre(int n) {
    	if(n<=1e7) {
    		return s[n];
    	}
    	if(mp[n]) {
    		return mp[n];
    	}
    	int j=2,ret=(n*(n+1)/2);
    	for(int i=2;i<=n;i=j+1) {
    		j=(n/(n/i));
    		
    		ret-=(j-i+1)*pre(n/i);
    	}
    	return mp[n]=ret;
    }
    
    signed main() {
    	sieve_euler(1e7);
    	int t=read();
    	
    	while(t--) {
    		printf("%lld
    ",pre(read()));
    	}
    	return 0;
    }
    
    

    这两个代码结合起来就可以通过洛谷上的模板题。

    求莫比乌斯函数平方前缀和

    (S(n)=sumlimits_{i=1}^n mu^2(i) (n<2^{31}))

    我们继续根据公式来求,我们令 (f=mu^2,g=1),由 Dirichlet 卷积知识,我们知道 (mu^2*1=mu)。那么:

    [egin{aligned} g(1)cdot S(n)=sumlimits_{i=1}^n (fast g)(i)-sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor)\1cdot S(n)=sumlimits_{i=1}^n mu(i)-sumlimits_{i=2}^n 1cdot S(lfloorfrac{n}{i} floor)\S(n)=sumlimits_{i=1}^n mu(i)-sumlimits_{i=2}^n S(lfloorfrac{n}{i} floor) end{aligned}]

    所以我们在求莫比乌斯函数平方和时,要求出莫比乌斯函数的前缀和,也是杜教筛,总时间复杂度 (O(n^{frac{2}{3}}))

    求莫比乌斯函数值平方前缀和

    (S(n)=sumlimits_{i=1}^n (mu(i))^2 (n<2^{31}))

    这个题和之前的就有所不同了。我们令 (f(i)=mu(i)^2),但是要找一个方便计算的 (g) 才行。

    我们选取函数 (g(x)=[x=k^2 kin N^+])。我们来算一下 (g)(f) 的狄利克雷卷积。我们会发现 (fast g=1)?!简要证明一下:

    ((fast g)(n)=sumlimits_{dmid n} g(d)cdot f(frac{n}{d}))。首先分类讨论,如果 (d) 是一个完全平方数,只有当 (d)(n) 的所有约数中最大的完全平方数时,乘积为 (1),否则 (mu(frac{n}{d})=0)。如果 (d) 不是完全平方数,那么 (g(d)=0)。得证。

    那我们的计算就简化了不少。

    [egin{aligned} g(1)cdot S(n)=sumlimits_{i=1}^n (fast g)(i)-sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor)\1cdot S(n)=sumlimits_{i=1}^n 1-sumlimits_{i=2}^{sqrt{n}} 1cdot S(lfloorfrac{n}{i^2} floor)\S(n)=n-sumlimits_{i=2}^{sqrt{n}} S(lfloorfrac{n}{i^2} floor) end{aligned}]

    我们甚至不用数论分块就可以解决此题。采用同样的优化方法,时间复杂度还是约为 (O(n^{frac{2}{3}}))

    相关练习题:

    1. 模板(简单版)
    2. 模板(强化版)
    3. 完全平方数
    4. Sum
    5. Lucas的数论
    6. 自己找
  • 相关阅读:
    Protected和Default的区别
    将数组中负数放在正数前面
    java.io包和杯子测楼
    hadoop基础
    极限编程和JUnit
    接口和抽象类
    C# 中窗口AutoScaleMode属性
    计算机的自启动管理
    labview中的移位寄存器、循环隧道,自动索引隧道的区别
    发现C#winform编程中不常用的控件(一)<FlowLayoutPanel控件><拆分器控件Splitcontainer >
  • 原文地址:https://www.cnblogs.com/huayucaiji/p/dujiaoshai.html
Copyright © 2011-2022 走看看