• 【数学】杜教筛


    前置知识:莫比乌斯反演,狄利克雷卷积等亿点数论知识。

    不会不要紧,这里顺便一提。

    积性函数

    对于一个数论函数 (mathbf f),若对于任意的 (aperp b),都有 (mathbf fleft(ab ight) = mathbf fleft(a ight)mathbf fleft(b ight)),则称 (f) 为一个积性函数。即使不满足 (aperp b) 时仍有 (mathbf fleft(ab ight) = mathbf fleft(a ight)mathbf fleft(b ight)) 的数论函数被称为完全积性函数

    一些积性函数的例子:

    • (1),1 函数。对于任意的 (n),有 (1left(n ight) = 1)。(完全积性函数)
    • (Id_k),幂函数。(Id_kleft(n ight) = n^k),当 (k = 1) 时下标可省略不写。(完全记性函数)
    • (varphi),欧拉函数。
    • (mu) 莫比乌斯函数。
    • (varepsilon)。单位元。(varepsilonleft(n ight) = [n = 1])。(完全积性函数)
    • (gcdleft(k, n ight))(k) 为常数)。

    狄利克雷卷积

    对于两个数论函数 (mathbf f)(mathbf g),定义其狄利克雷卷积为:

    [mathbf {f *g}left(n ight) = sum_{d|n}mathbf fleft(d ight)mathbf gleft(frac{n}{d} ight) ]

    可以证明,(*) 运算符合交换律和结合律。

    其单位元为 (varepsilon)

    两个积性函数的狄利克雷卷积仍旧是一个积性函数。

    一些例子

    • (mathbf f * varepsilon = mathbf f)
    • (mu * 1 = varepsilon)
    • (varphi * 1 = Id)
    • (Id * mu = varphi)
    • (1 * Id = sigma)
    • (1 * 1 = d)

    后面两个分别是约数和函数和约数个数函数。

    莫比乌斯反演

    若两个数论函数 (mathbf f)(mathbf g) 满足如下关系

    [mathbf fleft(n ight) = sum_{d|n}mathbf gleft(d ight) ]

    那么就有

    [mathbf gleft(n ight) = sum_{d|n}mathbf fleft(d ight)muleft(frac{n}{d} ight) ]

    证明:把上面的式子写成 (mathbf f = mathbf g * 1),然后两边同时卷上 (mu) 即可。

    杜教筛

    杜教筛可以用于在非线性时间内求数论函数前缀和。

    直接从具体例子入手,可能会好理解一些。

    思考一个问题:(varphi) 函数的前缀和怎么求?求出其 (1)(n) 的前缀和,(n)(10^9) 级别。

    当数据范围很小的时候,我们的做法是通过线性筛筛出所有的 (varphi) 函数值,然后求个前缀和。

    然而在这种数据范围下,这种做法显然行不通了。

    那我们来思考几个简单的问题:

    • (1) 的前缀和是什么?显然是 (n)
    • (Id) 的前缀和是什么?显然是 (frac{nleft(n + 1 ight)}{2})
    • (Idleft(1 ight)) 的值是多少?显然是 (1)

    而我们考虑到 (varphi * 1 = Id)

    考虑能不能用这些很好求前缀和的数论函数来求 (mu) 的前缀和。

    (sum_{i = 1}^nvarphileft(i ight) = sleft(n ight))

    那么就有

    [s * 1 = sum_{d|n}sum_{i = 1}^dvarphileft(i ight) ]

    改变枚举顺序,考虑每个 (1left(d ight) =1) 的贡献,那么原式等于

    [sum_{d = 1}^nsum_{i = 1}^{n / d}varphileft(i ight) ]

    [=sum_{d = 1}^nsleft(lfloorfrac{n}{d} floor ight) ]

    那么 (sleft(n ight)) 值是什么?答案是:

    [sleft(n ight) = sum_{d = 1}^nsleft(lfloorfrac{n}{d} floor ight) - sum_{d = 2}^nsleft(lfloorfrac{n}{d} floor ight) ]

    [= sum_{i = 2}^nleft(varphi * 1 ight)left(i ight) - sum_{d = 1}^nsleft(lfloorfrac{n}{d} floor ight) ]

    [= sum_{i = 2}^nIdleft(i ight) - sum_{d = 1}^nsleft(lfloorfrac{n}{d} floor ight) ]

    [= dfrac{nleft(n + 1 ight)}{2} - sum_{d = 1}^nsleft(lfloorfrac{n}{d} floor ight) ]

    哇,我们发现了什么?

    后面的那个 (s) 我们可以在一边数论分块,一边递归的求。

    这样我们就可以筛出一些比较小的时候的 (s) 的值,然后用个 map 记忆化一下计算。

    LL sum_phi(LL n) {
    	if(n <= N) return sphi0[(int)n];
    	if(sphi.count(n)) return sphi[n];
    	LL res = 1ll * n * (n + 1) / 2;
    	for(LL l = 2, r = 1; l <= n; l = r + 1) {
    		r = n / (n / l);
    		res -= 1ll * (r - l + 1) * sum_phi(n / l);
    	}
    	return sphi[n] = res;
    }
    

    推广到一般情况,则可以得出下面的式子

    对于一个数论函数 (mathbf f),设 (sleft(n ight) = sum_{i = 1}^nmathbf fleft(i ight)),另外有数论函数 (mathbf g)

    [sleft(n ight) = dfrac{sum_{i = 1}^nleft(mathbf {f * g} ight)left(i ight) - sum_{i = 2}^nmathbf gleft(i ight)sleft(lfloorfrac{n}{i} floor ight)}{mathbf gleft(1 ight)} ]

    如果您能够快速求出 (left(mathbf {f * g} ight)) 的前缀和,(mathbf g) 的前缀和,以及 (mathbf gleft(1 ight)) 的值的话,也就意味着您能够用类似于上面求 (sum_{i = 1}^n varphileft(i ight)) 的方式,快速求出 (mathbf f) 的前缀和!

    算法的复杂度是 (mathcal O(n^{frac{2}{3}})) 的。

    这就是杜教筛啦。

    完整代码 洛谷 P4213 【模板】杜教筛

    #include <bits/stdc++.h>
    #define LL long long
    
    template <typename Temp> inline void read(Temp & res) {
    	Temp fh = 1; res = 0; char ch = getchar();
    	for(; !isdigit(ch); ch = getchar()) if(ch == '-') fh = -1;
    	for(; isdigit(ch); ch = getchar()) res = (res << 3) + (res << 1) + (ch ^ '0');
    	res = res * fh;
    }
    
    namespace Math {
    	#define N (10000000)
    	const int Maxn = 1e7 + 5;
    	bool isprime[Maxn]; int prime[Maxn], cnt_prime; LL phi0[Maxn], mu0[Maxn];
    	void sieve() {
    		phi0[1] = mu0[1] = 1;
    		for(int i = 2; i <= N; ++i) {
    			if(!isprime[i]) {
    				prime[++cnt_prime] = i;
    				phi0[i] = i - 1; mu0[i] = -1;
    			}
    			for(int j = 1; j <= cnt_prime && prime[j] <= N / i; ++j) {
    				int cur = i * prime[j];
    				isprime[cur] = 1;
    				if(i % prime[j] == 0) {
    					phi0[cur] = phi0[i] * prime[j];
    					mu0[cur] = 0; break;
    				} else {
    					phi0[cur] = phi0[i] * (prime[j] - 1);
    					mu0[cur] = -mu0[i];
    				}
    			}
    		}
    		for(int i = 1; i <= N; ++i) phi0[i] += phi0[i - 1], mu0[i] += mu0[i - 1];
    	}
    	
    	std::map<LL, LL> mu, phi;
    	//mu * 1 = epsilon
    	LL sum_mu(LL n) {
    		if(n <= N) return mu0[(int)n];
    		if(mu.count(n)) return mu[n];
    		LL res = 1ll;
    		for(LL l = 2, r = 1; l <= n; l = r + 1) {
    			r = n / (n / l);
    			res -= 1ll * (r - l + 1) * sum_mu(n / l);
    		}
    		return mu[n] = res;
    	}
    	//phi * 1 = id 
    	LL sum_phi(LL n) {
    		if(n <= N) return phi0[(int)n];
    		if(phi.count(n)) return phi[n];
    		LL res = 1ll * n * (n + 1) / 2;
    		for(LL l = 2, r = 1; l <= n; l = r + 1) {
    			r = n / (n / l);
    			res -= 1ll * (r - l + 1) * sum_phi(n / l);
    		}
    		return phi[n] = res;
    	}
    	#undef N
    }
    
    int T, n;
    
    int main() {
    	using namespace Math;
    	sieve(); 
    	read(T);
    	while(T--) {
    		read(n);
    		printf("%lld %lld
    ", sum_phi(n), sum_mu(n));
    	}
    	return 0;
    }
    
  • 相关阅读:
    【转】c++继承中的内存布局
    Google 开源项目风格指南
    常见面试题
    PHP7.1中使用openssl替换mcrypt
    phpunit实践笔记
    PHP的错误处理
    CI的扩展机制
    #CI的MVC实现
    Laravel中的队列处理
    laravel的模块化是如何实现的
  • 原文地址:https://www.cnblogs.com/zimujun/p/14563547.html
Copyright © 2020-2023  润新知