「学习笔记」杜教筛
算法简介
这是一种通过构造函数 (g(x)) 来求一类积性函数前缀和的做法，方法比较精妙
考虑我们要求函数 (f) 的前缀和 (S(n) = sum_{i=1}^n f(i)) ，已经有构造好的积性函数 (g) 。
将 (f, g) 做狄利克雷卷积，此时推式子可以得到：
[sum_{k=1}^n(f*g)(k) = sum_{k=1}^nsum_{d|k} f(d)g(frac{k}{d}) \
= sum_{d=1}^n g(d)sum_{k=1}^{lfloorfrac{n}{d}
floor} f(k) \
= sum_{d=1}^n g(d)S(lfloorfrac{n}{d}
floor)
]
进一步观察发现：
[S(n) =sum_{k=1}^n(f*g)(k) - sum_{d=2}^n g(d)S(lfloorfrac{n}{d}
floor)
]
对于前面部分要求快速(O(1))得到 ((f*g)) 的前缀和。
对于后面部分需要计算的 (S(lfloorfrac{n}{d} floor)) 不超过 (2sqrt{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日的讲课