注：本篇博客是從我知乎搬過來的，一方面是blog的排版不知道比知乎高到哪里去了，另外感覺知乎也不太適合發這種較為理論的內容，遂轉戰博客啦。
這是我的第一篇博客（順便學習了各種格式和排版技巧），大家多多包涵！
最后，轉載請注明出處，謝謝。
1.寫在前面
自從學了杜教篩，突然對這一類問題產生了濃厚的興趣。於是花了好幾天時間，期間也學習了數論相關的許多理論，總算是努力沒有白費，把這一算法的時間復雜度分析出來了。由於網上並沒有這一篩法的相關介紹（僅有的一篇文獻已在附錄列出），在此寫一篇詳盡的分析，也便於以后忘記時查閱。
為了敘述方便，引入以下記號：
1. \(\left\lfloor x \right\rfloor \) 和 \(\left\lceil x \right\rceil \) 分別表示 \(x \) 的下取整和上取整，\([condition]\)表示條件 \(condition\)是否滿足，若滿足則值為1，不滿足為0；
2.設需要考慮的任意積性函數為 \( f(x) \) ，前綴和記為 \( S(n)=\sum\limits_{i = 1}^n {f(i)} \) ，並令 \(m = \left\lceil {\sqrt n } \right\rceil \)；
3.以下用 \( p \) 表示質數，並設 \( maxfactor_i \) 表示 \( i \) 的最大質因子，特別地，\( maxfactor_p=p \)；
4.用 \(i = \prod\limits_{j = 1}^s {p_j^{{k_j}}} \)表示正整數 \(i\) 的質因數分解，其中 \(p_j\) 按升序排列， \(k_j\ge 1\) 。
2.通用篩法原理
2.1 總體思想
我們按如下的算法（以后簡稱算法1）分類進行統計。
在該算法中，我們從2開始從小到大遍歷每一整數 \(i\) ，設 \(F=maxfactor_i\) ，我們進行如下2類操作：
操作一：對每一大於 \(F\) 的質數 \(p\) ，計算 \(f(ip)=f(i)f(p)\) ，並累加進和中，直到 \(ip>n\) 為止；
操作二：如果 \(i\) 含有至少2個質因子 \(F\) ，那么將 \(f(i)\) 累加進和中。
那么全部整數 \(i\) 遍歷完后，所得的總和再加上 \(\sum\limits_{p \le n} {f(p)}\) ，以及 \(f(1)=1\) ，便是 \(S(n)\) 。
命題1：上述算法求得的答案正好是 \(S(n)\) 。
證明：只需證明對於每個整數 \(i\) ，都不重不漏的正好統計了一次即可。
（1）對於 \(f(1)\) ，顯然只被統計了一次；
（2）所有 \(f(p)\) 均被統計了恰好一次，因為操作一以及操作二均不會統計 \(f(p)\) ；
（3）其它情況，設 \(i = \prod\limits_{j = 1}^s {p_j^{{k_j}}}\) ，其中 \(s \ge 2\) ，那么分兩種情況：
情況1：若 \(k_s=1\) ，那么操作一中 \(\prod\limits_{j = 1}^{s-1} {p_j^{{k_j}}}\) 必然會轉移到 \(i\) 而將 \(f(i)\) 統計，且其它數不會轉移到 \(i\) ；顯然操作二也不會統計 \(f(i)\) ；因而被恰好統計了一次；
情況2：若 \(k_s>1\) ，那么操作一不會統計 \(f(i)\) ，但操作二會統計，因而被恰好統計了一次。
證畢。
接來下來考慮一個關鍵問題：如何統計操作一中的 \(\sum\limits_{F < p,ip \le n} {f(p)}\) 。為此我們定義 \(S'(x) = \sum\limits_{p \le x} {f(p)}\) ，那么只要轉化為求 \(S(x)\) 即可。
命題2：上述算法只需要用到 \(S_1=\{ S'(x):1 \le x \le m\} \) 以及 \(S_2=\{ S'\left( {\left\lfloor {\dfrac{n}{x}} \right\rfloor } \right):1 \le x \le m\} \) 。
證明： \(\sum\limits_{F < p,ip \le n} {f(p)} =S'\left( {\left\lfloor {\dfrac{n}{i}} \right\rfloor } \right)-S'(F) \)。
對於 \(S'\left( {\left\lfloor {\dfrac{n}{i}} \right\rfloor } \right)\) ，若 \(i\le m\) ，則屬於 \(S_2\) ，否則屬於 \(S_1\) ；
對於 \(S'(F)\) ，由 \(F^2\le ip\le n\) 得 \(F\le m\) ，屬於 \(S_1\) 。
對於最后還需要加上的全體質數，顯然 \(S'(n)\) 屬於 \(S_2\) 。
證畢。
2.2 積性函數素數前綴和算法
現在問題已經歸結為求1到\(t\)內所有素數\(p\)的前綴和了。我們采用如下算法（以下簡稱算法2）進行求解。
對於求 \(S'(i)\) ，采用歐拉篩法的思想，初始時將每個數都視為質數，將 \(f(j),2 \le j \le i\) 全部加到 \(S'(i)\) 中。然后從小到大遍歷每一質數 \(p\) ，令 \(i\) 從 \(n\) 開始，從大到小執行以下操作： \[S'(i) - = {\rm{ }}\left( {S'\left( {\left\lfloor {\dfrac{i}{p}} \right\rfloor } \right){\rm{ }} - {\rm{ }}S'(p - 1)} \right)f\left( p \right)\] 直到 \(i < p^2\) 為止。
定理3：當 \(f(i)\) 是完全積性函數（即 \(f(ij) = f(i)f(j),\forall i,j\) ）時，上述算法求出的每一個 \(S'(i)\) 是正確的。
證明：我們證明一個循環不變式，即算法執行完當前質數 \(p\) 后，對於每一個 \(i\) ， \(S'(i)\) 均扣除且僅扣除了所有的 \(f(x)\) ，其中 \(x\) 是合數且 \(maxfactor_x \le p\) 。
初始時，顯然是正確的。考慮當前的質數 \(p\) ： \[{\rm{ }}\left( {S'\left( {\left\lfloor {\dfrac{i}{p}} \right\rfloor } \right){\rm{ }} - {\rm{ }}S'(p - 1)} \right)f\left( p \right)\\=f(p)\sum\limits_{j = p}^{\left\lfloor {i/p} \right\rfloor } {[maxfacto{r_j} \ge p]f(j)} \\=\sum\limits_{j = p}^{\left\lfloor {i/p} \right\rfloor } {[maxfacto{r_j} = p]f(jp)} \] 因而質數 \(p\) 對應的一輪循環執行完后，所有 \(x\) 是合數且\(maxfactor_x = p\) 的數會被篩去，且僅篩去了這些數。證畢。
然而，上述證明具有極為苛刻的條件， \(f(i)\) 是完全積性函數。如果這一條件不滿足，結論是不成立的。例如， \(f(p)=-1\) （即莫比烏斯函數 \(\mu\) ）就不滿足，可以令 \(n=4\) 執行上述流程進行驗證。
然而好消息是，函數 \(f(p)=p^a\) 滿足這一要求，於是我們可以采取將一般積性函數表示為 \(p^a\) 的線性組合的形式，對絕大多數積性函數均能這么做。例如： \[\mu(p)=-1\] \[\phi(p)=p-1\] \[\sigma_0(p)=2\] \[\sigma_1(p)=p+1\] 等等。我們對每一項 \(p^a\) 執行一次算法，並將結果保存起來，最后將所有項結果線性疊加即可。注意： \(S'(i)\) 僅統計質數的 \(f(i)\) 的和，與 \(f(p^k),k>1\) 完全無關。在這一步中完全不必關心它， \(f(p^k)\) 在算法1中才有用。
3.通用篩法實現
根據上述原理，這里仍分兩部分考慮實現。
3.1 算法2的實現
首先考慮計算 \(S'(i)\) 的部分。根據結論2，我們只需要兩個大小均為 \(m\) 的數組，其中 \(a[i]=S'(i)\) ， \(b[i]=S'\left( {\left\lfloor {\dfrac{n}{i}} \right\rfloor } \right)\) ，這只需要 \(\Theta (m)\) 的空間。
對於篩法的公式，不難發現以下三種情況：
（1）若 \(\dfrac{i}{p} > \left\lfloor {\sqrt n } \right\rfloor\) ，那么 \(S'(i)\) 便對應了 \(b[n/i]\)，而\(S'\left( {\left\lfloor {\dfrac{i}{p}} \right\rfloor } \right)\) 對應了 \(b[n/i*p]\) ；
（2）若 \(i > \left\lfloor {\sqrt n } \right\rfloor \) 且 \(\dfrac{i}{p} \le \left\lfloor {\sqrt n } \right\rfloor \) ，那么 \(S'(i)\) 便對應了 \(b[n/i]\)，而 \(S'\left( {\left\lfloor {\dfrac{i}{p}} \right\rfloor } \right)\) 便對應了 \(a[i/p]\) ；
（3）若 \(i \le \left\lfloor {\sqrt n } \right\rfloor \) ，則 \(S'(i)\) 便對應了 \(a[i]\)，而 \(S'\left( {\left\lfloor {\dfrac{i}{p}} \right\rfloor } \right)\) 便對應了 \(a[i/p]\) 。
由此可見，更新 \(a,b\) 數組不會用到 \(a,b\) 數組之外的元素 \(S'(i)\) 。
3.2 算法1的實現
計算完 \(S'(i)\) 后，根據之前的算法，我們要枚舉每一整數 \(i\) ，計算 \(f(i)\left(S'\left( {\left\lfloor {\dfrac{n}{i}} \right\rfloor } \right)-S'(maxfactor_i)\right)\) 。為了快速枚全部 \(i\) ，可以考慮遞歸實現。
設函數形式為 \(fun(x,p)\) ，表示 \(x\) 的最大質因子是 \(p\) ，那么我們先計算上面的式子，然后轉移到更大的 \(x\) ，即遞歸的調用 \(fun(y,q)\) ，其中 \(y=xq^k\) ， \(q>p\) 是更大質數。若 \(k\ge 2\) ，那么我們順便計算 \(f(y)\) 。這樣我們就能在 \(O(1)\) 的時間復雜度處理一個 \(i\) 了。
最后，采用一個重要的剪枝：若 \(i \times maxfactor_i>n\) ，顯然這樣的 \(i\) 是無需考慮的，因為不會對答案有貢獻。於是在遞歸調用時加上這一條件即可。顯然不加這一優化，我們需要遍歷1到 \(n\) 中的所有 \(i\) ，時間復雜度將退化為線性。
4.通用篩法的時間復雜度證明
終於到了最具有理論價值（研究了最久）的部分了。同樣地，我們將分別考慮兩個算法各自的時間復雜度。為此，需要引入許多數論中的定理，尤其是素數計數相關的定理。
4.1 素數冪和相關引理
引理4：1到 \(n\) 中質數的個數具有估計式 \(\pi(n)=\Theta \left( {\dfrac{n}{{\ln n}}} \right)\) ，同時隱含常數為1。
引理的證明從略。
引理5：\(\sum\limits_{p \le n} {{p^2}} = \Theta \left( {\dfrac{{{n^3}}}{{\ln n}}} \right)\) ，同時隱含常數為\(\dfrac 1 3\)。
引理的證明從略，以上兩引理可在最后參考文獻中查閱。
引理6：對任意整數 \(k\) ，當 \(n\rightarrow\infty\)時，\(\displaystyle \sum\limits_{p \le n} {\dfrac{1}{{{p^k}}}} = \int_1^n {\dfrac{{\mathrm{d}\pi (x)}}{{{x^k}}}} \) 。
該公式稱為Stieltjes積分公式，相關證明可查閱文獻。
推論7： \(\displaystyle \sum\limits_{ p > n} {\dfrac{1}{{{p^2}}}} =\Theta \left( {\dfrac{{1 }}{{n\ln n}}} \right)\)，同時隱含常數為1 。
由於並沒有查到相關的公式，因此以下對該結論證明。
證明：利用引理6得 \[\sum\limits_{p > n} {\dfrac{1}{{{p^2}}}} = \int_n^\infty {\dfrac{{\mathrm{d}\pi (x)}}{{{x^2}}}} = \dfrac{{\pi (x)}}{{{x^2}}}\bigg|_n^\infty + 2\int_n^\infty {\dfrac{{\mathrm{d}x}}{{{x^2}\ln x}}}\] 這一步是在 \(n\rightarrow\infty\)的前提下成立的，此時分部積分中 \(\pi(x)\) 可替換為 \(\dfrac{x}{\ln x}\) 。繼續計算可以得到 \[\sum\limits_{p > n} {\dfrac{1}{{{p^2}}}} =- \dfrac{1}{{n\ln n}} - 2\int_0^{\dfrac{1}{n}} {\dfrac{{\mathrm{d}t}}{{\ln t}}}=-\dfrac{1}{{n\ln n}} - 2Li\left(\dfrac 1 n\right)\] 其中右邊積分用 \(t=\dfrac 1 x\) 替換得到， \(\displaystyle Li(x)=\int_0^{x} {\dfrac{{\mathrm{d}x}}{{\ln x}}}\) 稱為對數積分。根據 \(Li(x)\) 的泰勒展開式， \[Li(x) = \dfrac{x}{{\ln x}}\sum\limits_{k = 0}^\infty {\dfrac{{k!}}{{{{(\ln x)}^k}}}}\] 因而 \[Li\left(\dfrac 1 n\right)\sim-\dfrac 1 {n\ln n}\] 故 \[\sum\limits_{ p > n} {\dfrac{1}{{{p^2}}}} =\Theta \left( {\dfrac{{1 }}{{n\ln n}}} \right)\] 證畢。
4.2 算法2的時間復雜度證明
前面列舉了一堆引理，現在終於到核心部分了。下面定理給出了算法2中的時間復雜度的精確估計。
定理8：計算 \(S'(i)\) 的時間復雜度為 \(\Theta \left( {\dfrac{{{n^{3/4}}}}{{\ln n}}} \right)\) ，精確到常數，其基本語句執行次數為 \( \dfrac{{{32n^{3/4}}}}{{3\ln n}}\) 。
證明：我們將1到 \(\sqrt n\) 中的質數分為兩部分，對應算法2的流程，分別考察其對時間復雜度的貢獻。
第一部分： \(2\le p \le \sqrt m\)
此時對於每一個 \(1\le i\le \sqrt m\) ， \(b[i]\) 均需要花 \(\Theta (1)\) 的時間更新。
而對於 \(p^2\le i\le \sqrt m\) ， \(a[i]\) 均需要花 \(\Theta (1)\) 的時間更新。
因而對於這個 \(p\) ，總時間復雜度為 \(\Theta (2m-p^2)\) 。
第二部分： \(\sqrt m < p \le m\)
此時 \(a[i]\) 不需要繼續更新；
而對於每一個 \(\dfrac n i\ge p^2\) ，即\(1\le i\le \dfrac n {p^2}\) ， \(b[i]\) 均需要花 \(\Theta (1)\) 的時間更新。
因而對於這個 \(p\) ，總時間復雜度為 \(\Theta \left( {\dfrac{n}{{{p^2}}}} \right)\) 。
將兩部分對 \(p\) 求和，得總時間復雜度為 \[T(n)=\sum\limits_{p \le \sqrt m } {(2m - {p^2})} + \sum\limits_{\sqrt m < p \le m} {\left( {\dfrac{n}{{{p^2}}}} \right)} \\= \left(2m\sum\limits_{p \le \sqrt m } 1\right) - \left(\sum\limits_{p \le \sqrt m } {{p^2}}\right) + \left(n\sum\limits_{\sqrt m < p \le m} {\dfrac{1}{{{p^2}}}} \right)\] 根據引理4、引理5和推論7，得 \[T(n)=\Theta\left(2m\dfrac{{\sqrt m }}{{\ln \sqrt m }} - \dfrac{{{{\left( {\sqrt m } \right)}^3}}}{{3\ln \sqrt m }} + n\dfrac{1}{{\sqrt m \ln \sqrt m }}\right) \\=\Theta\left(\dfrac {8m\sqrt m} {3\ln {\sqrt m}} \right)=\Theta\left(\dfrac {32n^{3/4}} {3\ln n} \right)\] 由於上述證明使用的均為等價形式，因而時間復雜度估計使用的是記號 \(\Theta\) 。證畢。
以下列出了一些典型 \(n\) 時，理論分析和實際測試的計算次數。
當 \(n=10^{11}\) 時，計算結果7489萬，測試結果8029萬；
當 \(n=10^{12}\) 時，計算結果3.86億，測試結果4.13億。
可以看出，兩者是十分接近的。差距主要在於上述引理是在 \(n\) 趨於無窮大時才成立。同時，由於已經考慮常數，故 \(n=10^{11}\) 時算法執行時間僅為300到500毫秒左右。
4.3 Dickman函數
接下來考慮算法1。為此，仍需要引入一些引理。
引理9：對於給定的實常數 \(0 < a < 1\) ，令 \(Q(n) = \{ i \le n:maxfactor_i \le {i^a}\}\) ，那么 \(\left| Q(n) \right| \sim n\rho \left(\dfrac 1 a\right)\) ，其中 \(\rho (x)\) 稱為迪克曼函數(Dickman function)。
性質10：對於任意 \(x>1\) ， \(\rho (x)>0\) 且隨 \(x\) 增大迅速衰減，具有近似式 \(\rho (x) \approx {x^{ - x}}\) 。
以上兩條引理證明從略，可在最后參考文獻中查閱。
上述引理表明，當 \(a\) 較小時， \(Q\) 集合中的元素個數急劇衰減；即對於一個大整數 \(i\)，它的最大質因子通常都比較大。
4.4 關於算法1的時間復雜度
現在回到算法1，之前已說明算法需要枚舉每一個 \(i\) ，而通過一個剪枝 \(i\times maxfactor_i\le n\) 可以大大減少枚舉的 \(i\) 的數量，這一原理便可由上面指出。
然而通過一番分析，仍然不能給出較好的時間復雜度證明。最終不幸的事情還是發生了，我們有如下定理：
定理11：令 \(M(n) = \{ i:i \times maxfacto{r_i} \le n\}\) ，對於任意的 \(\alpha<1\) ，有 \(\left| M(n) \right| = \Omega ({n^\alpha })\) 。
證明：根據引理9， 令 \(Q(n) = \{ i \le n:maxfactor_i \le {i^{1-\alpha}}\}\) ，那么 \(\left| Q(n^\alpha) \right| \sim n^\alpha\rho \left(\dfrac 1 {1-\alpha}\right)\) 。對於每一個數 \(x\in Q(n^\alpha)\) ，顯然 \(x \times maxfacto{r_x} \le n\) ，這就意味着1到 \(n\) 中至少有 \(n^\alpha\rho \left(\dfrac 1 {1-\alpha}\right)\) 個數在 \(M(n)\) 中，即 \(\left| M(n) \right| = \Omega ({n^\alpha })\) 。證畢。
盡管理論十分悲觀，但是不妨代入一個 \(\alpha=3/4\) 算一下。可以查閱到 \(\rho (4) \approx 5 \times 10^{-3}\) ，當 \(n=10^{12}\) 時， \(n^\alpha\rho \left(\dfrac 1 {1-\alpha}\right) \approx 5 \times 10^6\) ，仍是十分小的，這正是得益於\(\rho (x)\) 衰減速度非常快，使得在平常使用時該算法仍具有價值。
注意到定理11給出的明顯是一個下界，實際的值要大不少。下面列出了一些典型 \(n\) 時，實際測試結果。
當 \(n=10^{11}\) 時，測試結果3400萬；
當 \(n=10^{12}\) 時，測試結果1.88億。
可以看出，這一值甚至小於算法2的計算量。但是由於遞歸等原因，算法2往往跑的更快。
命題12：如果對每一質數 \(p\leq \sqrt m\) ， \(i \le n\) 中滿足 \(maxfactor_i=p\) 且 \(ip \le n\) 的 \(i\) 平均小於 \(\sqrt n\) 個，那么算法1就具有 \(O\left(\dfrac {n^{3/4}} {\log n} \right)\) 的時間復雜度。
證明：設 \(maxfactor_i=p\) ，將\(p\) 分為兩段分別考慮。
當 \(\sqrt m < p\leq m\) 時，算法1的時間復雜度必然小於 \(\displaystyle \sum\limits_{\sqrt m < p \le m} {\dfrac{n}{{{p^2}}}}\) ，因為對於每個 \(p\) 只有不超過 \(\left\lfloor {\dfrac{n}{{{p^2}}}} \right\rfloor\) 個 \(i\) 滿足 \(maxfactor_i=p\) 。根據推論7，這一求和式不超過 \(O\left(\dfrac {n^{3/4}} {\log n} \right)\) 。
當 \(p\le \sqrt m\) 時，根據引理5， \(p\) 的個數有 \(\Theta \left( {\dfrac{{{n^{1/4}}}}{{\ln n}}} \right)\) ，再根據命題假設，計算量不超過 \(O\left(\dfrac {n^{3/4}} {\log n} \right)\) 。
綜上可知，整個算法1的時間復雜度不超過 \(O\left(\dfrac {n^{3/4}} {\log n} \right)\) 。證畢。
通過打表知，在\(n \le 10^{12}\) 下命題2的假設滿足，因而可以“認為”整個求積性函數前綴和算法的時間復雜度為 \(\Theta \left( {\dfrac{{{n^{3/4}}}}{{\ln n}}} \right)\) 。
5.一些模板題目
HDU5091 求1到\(n\)中素數個數
LOJ6202 求素數的和
51nod1239 求歐拉函數前綴和
51nod1244 求莫比烏斯函數前綴和
以上是4個模板題。其中求歐拉函數的題使用本文的方法時間上完全可以和杜教篩持平。
SPOJ 求\(\sigma_0(i^3)\)之和
這道題采用本文方法便瞬間退化為一道毫無思維難度的題了。
ZOJ5340 求一般意義下的積性函數前綴和
值得注意的是，此題使用本文方法雖然直接退化為了一道模板題，但卻並不能過（若有人使用此方法能過麻煩聯系我~），其原因在於多組數據。而杜教篩在面對多組數據時是很靈活的，只需改變分塊大小（預處理更多的數）即可。因此面對實際問題要選用合適的方法。
6.后記
盡管最后的結果表明算法時間復雜度似乎是線性的（即對於任意\(\alpha < 1\)，該算法時間復雜度均高於\(n^{\alpha}\)），這已經將算法的理論價值大打折扣，但是：一方面，該算法仍具有巨大的應用價值。例如，用該算法求莫比烏斯函數前綴和，在oj上的提交結果為 \(n=10^{10}\) 下運行350ms，已接近杜教篩的運行時間；另一方面，算法2具有嚴格的理論時間復雜度證明（精確到常數），從而具有很強的理論價值。例如，可用來求1到 \(n\) 中的質數個數，或1到 \(n\) 中的全部質數之和等。
個人感覺，算法2是一個非常優美的算法，因為其時間復雜度的證明過程中，恰好最后的三項均為 \(\Theta \left( {\dfrac{{{n^{3/4}}}}{{\ln n}}} \right)\) ，即各個組成部分時間復雜度平衡。這就類似於杜教篩的 \(n^{2/3}\) 預處理來配平時間復雜度一樣，不同之處在於，這里沒有任何人為因素，完全是天然平衡的，而且是三項求和平衡。
總之，希望該算法的價值能夠在日后不斷體現出來，也更希望能有巨佬將算法1改進，使得理論時間復雜度真正達到 \(\Theta \left( {\dfrac{{{n^{3/4}}}}{{\ln n}}} \right)\) !
最后，由於本人才疏學淺，文章中如有任何不當之處，歡迎批評指正，謝謝！