杜教篩瞎推【學習筆記】

〇、前言

對於 bzoj3944 來說，和莫比烏斯反演等其他知識關系不大，但是 \(\mu\) 函數在自變量較大情況下的前綴和在反演題中也是會被用到的。

接下來通過 bzoj3944 Sum 一題引入杜教篩的思想。

參考資料

• 鈴懸的數學小講堂——杜教篩 by 鈴懸

一、例題

給定一個正整數 \(N(0\le N\le 2^{31}-1)\)，求 \(\sum_{i=1}^n\varphi(i)\) 和 \(\sum_{i=1}^n\mu(i)\)。多組詢問。

1.分析

對於多組詢問求積性函數前綴和，我們可以 \(O(n)\) 線篩預處理，\(O(1)\) 詢問。其中空間復雜度當然也是 \(O(n)\)

2.推導

其實地上本沒有路，走的人多了，也便成了路。—— 魯迅

對於一般的問題，我們要求 \(S(n)=\sum_{i=1}^nf(i)?\)。

構造 \(g(n)\)，\((f*g)(n)=\sum_{d|n}f(d)g\left(\frac nd\right)\)，我們要求的是 \(S(n)\)，式子中可以從右邊分離出 \(f(n)\)，再進一步推導，求出方程左邊的前綴和，有
\[ \begin{aligned} \sum_{i=1}^n(f*g)(i)&=\sum_{i=1}^n\sum_{d|i}f(d)g\left(\frac id\right)\&=\sum_{i=1}^n\sum_{x=1}^i\sum_{xy=i}f(x)g(y)\\end{aligned} \]

因為每對數 \((x,y)\) 最多只會做一次貢獻，所以我們可以直接枚舉 \(x\cdot y\) 的結果來統計貢獻。
\[ \begin{aligned} \sum_{i=1}^n(f*g)(i)&=\sum_{y=1}^n\sum_{xy\le n}f(x)g(y)\&=\sum_{y=1}^ng(y)\sum_{xy\le n}f(x)\&=\sum_{y=1}^ng(y)\sum_{i=1}^{\left\lfloor\frac ny\right\rfloor}f(i) \end{aligned} \]
我們在上面定義了 \(S(n)\)，所以式子繼續化成了
\[ \sum_{i=1}^n(f*g)(i)=\sum_{y=1}^ng(y)S\left(\left\lfloor\frac ny \right\rfloor\right) \]

移項得到
\[ g(1)S(n)=\sum_{i=1}^n(f*g)(i)-\sum_{y=2}^ng(y)S\left(\left\lfloor\frac ny\right\rfloor\right) \]
這時我們看到了可以數論分塊的 \(\left\lfloor\frac ny\right\rfloor\)，看上去像個遞歸過程。假定 \(g(n)\) 的前綴和可以在 \(O(1)\) 的復雜度內求出，且 \((f*g)(i)\) 的前綴和——即 \(\sum_{i=1}^n(f*g)(i)\) 可以在 \(O(1)\sim O(\sqrt n)\) 的時間內求出，則可以保證復雜度。

因此構造 \(g?\) 是一個重點。

3.復雜度

在構造之前，先考慮一下復雜度。

假定在理想狀態下，求出 \(g(n)\) 的前綴和和 \((f*g)(i)\) 的前綴和都是 \(O(1)\) 的。

復雜度分析為
\[ T(n)=\sum_{i=1}^\sqrt n \sqrt{\left\lfloor\frac ni\right\rfloor}=\int_1^\sqrt n\sqrt{\frac nx}\mathrm dx \]
_{式子可以化為 \(\sqrt n\times x^{-\frac 12} \mathrm dx\)，找到函數 \(\sqrt n\times 2x^{\frac 12}=2\sqrt {nx}\)?的導數是它，然後帶入}
\[ \begin{aligned} T(n)&=\int_1^\sqrt n\sqrt n\times x^{-\frac 12}\mathrm dx\&=2\sqrt{n\sqrt n}-2\sqrt n\&=O\left(n^{\frac 34}\right) \end{aligned} \]

同時，我們可以線篩預處理出前面一段的答案。此外，我們從 \(S(n)\) 進入整個遞歸過程，這之後也只會調用自變量為 \(\left\lfloor\frac ny\right\rfloor\) 的 \(S\) 函數。接著，由於 \(\left\lfloor\frac{\left\lfloor\frac ny\right\rfloor}{x}\right\rfloor=\left\lfloor\frac{n}{xy}\right\rfloor\)，自變量只可能是 \(\left\lfloor\frac ni\right\rfloor\) 的 \(O(\sqrt n)?\) 種取值中的一種，因此可以記憶化。

