杜教筛_牛客博客

介绍

我们在 OI 题目中，不难会遇到给定一个积性函数 $f(x)$ ，求

$\sum_{i=1}^N f(x)$

的题目。虽然对于积性函数我们已经有成熟的线性筛方法来说可以做到 $O(n)$ ，~~可是有的出题人就是毒瘤~~，可是有的题目的数据范围要求给出一个复杂度小于线性的做法，这时候神仙就搞出了杜教筛这一东西。

前置芝士

芝士1：线性筛

首先我们知道，基本上所有的积性函数都可以线性筛。
积性函数的特点：

$\forall x,y,(x,y) = 1 \Rightarrow f(xy) = f(x)f(y)$

如果我们想线性筛积性函数，首先这个过程一定是建立在欧拉筛素数上的，然后我们考虑一些性质：

当 $p$ 为质数时， $f(p)$ 的取值
当这个数不是第一次被筛到的时候，取值的特性

我们拿 $\mu$ 举例：
什么你还不知道 $\mu$ 是什么？快去看看莫比乌斯反演吧。
首先发现 $\mu(p) = -1$ ， $p \in prime$ 。
然后发现当 $x$ 能不止被一种乘积筛出来的时候， $\mu(x) = 0$ 。
如果是被最小的质因数筛出的时候，就可以依据积性函数的定义来更新了。
直接放一个线性筛 $\mu$ 和 $\phi$ 的代码吧。

inline void pre(int n=MAXN-5){
    phi[1] = mu[1] = 1;
    FOR(i,2,n){
        if(!p[i]){
            prime[++cnt] = i;
            mu[i] = -1;phi[i] = i-1;
        }
        for(Re int j = 1;j <= cnt && i*prime[j] <= n;++j){
            p[i*prime[j]] = true;
            if(i%prime[j]){
                phi[i*prime[j]] = phi[i]*phi[prime[j]];
                mu[i*prime[j]] = -mu[i];
            }
            else{
                phi[i*prime[j]] = phi[i]*prime[j];
                mu[i*prime[j]] = 0;
                break;
            }
        }
    }
}

芝士2：常见积性函数

纯粹是列举一下防止后面因为符号不同的问题而看不懂。 $\mu(n)$ ，莫比乌斯函数。 $\phi(n)$ ，欧拉函数。 $d(n)$ ，约数个数，即 $d(n) = \sum_{i=1}^n [i|n]$ 。 $\sigma(n)$ ，约数和，即 $\sigma(n) = \sum_{i=1}^n i[i|n]$ 。
完全积性函数： $\forall x,y \Rightarrow f(xy) = f(x)f(y)$
$I(n)$ ：常函数 $I(n) = 1$ 。 $\epsilon(n)$ ：狄利克雷卷积单位元。 $id(n)$ ：单位函数 $id(i) = i$ 。
如果你学过莫比乌斯反演，应该已经掌握了所有的前置芝士了，现在我们进入正题。

杜教筛过程

下文中定义狄利克雷卷积的运算符号为 $*$
我们现在要对于一个积性函数 $f(x)$ ，求

$S(n) = \sum_{i=1}^n f(i)$

考虑找到一个合适的函数 $g$ ，使得 $f*g$ 的前缀和易于计算。
我们考虑 $f*g$ 的前缀和，并加以推导，得
\begin{align}
&\sum_{i=1}^n (fg)(i) \
&= \sum_{i=1}^n \sum_{d|i} g(i)f(\frac{n}{i}) \
&= \sum_{d=1}^n g(d) \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} f(i) \
&= \sum_{i=1}^n g(i) S(\lfloor \frac{n}{i}\rfloor) \
&= g(1)S(n) + \sum_{i=2}^n g(i)S(\lfloor \frac{n}{i}\rfloor)
\end{align*}
考虑进行一些小的移项，最后可得我们一般使用的杜教筛的式子：

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

如果我们对于 $(f*g)$ 的前缀和有快速求法（e.g. $O(1)$ ），那么预处理出 $f$ 函数的前 $n^{\frac{2}{3}}$ ，后面数论分块+递归求解，就可以做到 $O(n^{\frac{2}{3}})$ 的优秀复杂度内求得 $S(n)$ 。
现在我们来尝试用杜教筛来解决求一些函数的前缀和问题。

BZOJ 3944 Sum

