杜教筛使用心得

现在杜教筛好像是个人都会了 orz。

所以我写篇博客复习一下，这里不讲原理只讲优化技巧的和如何使用

先贴两个学习原理的博客：

http://www.cnblogs.com/abclzr/p/6242020.html
http://www.cnblogs.com/candy99/p/dj_s.html

然后我们来讲一下，怎么使用，下面均以要求和的n为10⁹ 为讨论前提

常用公式

一般取用来卷积的辅助函数g为恒等函数：g(n)=1 ,n=1,2,3,…………（有时取g(n)=mu(n) 啥的也有奇效，反正是玄学）

这时杜教筛的公式就能化简为：

我们一般只用红框里的两个公式。

杜教筛不要求函数一定要有积性，而是要求两点

f(i)在10⁶以内的所有前缀和能在接近O(n)复杂度内全部处理出来。
狄利克雷卷积的前缀和——h(n) ，在1e9 内的任意一项h(m)，能在O(sqrt(m)) 以内的复杂度内求出来。

使用步骤一：预处理10⁶以内的s[i]

这一步如果f(i)是积性函数，可以直接用线性筛法直接搞定，如果不是的话，就要可能要注意分析性质瞎搞了(比如搞一些预处理，使得f(i)可以O(1)求之类的骚操作)。

注意这里只是以10⁶为例，因为后面的计算过程常数比较大，追求速度的话，最好要能处理到2*10⁶。

使用步骤二：写快速计算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) 即10⁶内的结果，我们可以再扔掉一部分状态，只留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]前10⁶项，这里我就不详细说明了，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
",cal(n),(n*(n+1)/2)%mod);
69     return 0;
70 }

测板子示例

这个这个必须一秒内出来，且答案正确，不然你就凉拉