杜教筛学习总结

杜教筛学习总结

前言

一直听说杜教筛非常nb，但关于它的学习一直被鸽= =但最近遇到太多数学题辣，所以不得不把这个坑填上了。

推荐博客：
传送门1
传送门2
传送门3

另外杜教筛可能需要一些前置知识，之前写过一篇关于莫比乌斯函数的，就顺便贴上来吧：传送门

正文

数论函数：我们平时遇到的一些特殊函数比如(varphi,mu)这种都属于数论函数。

积性函数：
定义：如果一个已知函数为数论函数，且(f(1)=1)，并且满足对于任意两个互质的(p,q)，有：(f(pcdot q)=f(p)cdot f(q))，那么则称(f)为积性函数。特殊地，若对于任意两个数(p,q)，都有(f(pcdot q)=f(p)cdot f(q))，那么就称这类函数为完全积性函数。

常见积性函数：

• (mu(n))——莫比乌斯函数，这是十分常见，需要对其性质有一定了解；

• (varphi(n))——欧拉函数，定义为(1)~(n)中与(n)互质的数的个数；

• (d(n))——约数个数，表示约数的个数。

• (完全积性函数):

• (varepsilon(n))——元函数，当(n=1)时其值为1，其余时候其值为0；

• (I(n))——恒等函数，不论(n)为何值，其值恒等于1；

• (id(n))——单位函数，(id(n)=n)；

这几个完全积性函数看似较为简单，其实有着很大的用处。它们完全积性的性质也是很好证明的。

狄利克雷卷积：
定义：两个函数(f(n),g(n))的狄利克雷卷积记为((f*g)(n))，表达式等价于(sum_{d|n}f(d)g(frac{n}{d}))。前面的括号代表卷积对象，后面的括号代表卷积的范围。
性质：
狄利克雷卷积运算具有一些很有用的性质：

• 交换率：(f*g=g*f);
• 结合律：((f*g)*h=f*(g*h));
• 分配律：(f*(g+h)=f*g+f*h).

另外还有比较重要的一个性质：
若(f,g)都为积性函数，那么(f*g)也为积性函数。
在杜教筛中，通过卷积来判断积性函数，是很常用的一个技巧。
另外对于我们上面所说的元函数(varepsilon)，有(varepsilon*f=f)，即元函数类似于单位1在乘法中的作用。

应用：

• 对于莫比乌斯函数(mu)来说，有(mu*I=varepsilon)；
• 对于欧拉函数(varphi)来说，有(varphi*I=id).

以上两个表达式也就是这两个函数最常用的性质之一。通过它们也可以来证明一些东西：
比如对于莫比乌斯反演形式一来说，已知(F=f*I)，我们将其两边同时卷上(mu)，那么等式就为(F*mu=f*varepsilon=f)，即(f=F*mu)。
是不是感觉证明一些就容易许多了？
还有对于(varphi*I=id)，两边同时卷上(mu)，那么就有(varphi=id*mu=sum_{d|n}dmu(frac{n}{d}))。这个式子也是很有用的。
总之，熟练掌握(mu,varphi)以及狄利克雷卷积的性质是很有好处的。

杜教筛

前面都是一些前置知识，接下来就步入正题——杜教筛啦。
杜教筛主要就是用来在压线性时间(O(n^{frac{2}{3}}))内求出积性函数的前缀和。是不是感觉很厉害？好吧，其实杜教筛也并非那么的神奇，往下看就知道了~

现在对于积性函数(f)，要求(sum_{i=1}^nf(i))，我们将其记为(S(n))。
构造(h=f*g)，那么有：

[egin{aligned} sum_{i=1}^nh(i)=&sum_{i=1}^nf*g\ =&sum_{i=1}^nsum_{d|i}f(frac{i}{d})g(d)\ =&sum_dg(d)sum_{d|i}f(frac{i}{d})\ =&sum_dg(d)S(lfloorfrac{n}{d} floor)\ end{aligned} ]

现在把第一项单独提出来，那么就有：

[sum_{i=1}^nh(i)=g(1)cdot S(n)+sum_{d=2}^ng(d)S(lfloorfrac{n}{d} floor) ]

