杜教筛学习笔记

发布于 2019-03-31  306 次阅读


介绍

我们在 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)$$
如果我们想线性筛积性函数,首先这个过程一定是建立在欧拉筛素数上的,然后我们考虑一些性质:
1. 当 $p$ 为质数时,$f(p)$ 的取值
2. 当这个数不是第一次被筛到的时候,取值的特性

我们拿 $\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 (f * g)(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 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#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 S1;
unordered_map 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 \times \phi(i))$$
求上式的值,其中 $n$ 是线性筛跑不过的量级。
我们一眼看不出来 $g$ 用什么好,不妨先假定有了这个函数并代入狄利克雷卷积的形式:$\sum_{d|n} (d * \phi(d)) \times 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 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#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 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;
}

一个OIer。