杜教筛算法学习&bzoj3944题解

转载注明出处：http://xxccxcxc.top/post/15/

预备知识：积性函数基础总结（只用到前半部分）

bzoj3944 Sum

题目大意

给定$n<2^{31}$，求$ans1=\sum_{i=1}^{n}\varphi(i)$和$ans2=\sum_{i=1}^{n}\mu(i)$，多组询问。

分析

对于积性函数的前缀和，常规做法是线性筛预处理，查询是$O(1)$，不过需要的空间是$O(n)$的，像这题$n$一大空间就会爆。对此，需要引入新的算法。

杜教筛

这题给的两个积性函数，有两个经典的性质：

$$\sum_{d\mid n}\varphi(d) = n，即\varphi*1=id$$

$$\sum_{d\mid n}\mu(d) = [n=1]，即\mu*1=\varepsilon$$

杜教筛就是利用这类式子，来快速求出目标函数的前缀和的。

先来看算法原理：

设所求函数$S(n)=\sum_{i=1}^{n}f(i)$，构造一个函数$g$：

$$\begin{align} \sum_{i=1}^{n}(f*g)(i)&=\sum_{i=1}^{n}\sum_{d\mid i}f(d)g(\frac{n}{d})\\\\ &=\sum_{d=1}^{n}g(d)\sum_{1\le i\le n,d\mid i}f(\frac{i}{d})\\\\ &=\sum_{d=1}^{n}g(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)\\\\ &=\sum_{d=1}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)\\\\ &=g(1)S(n)+\sum_{d=2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor) \end{align}$$

移项得：

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

这样的话，如果我们找了一个合适的函数$g$，使得$g$和$f*g$的前缀和都可以在$O(1)$的时间内求出，那么问题就转化成了求$S(\lfloor\frac{n}{d}\rfloor)(2\le d\le n)$，大大减小了规模。可以证明，先预处理出$n^{\frac{2}{3}}$内的所有$S$值，再进行记忆化搜索，复杂度为$O(n^{\frac{2}{3}})$。

回到例题，将$f=\mu,g=1$代入，此时$f*g=\varepsilon$，可以得到：

$$S(n)=1-\sum_{d=2}^{n}S(\lfloor\frac{n}{d}\rfloor)$$

同理，将$f=\varphi,g=1$代入，此时$f*g=id$，可以得到：

$$n\cdot S(n)=\frac{n(n+1)}{2}-\sum_{d=2}^{n}S(\lfloor\frac{n}{d}\rfloor)$$

问题得解。

代码

#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <map>
using namespace std;
typedef long long LL;
const int M = 3000000;

LL phi[M+5], mu[M+5];
int prime[M], pcnt;
bool notPri[M+5];
map<int, LL> muAns, phiAns;

void Sieve() {
    phi[1] = mu[1] = 1;
    for (int i = 2; i < M; i++) {
        if (!notPri[i]) {
            prime[pcnt++] = i;
            phi[i] = i - 1;
            mu[i] = -1;
        }
        for (int j = 0; j < pcnt && 1LL * i * prime[j] < M; j++) {
            int t = 1LL * i * prime[j];
            notPri[t] = true;
            if (i % prime[j]) {
                phi[t] = phi[i] * phi[prime[j]];
                mu[t] = mu[i] * mu[prime[j]];
            }
            else {
                phi[t] = phi[i] * prime[j];
                mu[t] = 0;
                break;
            }
        }
    }
    for (int i = 2; i < M; i++) {
        phi[i] += phi[i - 1];
        mu[i] += mu[i - 1];
    }
}

LL CalcMu(LL n) {
    if (n < M) return mu[n];
    if (muAns.count(n)) return muAns[n];
    LL last; LL ans = 1;
    for (LL i = 2; i <= n; i = last + 1) {
        last = n / (n / i);
        ans -= (last - i + 1) * CalcMu(n / i);
    }
    return muAns[n] = ans;
}

LL CalcPhi(LL n) {
    if (n < M) return phi[n];
    if (phiAns.count(n)) return phiAns[n];
    LL last; LL ans = 1LL * n * (n + 1) >> 1;
    for (LL i = 2; i <= n; i = last + 1) {
        last = n / (n / i);
        ans -= (last - i + 1) * CalcPhi(n / i);
    }
    return phiAns[n] = ans;
}

int main() {
    Sieve();
    int T; scanf("%d", &T);
    for (int n; T--; ) {
        scanf("%d", &n);
        printf("%lld %lld\n", CalcPhi(n), CalcMu(n));
    }
    return 0;
}

后记

不想单开一个分类，就放在莫比乌斯反演里了。

关于杜教筛的其他应用和套路，以后会再提到。