题目链接
题目要让你求 $\sum_{i=1}^n \mu(i)$ 和 $\sum_{i=1}^n \phi(i)$ ，其中 $n \leq 2^{31}-1$ 。
首先考虑 $\mu$ ：
首先我们需要找到一个函数，使得它与 $\mu$ 的狄利克雷卷积的前缀和是容易计算的。注意到关于 $\mu$ 的一个性质： $\mu * I = \epsilon$ 。
发现 $\epsilon$ 的前缀和十分好算（就是 $1$ ），并且 $I$ 非常有利于乘法分块，于是我们取 $f(n) = \mu(n),g(n) = I(n)$ ，代入上式，得：

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

类比我们来考虑 $\phi$ ：
发现一个性质： $\phi * I = id$
$id$ 的前缀和可以用等差数列求和公式求得，代入上式可以发现：

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

也就是

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

现在我们来证明一下这东西的复杂度是什么。
设 $T(n)$ 表示计算出 $S(n)$ 的复杂度，可以得到：

$T(n) = O(\sqrt{n}) + \sum_{i=2}^{n}(T(i)+T(\frac{n}{i}))$

展开一层就可以了，因为往下都是高阶小量可以不考虑。

$T(n) =O(\sqrt{n})+\sum_{i=2}^n(O(\sqrt{i})+O(\sqrt{\lfloor \frac{n}{i} \rfloor})) = O(n^{\frac{3}{4}})$

我们可以用筛法处理前 $k$ 个，如果 $k \geq \sqrt{n}$ ，则

$T(n) = \sum_{i=1}^{\frac{n}{k}}O(\sqrt{\lfloor \frac{n}{i} \rfloor}) = O(\frac{n}{\sqrt{k}})$

一般我们 $k$ 取 $n^{\frac{2}{3}}$ ，可以得到复杂度为 $O(n^{\frac{2}{3}})$ 的算法。

#include <tr1/unordered_map>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <climits>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <cmath>
#include <ctime>
#include <queue>
#include <stack>
#include <map>
#include <set>

#define fi first
#define se second
#define U unsigned
#define P std::pair
#define Re register
#define LL long long
#define pb push_back
#define MP std::make_pair
#define all(x) x.begin(),x.end()
#define CLR(i,a) memset(i,a,sizeof(i))
#define FOR(i,a,b) for(Re int i = a;i <= b;++i)
#define ROF(i,a,b) for(Re int i = a;i >= b;--i)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl
using namespace std::tr1;
const int MAXN = 6000000+5;
bool p[MAXN];
int mu[MAXN],phi[MAXN],prime[10001],cnt;
LL sum1[MAXN];int sum2[MAXN];

inline void pre(int n=MAXN-5){
    phi[1] = mu[1] = 1;
    FOR(i,2,n){
        if(!p[i]){
            prime[++cnt] = i;
            mu[i] = -1;phi[i] = i-1;
        }
        for(Re int j = 1;j <= cnt && i*prime[j] <= n;++j){
            p[i*prime[j]] = true;
            if(i%prime[j]){
                phi[i*prime[j]] = phi[i]*phi[prime[j]];
                mu[i*prime[j]] = -mu[i];
            }
            else{
                phi[i*prime[j]] = phi[i]*prime[j];
                mu[i*prime[j]] = 0;
                break;
            }
        }
    }
    FOR(i,1,n) sum1[i] = sum1[i-1]+phi[i],sum2[i] = sum2[i-1]+mu[i];
}

unordered_map<LL,LL> S1;
unordered_map<int,int> S2;

inline LL calc1(LL x){
    if(x <= MAXN-5) return sum1[x];
    if(S1[x]) return S1[x];
    LL ans = x*(x+1)/2;
    for(Re LL l=2,r;l <= x;l = r+1){
        r = x/(x/l);
        ans -= calc1(x/l)*(r-l+1);
    }
    return S1[x] = ans;
}

inline int calc2(int x){
    if(x <= MAXN-5) return sum2[x];
    if(S2[x]) return S2[x];
    int ans = 1;
    for(Re int l = 2,r;l >= 0 && l <= x;l = r+1){
        r = x/(x/l);
        ans -= calc2(x/l)*(r-l+1);
    }
    return S2[x] = ans;
}

inline void Solve(){
    int n;scanf("%d",&n);
    printf("%lld %d\n",calc1(n),calc2(n));
}

int main(){
    pre();
    int T;scanf("%d",&T);
    while(T--) Solve();
    return 0;
}

两道小题目

Sum 加强版

$\sum_{i=1}^n (i*\phi(i))$

求上式的值，其中 $n$ 是线性筛跑不过的量级。
我们一眼看不出来 $g$ 用什么好，不妨先假定有了这个函数并代入狄利克雷卷积的形式： $\sum_{d|n} (d*\phi(d))*g(\frac{n}{d})$ 。
这个常数 d 太烦人了，不如我们想办法消去它。我们尝试取 $g = I$ ，代入化简得到

