「学习笔记」杜教筛

「学习笔记」杜教筛

算法简介

这是一种通过构造函数 \(g(x)\) 来求一类积性函数前缀和的做法，方法比较精妙

考虑我们要求函数 \(f\) 的前缀和 \(S(n) = \sum_{i=1}^n f(i)\) ，已经有构造好的积性函数 \(g\) 。

将 \(f, g\) 做狄利克雷卷积，此时推式子可以得到：

\[\sum_{k=1}^n(f*g)(k) = \sum_{k=1}^n\sum_{d|k} f(d)g(\frac{k}{d}) \\ = \sum_{d=1}^n g(d)\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor} f(k) \\ = \sum_{d=1}^n g(d)S(\lfloor\frac{n}{d}\rfloor) \]

进一步观察发现：

\[S(n) =\sum_{k=1}^n(f*g)(k) - \sum_{d=2}^n g(d)S(\lfloor\frac{n}{d}\rfloor) \]

对于前面部分要求快速\(O(1)\)得到 \((f*g)\) 的前缀和。

对于后面部分需要计算的 \(S(\lfloor\frac{n}{d}\rfloor)\) 不超过 \(2\sqrt{n}\) 个，套用除数函数的那个证明即可。

先套用主定理再积分一下会发现这样做的复杂度是 \(O(n^{\frac{3}{4}})\) 的，如果能线性预处理出前 \(n^{\frac{2}{3}}\) 项的 \(S\) ，那么复杂度就是 \(O(n^{\frac{2}{3}})\) 。

\(g\) 的常用构造

对于欧拉函数 \(\phi\) 求前缀和 \(S(n) = \sum_{i=1}^n \phi(i)\) ，以及莫比乌斯函数 \(\mu\) 求前缀和，\(S(n) = \sum_{i=1}^n \mu(i)\) 均可以用恒等函数 \(I(n)=1\) 来做 \(g\)。

因为可以证明得到以下等式：

\[\phi * 1 = Id \\ \mu * 1 = e=[n=1] \]

这两者的前缀和都可以 \(O(1)\) 计算。

code

/*program by mangoyang*/
#include<bits/stdc++.h>
#define inf (0x7f7f7f7f)
#define Max(a, b) ((a) > (b) ? (a) : (b))
#define Min(a, b) ((a) < (b) ? (a) : (b))
typedef long long ll;
using namespace std;
template <class T>
inline void read(T &x){
    int f = 0, ch = 0; x = 0;
    for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = 1;
    for(; isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
    if(f) x = -x;
}
const int Lim = 10000005;
map<int, int> ansmiu;
map<int, ll> ansphi;
ll sphi[Lim];
int prime[Lim], smiu[Lim], tot, m;

inline int getmiu(int n){
	if(n < Lim) return smiu[n];
	if(ansmiu[n]) return ansmiu[n];
	int ans = 1;
	for(int l = 2, r; l <= n; l = r + 1){
		r = n / (n / l);
		ans -= (r - l + 1) * getmiu(n / l); 
	}
	return ansmiu[n] = ans;
}
inline ll getphi(int n){
	if(n < Lim) return sphi[n];
	if(ansphi[n]) return ansphi[n];
	ll ans = 1ll * n * (n + 1) / 2;
	for(int l = 2, r; l <= n; l = r + 1){
		r = n / (n / l);
		ans -= 1ll * (r - l + 1) * getphi(n / l);
	}
	return ansphi[n] = ans;
}
signed main(){
	smiu[1] = sphi[1] = 1;
	for(int i = 2; i < Lim; i++){
		if(!sphi[i]) sphi[i] = i - 1, smiu[i] = -1, prime[++tot] = i;
		for(int j = 1; j <= tot && i * prime[j] < Lim; j++){
			if(i % prime[j] == 0){
				sphi[i*prime[j]] = sphi[i] * prime[j];
				smiu[i*prime[j]] = 0;
				break;
			}
			sphi[i*prime[j]] = sphi[i] * (prime[j] - 1);
			smiu[i*prime[j]] = -smiu[i];
		}
	}
	for(int i = 2; i < Lim; i++) 
		sphi[i] += sphi[i-1], smiu[i] += smiu[i-1];
	read(m);
	for(int i = 1, x; i <= m; i++)
		read(x), cout << getphi(x) << " " << getmiu(x) << endl;
	return 0;
}

参考资料：cly-none与12.15日的讲课