假設我們預處理出了前 \(n^c\) 的答案，那麽後面一部分積分的上界就是 \(n^{1-c}\) 了，因為只有當 \(k<c\) 時需要計算 \(S\left(\frac{n}{n^{k}}\right)=S\left(n^{1-k}\right)\)。

則
\[ \begin{aligned} T(n)&=O(n^c)+\int_1^{n^{1-c}}\sqrt n\times x^{-\frac 12}\mathrm dx\&=O(n^c)+2\sqrt{n\cdot n^{1-c}}-2\sqrt n\&=O(n^c)+O\left(n^{1-\frac c2}\right) \end{aligned} \]
由基本不等式得二者取等時 \(T(n)\) 有最小值，此時 \(c=1-\frac c2,\ c=\frac 23\)，\(T(n)_\min=O\left(n^\frac 23\right)\)。

4.構造

對於一些常見的函數如 \(id(n),1(n),\epsilon(n)\)，\(g(n)\) 的前綴和比較好求，因此在杜教篩時優先考慮這些函數。

對於 \(\varphi(n)\)，由於 \((\varphi*1)(n)=n\)，\((\varphi*1)(n)\) 和 \(1(n)\) 的前綴和都比較好求，所以可以利用它。

把它帶入移項後的式子：
\[ \begin{aligned} g(1)S(n)&=\sum_{i=1}^n(f*g)(i)-\sum_{y=2}^ng(y)S\left(\left\lfloor\frac ny\right\rfloor\right)\1\times S(n)&=\sum_{i=1}^n(\varphi*1)(i)-\sum_{y=2}^n1\times S\left(\left\lfloor\frac ny\right\rfloor\right)\S(n)&=\frac{n(n+1)}2-\sum_{y=2}^nS\left(\left\lfloor\frac ny\right\rfloor\right) \end{aligned} \]
然後遞歸就可以了。

對於 \(S(x)\) 的計算（\(n\) 是常數）

當 \(x\le \left(2^{31}-1\right)^{\frac 23}\) 時，可以直接返回已經存下的值，其余情況根據關於 \(n\) 的分母，即 \(\left\lfloor\frac nx\right\rfloor\) 的值查詢是否記憶過這個答案，此時只需要存 \(O\left(n^\frac 13\right)\) 種狀態。

5.代碼

#include<cstdio>
#include<cstring>
#define ll long long
const int mxm=2000000;
int pri[1<<20],cnt=0,n;
ll phi[2000100],mu[2000100];
bool is[2000100];
ll f[2000100],g[2000100];
bool used[2000];
ll calc(int x)
{
    if(x<=mxm)
        return phi[x];
    int t=n/x;
    if(used[t])
        return g[t];
    used[t]=1;
    ll ans=(ll)(1ll+x)*x/2;
    for(int i=2;i<=x;++i)
    {
        int nxt=x/(x/i);
        ans-=calc(x/i)*(nxt-i+1);
        i=nxt;
        if(i==2147483647)
            break;
    }
    return g[t]=ans;
}
ll calc1(int x)
{
    if(x<=mxm)
        return mu[x];
    int t=n/x;
    if(used[t])
        return g[t];
    used[t]=1;
    ll ans=1;
    for(int i=2;i<=x;++i)
    {
        int nxt=x/(x/i);
        ans-=calc1(x/i)*(nxt-i+1);
        i=nxt;
        if(i==2147483647)
            break;
    }
    return g[t]=ans;
}
int main()
{
    is[0]=is[1]=1;
    phi[1]=1;
    mu[1]=1;
    for(int i=2;i<=mxm;++i)
    {
        if(!is[i])
        {
            pri[++cnt]=i;
            mu[i]=-1;
            phi[i]=i-1;
        }
        for(int j=1;j<=cnt&&i*pri[j]<=mxm;++j)
        {
            is[i*pri[j]]=1;
            if(i%pri[j]==0)
            {
                mu[i*pri[j]]=0;
                phi[i*pri[j]]=pri[j]*phi[i];
                break;
            }
            else
            {
                mu[i*pri[j]]=-mu[i];
                phi[i*pri[j]]=(pri[j]-1)*phi[i];
            }
        }
        mu[i]+=mu[i-1];
        phi[i]+=phi[i-1];
    }
    int T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d",&n);
        memset(used,0,sizeof(used));
        printf("%lld ",calc(n));
        memset(used,0,sizeof(used));
        printf("%lld\n",calc1(n));
    }
    return 0;
}

在數論題中最好通過帶入最大的數（在本題中是 \(2147483647\)）來檢驗程序的運行時間和數據規模。在本題中，當循環到 \(i=2147483647\) 時，++i 就會爆 int，在這裏需要特判。

數組不要開小了！！

二、練習

感覺把反演給忘了，先復習反演。

??????

杜教篩瞎推 學習筆記【杜教篩】【數學】