式子等价于：

[g(1)cdot S(n)=sum_{i=1}^n(f*g)(i)-sum_{d=2}^ng(d)S(lfloorfrac{n}{d} floor) ]

如果我们能够找到一个积性函数(g)，使得(sum_{i=1}^n(f*g)(i))这部分的前缀和能够快速求，而后面可以直接整除分块，那就能够起到加速的作用的。

举几个常见例子来说明吧：

• 求(sum_{i=1}^nmu(i))。

我们知道:(mu*I=varepsilon)，而(varepsilon)这个函数足够简单！是能够快速求和的，所以我们令(g=varepsilon)，那么求和式子就能化为：

[S(n)=1-sum_{d=2}^nS(lfloorfrac{n}{d} floor) ]

对于后部分，整除分块处理即可。

• 求(sum_{i=1}^nvarphi(i))。

还是利用(varphi)的性质，我们知道当它卷上(I)时，能够变成一个较简单的函数：(id)。这个是支持快速求和的，那么我们直接取(g=I)即可。

• 求(sum_{i=1}^nivarphi(i))。

这个貌似就没有那么好处理了，假设现在有个(g)，我们将其卷积的形式表示出来：(sum_{d|n}(dvarphi(d))g(frac{n}{d}))。
注意观察这个式子，当(g=id)时，多处来的一个(d)能够被消去，并且剩余一个(n)可以提到外面去。
那么自然就会想到取(g=id)啦，这时((f*g)(i)=isum_{d|n}varphi(d)=i)。
是不是感觉挺奇妙的，式子也变的很简单了。

• 求(sum_{i=1}^ni^2varphi(i))。

跟上面同样的思路，先把卷积形式写出来，有：(sum_{d|n}d^2varphi(d)g(frac{n}{d}))。
同样的，想到令(g=id)，就有：(sum_{i=1}^nf*g=sum_{i=1}^ni^2)，其实就是和上面同样的套路~

反正对于寻找(g)这一步，需要一定的观察力以及耐心，实在不行一个一个试就行了。

具体实现的话可以递归来实现，(getsum(n))能够得到(S(n))，然后递归处理即可。
证明的话可以用主定理来证，直接上杜教筛的话复杂度是(O( n^{frac{3}{4}}))的，但我们可以先预处理出(n^{frac{2}{3}})的前缀和，那样复杂度就是(O(n^{frac{2}{3}}))了。
可以加上一个记忆化的操作，可以更快一点，尤其是对于多组数据来说。
下面附上一个模板题的代码：
洛谷P4213

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 5e6 + 5;
int mu[N], p[N];
ll phi[N];
bool chk[N];
unordered_map <int, ll> mp1, mp2;
void init() {
    mu[1] = phi[1] = 1;
    int cnt = 0, k = N - 1;
    for(int i = 2; i <= k; i++) {
        if(!chk[i]) p[++cnt] = i, mu[i] = -1, phi[i] = i - 1;
        for(int j = 1; j <= cnt && i * p[j] <= k; j++) {
            chk[i * p[j]] = 1;
            if(i % p[j] == 0) {mu[i * p[j]] = 0; phi[i * p[j]] = phi[i] * p[j]; break;}
            mu[i * p[j]] = -mu[i]; phi[i * p[j]] = phi[i] * (p[j] - 1);
        }
    }
    for(int i = 1; i < N; i++) mu[i] += mu[i - 1], phi[i] += phi[i - 1];
}
ll djs_mu(int n) {
    if(n <= 5000000) return mu[n];
    if(mp1[n]) return mp1[n];
    ll ans = 1;
    for(int i = 2, j; i <= n; i = j + 1) {
        j = n / (n / i);
        ans -= (j - i + 1) * djs_mu(n / i);
    }
    return mp1[n] = ans;
}
ll djs_phi(int n) {
    if(n <= 5000000) return phi[n];
    if(mp2[n]) return mp2[n];
    ll ans = 1ll * (n + 1) * n / 2;
    for(int i = 2, j; i <= n; i = j + 1) {
        j = n / (n / i);
        ans -= (j - i + 1) * djs_phi(n / i);
    }
    return mp2[n] = ans;
}
int n, T;
int main() {
    init();
    cin >> T;
    while(T--) {
        cin >> n;
        ll ans1 = djs_mu(n);
        ll ans2 = djs_phi(n);
        cout << ans2 << ' ' << ans1 << '
';
    }
    return 0;
}

