現在杜教篩好像是個人都會了 orz。
所以我寫篇博客復習一下，這里不講原理只講優化技巧的和如何使用
先貼兩個學習原理的博客：
http://www.cnblogs.com/abclzr/p/6242020.html
http://www.cnblogs.com/candy99/p/dj_s.html
然后我們來講一下，怎么使用，下面均以要 求和的n為109 為討論前提
常用公式
一般取用來卷積的輔助函數g為恆等函數：g(n)=1 ,n=1,2,3,…………（有時取g(n)=mu(n) 啥的也有奇效，反正是玄學）
這時杜教篩的公式就能化簡為：
我們一般只用紅框里的兩個公式。
杜教篩不要求函數一定要有積性，而是要求兩點
- f(i)在106以內的所有前綴和能在接近O(n)復雜度內全部處理出來。
- 狄利克雷卷積的前綴和——h(n) ，在1e9 內的任意一項h(m)，能在O(sqrt(m)) 以內的復雜度內求出來。
使用步驟一：預處理106 以內的s[i]
這一步如果f(i)是積性函數，可以直接用線性篩法直接搞定，如果不是的話，就要可能要注意分析性質瞎搞了(比如搞一些預處理，使得f(i)可以O(1)求之類的騷操作)。
注意這里只是以106為例，因為后面的計算過程常數比較大，追求速度的話，最好要能處理到2*106。
使用步驟二：寫快速計算h函數第k項的函數 h(int k)
這里求h，有以下幾種方法
-
用公式推出
的通項看是否有優良的性質。可以快速搞出前n項和 - 對打表找規律
- 對
打表找規律
使用步驟三：寫個記憶化搜索計算s(n)
設我們預處理出s[i]的范圍是i∈[1,UB】。
通用的寫法
若之前算過n 我們則返回我們記下的值，因為n比較大不能用數組直接存，最好用哈希表。
若n<=UB,我們返回預處理好的s[i]
否則使用公式
遞歸計算
要注意到(n/i) 只有2*sqrt(n) 不同的取值，直接用一下數論分塊就行了（如果不知道啥是數論分塊的話，你百度一下就行 了）
然后注意把值存到哈希表內。
代碼如下
1 long long S(long long n) 2 { 3 if(n<=UB) 4 { 5 return s[n]; 6 } 7 long long ans=0; 8 ans=table.find(n);///查詢哈希表中是否有n 9 if(ans==-1) 10 { 11 long long i,j; 12 ans=0; 13 for(i=2; i<=n; i++) 14 { 15 j=i; 16 i=n/(n/i); 17 ans+=(i-j+1)*S(n/j); 18 } 19 ans=(h(n)-ans); 20 table.insert(n,ans);/// 將n對應的值存入哈希表 21 } 22 return ans; 23 }
更方便的寫法
但是如果只有單組，或者你不最求速度，這里有通過構造一個完美哈希把哈希表扔掉，且非遞歸的寫法。
首先我發現遞歸的時候是從大遞歸到小的，跟dp原理遞推的方法一樣，我可以反過來從小算到大，這樣就能保證每次轉移時子狀態已經算好了。
然后問題的關鍵在於有多少種不同狀態，可以發現遞歸時，n是不斷被除的，當前狀態即x=n /a1/a2/a3/……/ak 這到底會有多少種不同的取值呢？ 答案是2*sqrt(n)
因為可以證明在取整情況下也滿足 x=n /a1/a2/a3/……/ak = n/(a1*a2*a3*……*ak) 即x=n/m 而n/m在取整的情況下 就2*sqrt(n)種不同的取值.
接着因為我們預處理了n^(2/3) 即106內的結果， 我們可以再扔掉一部分狀態，只留n^(1/3)的狀態用公式計算出來就行了。
到這里,因為每次轉移復雜度是2*sqrt(x)的，可以證明出杜教篩的復雜度是
※提取要計算的狀態
好廢話講了那么多，我們先預處理出，所要計算的n^(1/3)的狀態的位置（即（n/1,n/2,n/3,……n/(1e3)））
這里直接用暴力算出來就行，代碼如下，注意要按大小順序存起來方便后面轉移
1 vector<ll>q; 2 for(int i=1; i<=n; i++) 3 { 4 i=n/(n/i); 5 if(n/i<=UB) 6 break; 7 q.push_back(n/i); 8 }
※構造完美哈希算法
這里只要一行代碼，就可以構造出完美哈希從而扔掉哈希表，然后把s數組多開大sqrt(n)的大小就行了，其中的原理和精髓可看最終代碼
#define id(x) x<=UB?x:(n/(x)+UB)
注：這里可以證明若UB>sqrt(n),對於所有大於UB的x。n/x 的取值一對一映射到[1,n/UB] 區間內
※最終的實現方式
最后總體代碼如下
1 ll cal(ll n) 2 { 3 vector<ll>q; 4 for(int i=1; i<=n; i++) 5 { 6 i=n/(n/i); 7 if(n/i<=UB) 8 break; 9 q.push_back(n/i); 10 } 11 int k; 12 ll ans,i,j,m,t; 13 for(k=q.size()-1; k>=0; k--)///從小到大枚舉狀態 14 { 15 ans=0; 16 m=q[k]; 17 for(i=2; i<=m; i++)///遞推公式，和文章上面哈希表的寫法對比一下就知道寫的是什么了。 18 { 19 j=i; 20 t=(m/i); 21 i=m/t; 22 ans+=(i-j+1)*s[id(t)]; 23 } 24 ans=(h(m)-ans); 25 s[id(m)]=ans; 26 } 27 return s[id(n)]; 28 }
使用步驟四：驗板子
你要想要知道自己板子有沒寫錯，或者寫搓怎么辦？。上OJ驗板子， 一直wa也不知道wa在哪，題目函數太復雜怎么辦？
答：這里要用驗篩法板子的神級函數，f(i)=i. 這個函數前n項和對1e9+7 取模的結果，用腳指頭算都知道是 n*(n+1)/2%（1e9+7） ，很容易知道你板子算出來的結果對不對。中間過程也比較好處理。
- 首先步驟一：求s[i]前106項，這里我就不詳細說明了，s[i]=s[i-1]+i
- 接着步驟二：寫求h的函數,
，emmmmmmmmmmmm這東西能快速算？？？。 這里有一個只有s(n)很快速計算第n項時（要是能快速計算還寫毛杜教篩，所以這個技巧沒卵用，不過能用來驗板子）才能使用的技巧——枚舉d的取值，則有
，在配合數論分塊就能在根號復雜度內算出h(n)了 - 套上面步驟三寫好的代碼，注意加個取模
- 運行代碼，從文件中讀入1000000000。看程序的運行時間和計算結果
示例代碼如下:
1 #include<stdio.h> 2 #include<string.h> 3 #include<vector> 4 #include<algorithm> 5 using namespace std; 6 #define ll long long 7 #define id(x) x<=UB?x:(n/(x)+UB) 8 const int UB=1e6; 9 const int N=100000; 10 const ll mod=1e9+7; 11 ll s[UB+N]; 12 void init()///預處理s的前UB項 13 { 14 for(int i=1;i<=UB;i++) 15 { 16 s[i]=s[i-1]+i; 17 s[i]%=mod; 18 } 19 } 20 ll h(ll n)///h=Σ Σf(d)*[i%d==0] h為狄利克雷卷積的前綴和 21 { 22 ll i,j; 23 ll sum=0; 24 for(i=1;i<=n;i++) 25 { 26 j=i; 27 i=n/(n/i); 28 sum+=(j+i)*(i-j+1)/2%mod*(n/i)%mod; 29 } 30 return sum%mod; 31 } 32 ll cal(ll n) 33 { 34 vector<ll>q; 35 for(int i=1; i<=n; i++) 36 { 37 i=n/(n/i); 38 if(n/i<=UB) 39 break; 40 q.push_back(n/i); 41 } 42 int k; 43 ll ans,i,j,m,t; 44 for(k=q.size()-1; k>=0; k--)///從小到大枚舉狀態 45 { 46 ans=0; 47 m=q[k]; 48 for(i=2; i<=m; i++)///遞推公式，和文章上面哈希表的寫法對比一下就知道寫的是什么了。 49 { 50 j=i; 51 t=(m/i); 52 i=m/t; 53 ans+=(i-j+1)*s[id(t)]%mod; 54 } 55 ans=(h(m)-ans)%mod; 56 if(ans<0) 57 ans+=mod; 58 s[id(m)]=ans; 59 } 60 return s[id(n)]; 61 } 62 int main() 63 { 64 int i,j; 65 ll n=1e9; 66 // scanf("%lld",&n); 67 init(); 68 printf("out=%lld,ans=%lld\n",cal(n),(n*(n+1)/2)%mod); 69 return 0; 70 }
這個這個必須一秒內出來，且答案正確，不然你就涼拉