杜教筛与 Powerful Number 筛杂谈

m256i · 2022-07-24 18:03:46 · 算法·理论

本文的所有函数均指数论函数，即定义域为正整数的函数。

前置

你需要先了解的知识。定理只放式子，不给证明。

狄利克雷卷积：(f*g)(n)=\displaystyle\sum_{ij=n}f(i)g(j)=\sum_{d \mid n}f(d)g\left(\dfrac{n}{d}\right)。容易证明其满足交换律和结合律。

杜教筛

原先是写在这篇题解的，直接照搬过来稍微修改一下。

杜教筛用于求满足 f*g=h 的数论函数 f，其中 g 与 h 是另外两个方便计算的数论函数。

我们定义 S(n)=\displaystyle\sum_{i=1}^nf(i)。

杜教筛基于这个基本式子：

\sum_{i=1}^nh(i)=\sum_{i=1}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)

证明：

\begin{aligned} \sum_{i=1}^nh(i)&=\sum_{i=1}^n(f*g)(i)\\ &=\sum_{i=1}^n(g*f)(i)\\ &=\sum_{i=1}^n\sum_{d \mid i}g(d)f\left(\dfrac{i}{d}\right)\\ &=\sum_{d=1}^n\sum_{i=1}^n[d \mid i]g(d)f\left(\dfrac{i}{d}\right)\\ &=\sum_{d=1}^n\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}g(d)f(i)\\ &=\sum_{d=1}^ng(d)\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}f(i)\\ &=\sum_{d=1}^ng(d)S\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)\\ &=\sum_{i=1}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right) \end{aligned}

那么在这个等式两边同时减去 \displaystyle\sum_{i=2}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)，则有：

\begin{aligned} \left(\sum_{i=1}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\right)-\sum_{i=2}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)&=\left(\sum_{i=1}^nh(i)\right)-\sum_{i=2}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\\ g(1)S(n)&=\left(\sum_{i=1}^nh(i)\right)-\sum_{i=2}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\\ S(n)&=\dfrac{1}{g(1)}\left(\left(\sum_{i=1}^nh(i)\right)-\sum_{i=2}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\right) \end{aligned}

这个式子可以直接数论分块，递归解决即可，边界条件为 n=1，此时结果为 f(1)。我们需要在计算过程中将 S(n) 存储下来。

时间复杂度

由于记忆化的存在，杜教筛的时间复杂度不是很方便分析。

我们考虑画出递归树。以 n=12 为例：

加上记忆化后，一些递归树被砍掉了：

尝试从右往左计算，而不是从左往右：

可以看到，这时的递归树只剩一层了！事实上，由于 \left\lfloor\dfrac{\left\lfloor\frac{a}{b}\right\rfloor}{c}\right\rfloor=\left\lfloor\dfrac{a}{bc}\right\rfloor，那么在深层计算的一定在浅层出现。

并且容易发现，改变计算顺序后复杂度不变（因为每个 n 始终会递归恰好一次）。

所以复杂度为：

\sum_{i \in S}\sqrt i

其中 S 是由所有不同的 \left\lfloor\dfrac{n}{i}\right\rfloor 构成的集合。

用和数论分块一样的方法，它渐进等于：

\sum_{i=1}^{\sqrt n}\sqrt i+\sqrt{\dfrac{n}{i}}

显然 \displaystyle\sum_{i=1}^{\sqrt n}\sqrt i 在这里是可以忽略的（因为它不会大于右边）。我们只需要计算：

\sum_{i=1}^{\sqrt n}\sqrt{\dfrac{n}{i}}

根据这篇文章中的技巧，我们将和式转换成定积分，在渐进意义上相等。因此计算：

\begin{aligned} T(n)&\asymp\sum_{i=1}^{\sqrt n}\sqrt{\dfrac{n}{i}}\\ &\sim\int_1^{\sqrt n}\sqrt{\dfrac{n}{x}}\mathrm dx\\ &\sim\int_0^{\sqrt n}\sqrt{\dfrac{n}{x}}\mathrm dx\\ &=\int_0^{\sqrt n}\dfrac{\sqrt n}{\sqrt x}\mathrm dx\\ &=\sqrt n\int_0^{\sqrt n}x^{-0.5}\mathrm dx\\ &=\sqrt n\dfrac{(\sqrt n)^{0.5}}{0.5}\\ &=2n^{0.75} \end{aligned}