再举个例子吧~

HDU6706 huntian oy
题意：
定义(f(n,a,b)=sum_{i=1}^nsum_{j=1}^igcd(i^a-j^a,i^b-j^b)[gcd(i,j)=1]\% 10^9+7)。
先有(T)组询问，每组询问给出(n,a,b)，求(f(n,a,b))。

思路：
先有一个公式如下：
若(ageq b)，且(a,b)互质，则有：

[egin{aligned} &gcd(a^n-b^n,a^m-b^m)\ =&gcd(a^m-b^m,a^{n\%m}-b^{n\%m})=a^{gcd(n,m)}-b^{gcd(n,m)} end{aligned} ]

我也不知道这咋证= =
反正有了这个式子之后就好推了：

[egin{aligned} f(n,a,b)=&sum_{i=1}^n sum_{j=1}^i i-j[gcd(i,j)=1]\ =&sum_{i=1}^n sum_{j=1}^i i[gcd(i,j)=1]-sum_{i=1}^n sum_{j=1}^i j[gcd(i,j)=1]\ =&sum_{i=1}^nivarphi(i)-sum_{i=1}^nfrac{ivarphi(i)+[i=1]}{2} end{aligned} ]

至于为什么后面那部分可以化成这样，这涉及到(gcd(n,m)=gcd(n-m,n))，也就是说，与之互素的数是成对出现的，并且其和为(n)。当然(n=1)的情况除外。
之后式子就可以化为这样：

[f(n,a,b)=frac{1}{2}sum_{i=1}^nivarphi(i)-1 ]

之后直接上杜教筛就是了。

Code

#include <bits/stdc++.h>
#define heyuhhh ok
using namespace std;
typedef long long ll;
const int N = 5e6 + 5, MOD = 1e9 + 7, inv6 = 166666668;
int p[N];
ll phi[N];
bool chk[N];
unordered_map <int, int> mp;
inline void init() {
    phi[1] = 1;
    int cnt = 0, k = N - 1, i;
    for(i = 2; i <= k; ++i) {
        if(!chk[i]) p[++cnt] = i, phi[i] = i - 1;
        for(int j = 1; j <= cnt && i * p[j] <= k; j++) {
            chk[i * p[j]] = 1;
            if(i % p[j] == 0) {phi[i * p[j]] = phi[i] * p[j]; break;}
            phi[i * p[j]] = phi[i] * (p[j] - 1);
        }
    }
    for(i = 1; i < N; ++i) {
        phi[i] = phi[i - 1] + 1ll * i * phi[i];
        if(phi[i] >= MOD) phi[i] %= MOD;
    }
}

int djs_phi(int n) {
    if(n <= 5000000) return phi[n];
    if(mp[n]) return mp[n];
    int ans = 1ll * (n + 1) * n % MOD * (2ll * n + 1) % MOD * inv6 % MOD;
    for(register int i = 2, j; i <= n; i = j + 1) {
        j = n / (n / i);
        ans -= ((ll(j - i + 1) * (j + i)) >> 1) % MOD * (ll)djs_phi(n / i) % MOD;
        ans %= MOD;
    }
    if(ans < 0) ans += MOD;
    return mp[n] = ans;
}
int n, a, b, T;
int main() {
#ifdef heyuhhh
    freopen("input.in", "r", stdin);
#else
    ios::sync_with_stdio(false); cin.tie(0);
#endif
    init();
    cin >> T;
    while(T--) {
        cin >> n >> a >> b;
        int ans = djs_phi(n) - 1;
        ans = (1ll * ans * (MOD + 1) / 2) % MOD;
        cout << ans << '
';
    }
    return 0;
}

总结

昨晚匆匆忙忙把这篇博客写完了=，=，感觉还是有很多东西没有写好，emmm以后再慢慢改吧。