杜教筛

杜教筛

Ⅰ.杜教筛有什么用

杜教筛可以快速(在低于线性的时间内)求出积性函数的前缀和,比如莫比乌斯函数(mu(i))的前缀和、欧拉函数(phi(i))的前缀和

Ⅱ.如何计算

假设我们需要计算积性函数(g(i))的前缀和(S(i))
可以构造出一个卷积和:

[sum_{i=1}^{n}(f*g)(i) = sum_{i=1}^{n}sum_{d|i}f(i)g(frac{i}{d}) ]

变换一下求和顺序,先枚举因子得到:

[sum_{i=1}^{n}sum_{d|i}f(i)g(frac{i}{d}) = sum_{d=1}^{n}f(d)sum_{i=1}^{lfloorfrac{n}{d} floor}g(i)=sum_{d=1}^{n}f(d)S(lfloor frac{n}{d} floor) ]

所以:

[sum_{i=1}^{n}(f*g)(i) = sum_{d=1}^{n}f(d)S(lfloor frac{n}{d} floor) ]

[Downarrow ]

[f(1)S(n)=sum_{i=1}^{n}(f*g)(i)-sum_{d=2}^{n}f(d)S(lfloor frac{n}{d} floor) ]

只要能够快速计算出卷积和(sum_{i=1}^{n}(f*g)(i)),那么就能够通过数论分块来快速计算(sum_{d=2}^{n}f(d)S(lfloor frac{n}{d} floor))从而快速计算出(f(1)S(n))

一些性质:

[mu * 1 = epsilon ]

[phi * 1 = ID ]

其中(epsilon(n) = [n==1])(ID(n)=n)(1(n)=1)
洛谷P4213 模板

//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
typedef long long int LL;
const int MAXN = 1e6+7;
LL phi[MAXN],mu[MAXN];
vector<int> prime;
map<int,LL> MU,PHI;
void preprocess(){
    mu[1] = 1;  phi[1] = 1;
    for(int i = 2; i < MAXN; i++){
        if(!phi[i]){
            phi[i] = i - 1;
            mu[i] = -1;
            prime.emplace_back(i);
        }
        for(int j = 0; j < (int)prime.size(); j++){
            if(i*prime[j]>=MAXN) break;
            mu[i*prime[j]] = -mu[i];
            phi[i*prime[j]] = phi[i] * (prime[j] - 1);
            if(i%prime[j]==0){
                mu[i*prime[j]] = 0;
                phi[i*prime[j]] = phi[i] * prime[j];
                break;
            }
        }
    }
    for(int i = 1; i < MAXN; i++){
        phi[i] += phi[i-1];
        mu[i] += mu[i-1];
    }
}
LL calphi(int x){
    if(x<MAXN) return phi[x];
    if(PHI.count(x)) return PHI[x];
    LL tot = (1ll+x)*x/2;
    for(int i = 2; i <= x; i++){
        int j = x / (x / i);
        tot -= (j-i+1ll) * calphi(x/i);
        i = j;
    }
    PHI[x] = tot;
    return PHI[x];
}
LL calmu(int x){
    if(x<MAXN) return mu[x];
    if(MU.count(x)) return MU[x];
    LL tot = 1;
    for(int i = 2; i <= x; i++){
        int j = x/(x/i);
        tot -= (j-i+1ll) * calmu(x/i);
        i = j;
    }
    MU[x] = tot;
    return MU[x];
}
void solve(){
    int x; cin >> x;
    cout << calphi(x) << ' ' << calmu(x) << endl;
}
int main(){
    int T; cin >> T;
    preprocess();
    while(T--) solve();
    return 0;
}
原文地址:https://www.cnblogs.com/kikokiko/p/12857799.html