因此不加优化的杜教筛复杂度为 \Theta\left(n^{\frac{3}{4}}\right)。

优化 1：预处理

线性筛预处理 S(1),S(2),\cdots,S(c)\;[c \ge \sqrt n]，复杂度为：

\begin{aligned} T(n)&\asymp\Theta(c)+\sum_{i=1}^{\sqrt n}\left[\dfrac{n}{i} \ge c\right]\sqrt{\dfrac{n}{i}}\\ &\asymp c+\sum_{i=1}^{\sqrt n}\left[i \le \dfrac{n}{c}\right]\sqrt{\dfrac{n}{i}}\\ &=c+\sum_{i=1}^{\frac{n}{c}}\sqrt{\dfrac{n}{i}}\\ &\sim c+\int_1^{\frac{n}{c}}\sqrt{\dfrac{n}{x}}\mathrm dx\\ &\sim c+\int_0^{\frac{n}{c}}\sqrt{\dfrac{n}{x}}\mathrm dx\\ &=c+\sqrt n\int_0^{\frac{n}{c}}x^{-0.5}\mathrm dx\\ &=c+\sqrt n\dfrac{\left(\frac{n}{c}\right)^{0.5}}{0.5}\\ &=c+\dfrac{2n}{\sqrt c}\\ &\asymp c+\dfrac{n}{\sqrt c} \end{aligned}

显然，当加号两边相等时结果最小。

\begin{aligned} c&=nc^{-\frac{1}{2}}\\ c^{\frac{3}{2}}&=n\\ c&=n^{\frac{2}{3}} \end{aligned}

因此，当预处理 n^{\frac{2}{3}} 以内的结果时复杂度最小，为 \Theta\left(n^{\frac{2}{3}}\right)。

优化 2：记忆化

我们的记忆化需要用到 map 或哈希，而这会增大常数（如果用 map 还会使复杂度增加一个 \log），这是我们不希望看到的。于是有了一种更简便的方法：加入一个数组 h，记录 h_x=S\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right)，由于当 x \le n^{\frac{2}{3}} 时会直接访问线性筛，h 数组的大小开到 n^{\frac{1}{3}} 即可。

这个 x 如何确定？可以证明，如果 a=\left\lfloor\dfrac{n}{b}\right\rfloor，那么 \left\lfloor\dfrac{n}{a}\right\rfloor 一定可以作为 x。

请注意，这种方法对多组重复数据不友好（因为每次 h 数组都要清空）。

代码

这是 \varphi 函数的代码。

constexpr int MAXN = 1000006; // n^2/3
constexpr ll P = 1e9+7;
bool isp[MAXN];
int primes[MAXN];
ll phi[MAXN];
ll h[1003]; // n^1/3

template <typename T>
inline T mod(T n){
    T ans = n%P;
    return (ans<0)?ans+P:ans;
}

template <typename T>
inline T pow(T x, T y){
    if (x >= P) x %= P;
    if (y >= P-1) y %= P-1;
    T ans = 1;
    while (y){
        if (y&1) ans = (ans*x)%P;
        x = (x*x)%P;
        y >>= 1;
    }
    return ans;
}

void init(){
    memset(isp, 1, sizeof isp);
    int cnt = 0;
    phi[1] = 1;
    for (int i = 2; i < MAXN; ++i){
        if (isp[i]) primes[cnt++] = i, phi[i] = i-1;
        for (int j = 0; j < cnt && i*primes[j] < MAXN; ++j){
            isp[i*primes[j]] = 0;
            if (!(i%primes[j])){
                phi[i*primes[j]] = phi[i]*primes[j]; break;
            }
            phi[i*primes[j]] = phi[i]*phi[primes[j]];
        }
    }
    for (int i = 2; i < MAXN; ++i) phi[i] = (phi[i]+phi[i-1])%P;
}

ll N;

ll sum_phi(ll n){
    if (n < MAXN) return phi[n];
    int x = N/n;
    if (h[x]) return h[x];
    ll ans = (n*(n+1)>>1)%P;
    for (ll l = 2, r; l <= n; l = r+1){
        r = n/(n/l);
        ans = mod(ans-(r-l+1)*sum_phi(n/l));
    }
    return h[x]=ans;
}

（待更新）