杜教筛入门

Posted 2021-03-26 pupuvovovovo

以下主要的话都用无序列表表示。

诶，是不是应该先讲背景

有什么好讲的？

问一个积性函数的前缀和，项数到1e10。

前置知识

线性筛，积性函数。

正文

钦定你已经可以再(O(sqrt{n}))的复杂度内求出：

[sum_{i=1}^n{lfloor frac{n}{i} floor}]

[sum_{i=1}^n{i imes lfloor frac{n}{i} floor}]

对于第一个，先枚举小于(sqrt n)的i，得出这段的值；又因为(i)在一段区间内(lfloor frac{n}{i} floor)都是(i_0leq sqrt n)，因此可以求前式和后式。

• 形式化地讲，若(sum_{i=1}^n{f(i) imes g(lfloor frac{n}{i} floor)})中，若能方便的求得(f(i))的前缀和以及(g(x))，就能方便的求得原式。

然后对于积性函数(f(i))，我们想求(S_n=sum_{i=1}^n{f(i)})。

那么，我们找一个积性函数(g)，令(T_n=sum_{dmid n}g(d) imes f(frac{n}{d}))（就是狄利克雷卷积）。则有：
[ egin{aligned} sum_{i=1}^n{T_i} & =sum_{i=1}^n sum_{dmid i}{g(d) imes f(frac{i}{d})} & =sum_{d=1}^n{g(d) imes sum_{dmid i, ileq n}{f(frac{i}{d})}} & =sum_{d=1}^n{g(d) imes S_{lfloor frac{n}{d} floor}} end{aligned} ]

然后，钦定(g(1)=1)，那么就有
[ egin{aligned} S_n & =sum_{d=1}^n{g(d) imes S_{lfloor frac{n}{d} floor}}-sum_{d=2}^n{g(d) imes S_{lfloor frac{n}{d} floor}} & =sum_{i=1}^n{T_i}-sum_{d=2}^n{g(d) imes S_{lfloor frac{n}{d} floor}} end{aligned} ]

如果不管怎么求(S_{lfloor frac{n}{d} floor})的话，发现满足以下两个条件就可以求(S_n)了。

• 可以求g的前缀和
• 可以求T的前缀和

考虑怎么求上式后面的(S_{lfloor frac{n}{d} floor})，由于f是积性函数，必定可以用线性筛筛出前n项，于是可以用线性筛筛了f求出(S_{1..sqrt n})，至于大于(sqrt{n})的下标，可以记搜：因为(lfloor frac{lfloor frac{n}{d} floor}{e} floor=lfloor frac{n}{d imes e} floor)，于是可以记录n/d的d，当(d>sqrt{n})时直接返回结果就行。

• 于是，求前缀和成功转化成用人类智慧求一个g。

模板题

对于(phi)，g为1（常值函数），T为i（自然数序列），用到的结论是(x=sum_{dmid x}{phi_d})。

证明有这样几个方向：

1. 证明(forall imid n,1leq j<i,(i,j)=1)，都有(i_1 imes j_1 eq i_2 imes j_2)。显然是要被叉翻的（其实是我在用3时不知道干嘛了）。
2. 直接通过(phi)的计算公式和积性通过一波推理得到一些式子，化简得到n。（巨佬做法）
3. 证明(forall imid n,1leq j<i,(i,j)=1)，都有(n/i_1 imes j_1 eq n/i_2 imes j_2)。显然，因为互质，所以两个都是最简分数，于是当i,j不同时，分数不可能相等，于是证明了任意一个(xin [1,n])都只能由一个(i,j)对转移而来，即一一对应，证毕。

以上2的证明：
[ egin{aligned} sum_{dmid n}{phi(d)} & =sum{phi(prod_{i=1}^m{p_i^j})} & =sum{prod_{i=1}^m{phi(p_i^j)}} & =prod_{i=1}^m{sum_{j=0}^{a_i}{phi(p_i^j)}} & =prod_{i=1}^m{(sum_{j=1}^{a_i}{(p_i^j-p_i^{j-1})}+1)} & =prod_{i=1}^m{p_i^{a_i}}=n end{aligned} ]

然后写出来就可以了。

以下是一份跑的非常慢的代码模板（洛谷模板题的关键代码）。

const int N = 3000005, nn = 3000000;
struct getSum{
    ll presum[N], aftersum[N];
    bool calced[N];
    inline ll sum(int n,int d,ll sT(int),ll sg(int),ll g(int)){
        if (n / d <= nn) return presum[n / d];
        if (calced[d]) return aftersum[d];
        int nn = n / d;
        ll ans = sT(nn);
        ans -= (sg(nn) - sg(nn / 2)) * presum[1];
        REP(i, 2, floor(sqrt(nn))){
            ans -= g(i) * sum(n, d * i, T, b, g);
            ans -= (sg(nn / i) - sg(Max(nn / (i + 1), i))) * presum[i];
        }
        calced[d] = 1;
        aftersum[d] = ans;
        return ans;
    }
};
int T, n;
int b[N];
ll phi[N], miu[N];
int temp[N / 10], top;
int main(){
    read(T);
    phi[1] = miu[1] = 1;
    REP(i, 2, nn){
        if (!b[i]){
            temp[++top] = i;
            phi[i] = i - 1;
            miu[i] = -1;
        }
        for (int j = 1; j <= top && i * temp[j] <= nn; ++j){
            b[i * temp[j]] = 1;
            if (i % temp[j]){
                phi[i * temp[j]] = phi[i] * (temp[j] - 1);
                miu[i * temp[j]] = -miu[i];
            }
            else{
                phi[i * temp[j]] = phi[i] * temp[j];
                miu[i * temp[j]] = 0;
                break;
            }
        }
    }
    REP(i, 2, nn){
        phi[i] = phi[i] + phi[i - 1];
        miu[i] = miu[i] + miu[i - 1];
    }
    memcpy(van.presum, phi, 8 * (nn + 1));
    memcpy(deep.presum, miu, 8 * (nn + 1));
    while (T--){
        read(n);
        mem(van.calced);
        printf("%lld ", van.sum(n, 1, [](int n) { return 1LL * n * (n + 1) / 2; }, [](int n) { return 1LL * n; }, [](int n) { return 1LL; }));
        mem(deep.calced);
        printf("%lld
", deep.sum(n, 1, [](int n) { return 1LL; }, [](int n) { return 1LL * n; }, [](int n) { return 1LL; }));
    }
    return 0;
}