$\sum_{d|n} (d*\phi(d))*\frac{n}{d} = n\sum_{d|n} \phi(d) = n^2$

这样的话前缀和就很好求了，带回原式可得：

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

大家都会 $O(1)$ 求 $\sum_{i=1}^n i^2$ ，也就是

$\sum_{i=1}^n i^2 = \frac{n(n+1)(2n+1)}{6}$

Hdu 5608 Function

题目链接
告诉你一个数论函数 $f(d)$ ，满足

$\sum_{d|n} f(d) = n^2-3n+2$

多次询问 $f$ 的前缀和。 $T \leq 500,N \leq 10^9$ ，只有五组 $N > 10^6$ 。
长得十分像莫比乌斯反演的经典形式。
我们设 $g(n) = n^2-3n+2$ ，原式转化为

$g(n) = \sum_{d|n} f(d)$

反演一波可以得到：

$f(n) = \sum_{d|n} g(d)\mu(\frac{n}{d})$

我们要求的答案变成了：

$\sum_{i=1}^n \sum_{d|n} g(d) \mu(\frac{n}{d})$

进行一些微小的变换，可以推导出
\begin{align}
& \sum_{i=1}^n \sum_{d|n} g(d) \mu(\frac{n}{d})\
&= \sum_{d=1}^n g(d) \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \mu(i)
\end{align}
发现 $g$ 的前缀和十分好求，这样就可以通过数论分块+杜教筛 $\mu$ 来做了。

#include <unordered_map>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <climits>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <cmath>
#include <ctime>
#include <queue>
#include <stack>
#include <map>
#include <set>

#define fi first
#define se second
#define U unsigned
#define P std::pair
#define Re register
#define LL long long
#define pb push_back
#define MP std::make_pair
#define all(x) x.begin(),x.end()
#define CLR(i,a) memset(i,a,sizeof(i))
#define FOR(i,a,b) for(Re int i = a;i <= b;++i)
#define ROF(i,a,b) for(Re int i = a;i >= b;--i)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl

const int MAXN = 2000000+5;
const int ha = 1e9 + 7;
bool p[MAXN];
int prime[MAXN],mu[MAXN],cnt;

inline int qpow(int a,int n=ha-2){
    int res = 1;
    while(n){
        if(n & 1) res = 1ll*res*a%ha;
        a = 1ll*a*a%ha;
        n >>= 1;
    }
    return res;
}

inline void pre(int n=MAXN-5){
    mu[1] = 1;p[1] = 1;
    FOR(i,2,n){
        if(!p[i]) prime[++cnt] = i,mu[i] = -1;
        for(int j = 1;j <= cnt && i*prime[j] <= n;++j){
            p[i*prime[j]] = true;
            if(i%prime[j]) mu[i*prime[j]] = -mu[i];
            else{
                mu[i*prime[j]] = 0;
                break;
            }
        }
    }
    FOR(i,1,n) mu[i] = (mu[i-1]+mu[i]+ha)%ha;
}

std::unordered_map<int,int> S;

inline int calc(int x){
    if(x <= MAXN-10) return mu[x];
    if(S[x]) return S[x];
    int ans = 1;
    for(int l = 2,r;l > 0 && l <= x;l = r+1){
        r = x/(x/l);
        ans = (ans-1ll*calc(x/l)*(r-l+1)%ha+ha)%ha;
    }
    return S[x] = ans;
}

int inv6 = qpow(6),inv2 = qpow(2);

inline int get(LL x){
    int res = 1ll*x*(x+1)%ha*(2*x+1)%ha;
    res = 1ll*res*inv6%ha;
    res = (1ll*res+2ll*x)%ha;
    int t = 3ll*x%ha*(x+1)%ha*inv2%ha;
    res = (res-t+ha)%ha;
    return res;
}

inline int get(int l,int r){
    return (get(r)-get(l-1)+ha)%ha;
}

inline void Solve(){
    int ans = 0,x;scanf("%d",&x);
    for(int l = 1,r;l <= x;l = r+1){
        r = x/(x/l);
        int t = 1ll*calc(x/l)*get(l,r)%ha;
       // DEBUG(l);DEBUG(r);DEBUG(t);
        ans = ((ans+t)%ha+ha)%ha;
    }
    printf("%d\n",ans);
}

int main(){
    pre();
    int T;scanf("%d",&T);
    while(T--) Solve();
    return 0;
}