1. 程式人生 > >Luogu 4213 【模板】杜教篩(Sum)

Luogu 4213 【模板】杜教篩(Sum)

當作杜教篩的筆記吧。

杜教篩

要求一個積性函式$f(i)$的字首和,現在這個東西並不是很好算,那麼我們考慮讓它捲上另外一個積性函式$g(i)$,使$(f * g)$的字首和變得方便計算,然後再反推出這個$f$函式的字首和。

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

把$g(1)S(n)$移過來

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

這個式子就是杜教篩的精髓了。

我們可以先篩出$[1, \sqrt{n}]$區間內的該積性函式的字首和,然後再分塊遞迴求解$(\sqrt{n}, n]$中的該函式的字首和,可以做到$O(n^{\frac{2}{3}})$的優秀的複雜度(並不會這個複雜度的證明)。

應該用一個雜湊表存一下已經計算過的各個$S(n)$的值($unordered\_map$)。

接下來的問題就是考慮如何搭配出這個積性函式$g$了。

模板題

考慮如何計算$\mu$和$\varphi$。

我們知道$\mu * I = \epsilon$,那麼有

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

滑稽吧,$\epsilon$的字首和還不是$1$。

我們又知道$\varphi * I = id$,那麼又有

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

而$\sum_{i = 1}^{n}id(i) = \sum_{i = 1}^{n}i = \frac{i(i + 1)}{2}$。

解決了!

Code:

#include <cstdio>
#include <cstring>
#include <unordered_map>
using namespace std;
typedef long long ll;

const int N = 5e6 + 5;
const int Maxn = 5e6;

int testCase, pCnt = 0, pri[N], mu[N], phi[N];
ll sumMu[N], sumPhi[N];
bool np[N];
unordered_map <int, ll> sMu, sPhi;

template <typename T>
inline void read(T &X) {
    X = 0; char ch = 0; T op = 1;
    for (; ch > '9'|| ch < '0'; ch = getchar())
        if (ch == '-') op = -1;
    for (; ch >= '0' && ch <= '9'; ch = getchar())
        X = (X << 3) + (X << 1) + ch - 48;
    X *= op;
}

inline void sieve() {
    mu[1] = 1, phi[1] = 1;
    for (int i = 2; i <= Maxn; i++) {
        if (!np[i]) pri[++pCnt] = i, phi[i] = i - 1, mu[i] = -1;
        for (int j = 1; j <= pCnt && i * pri[j] <= Maxn; j++) {
            np[i * pri[j]] = 1;
            if (i % pri[j] == 0) {
                phi[i * pri[j]] = phi[i] * pri[j];
                mu[i * pri[j]] = 0;
                break;
            }
            phi[i * pri[j]] = phi[i] * phi[pri[j]];
            mu[i * pri[j]] = -mu[i];
        }
    }
    
    for (int i = 1; i <= Maxn; i++) {
        sumMu[i] = sumMu[i - 1] + mu[i];
        sumPhi[i] = sumPhi[i - 1] + phi[i];
     }
}

ll getPhi(int n) {
    if (n <= Maxn) return sumPhi[n];
    if (sPhi[n]) return sPhi[n];
    ll res = 1LL * n * (n + 1) / 2;
    for (int l = 2, r; l <= n; l = r + 1) {
        r = (n / (n / l));
        res -= 1LL * (r - l + 1) * getPhi(n / l);
    }
    return sPhi[n] = res;
}

ll getMu(int n) {
    if (n <= Maxn) return sumMu[n];
    if (sMu[n]) return sMu[n];
    ll res = 1LL;
    for (int l = 2, r; l <= n; l = r + 1) {
        r = (n / (n / l));
        res -= 1LL * (r - l + 1) * getMu(n / l);
    }
    return sMu[n] = res;
}

int main() {
    sieve();
    read(testCase);
    for (int n; testCase--; ) {
        read(n);
        printf("%lld %lld\n", getPhi(n), getMu(n));
    }
    return 0;
}
View Code

感覺時限特別急,能別開$long \ long$就別開。