线性筛（欧拉筛）
昨天的考试跪的一塌糊涂：第一题水过，第二题带WA的朴素，最后题忘了特判左端点全跪，分数比起预计得分整整打了个对折啊！
步入正题：线性筛（欧拉筛）
一般的筛法（PPT里叫埃拉托斯特尼筛法，名字异常高贵）的效率是O(NlglgN)（其实很接近O(n)啊！），对于一些例如N=10000000的残暴数据会跪，于是，线性筛登场了…
1 #include <cstring> 2 using namespace std; 3 int prime[1100000],primesize,phi[11000000]; 4 bool isprime[11000000]; 5 void getlist(int listsize) 6 { 7 memset(isprime,1,sizeof(isprime)); 8 isprime[1]=false; 9 for(int i=2;i<=listsize;i++) 10 { 11 if(isprime[i])prime[++primesize]=i; 12 for(int j=1;j<=primesize&&i*prime[j]<=listsize;j++) 13 { 14 isprime[i*prime[j]]=false; 15 if(i%prime[j]==0)break; 16 } 17 } 18 }
以上是线性筛代码。
就我的理解，线性筛有两个地方与一般筛不同：
1.两层循环的顺序不同（一般筛是第一维prime[i] 第二维j，欧拉筛是第一维i 第二位prime[j]）
2.一行神奇的代码：
14 if(i%prime[j]==0)break;
这行代码神奇地保证了每个合数只会被它的最小素因子筛掉，就把复杂度降到了O(N)。
接下来是证明这个算法正确性的说明：
P.S.因为自己很蠢，从下午开始研究，吃完饭之后看了部电影，百度了下，找到一个讲的超好的网站，它的说法如下：
prime[]数组中的素数是递增的,当i能整除prime[j]，那么i*prime[j+1]这个合数肯定被prime[j]乘以某个数筛掉。
因为i中含有prime[j],prime[j]比prime[j+1]小，即i=k*prime[j]，那么i*prime[j+1]=(k*prime[j])*prime
[j+1]=k’*prime[j]，接下去的素数同理。所以不用筛下去了。因此，在满足i%prime[j]==0这个条件之前以及第一次
满足改条件时,prime[j]必定是prime[j]*i的最小因子。
之前看了网上各种各样的说法，整个理解都有问题，总算是明白了。
显然，线性筛只拿来筛筛素数是很不科学的，它的速度大约是一般筛的3~4倍，在数据量小的时候甚至慢些（用到了mod运算）。
接下来就是它的重量级应用了：求解积性函数
P.S.睡了一觉，精神抖擞！
积性函数f(x)应满足f(a×b)=f(a)×f(b)，a与b应互素。而完全积性函数则应满足对于任意的a与b，前面的等式应成立。
先是欧拉phi函数（反演不怎么懂）：
1 void getlist(int listsize) 2 { 3 memset(isprime,1,sizeof(isprime)); 4 isprime[1]=false; 5 for(int i=2;i<=listsize;i++) 6 { 7 if(isprime[i]) 8 { 9 prime[++primesize]=i; 10 phi[i]=i-1; 11 } 12 for(int j=1;j<=primesize&&i*prime[j]<=listsize;j++) 13 { 14 isprime[i*prime[j]]=false; 15 if(i%prime[j]==0) 16 { 17 phi[i*prime[j]]=phi[i]*prime[j]; 18 break; 19 } 20 phi[i*prime[j]]=phi[i]*(prime[j]-1); 21 } 22 } 23 }
至于为什么这样的递推是成立的，有如下的几个原因：
1.每个合数只会被筛到一次（前面说明过）。
2.当正整数p是素数时，phi[p]=p-1。
3.phi函数是一个积性函数，当a与b互素时，满足phi(a×b)=phi(a)×phi(b)。
4.上述代码中的i永远大于等于prime[j]，又因为prime[j]必然是素数，所以i、prime[j]互素当且仅当i=0 (mod prime[j])。
5.当p是素数时，phi(pk)=(p-1)×pk-1
原因1保证了每个数的phi只会被计算一次，第10行的代码可以由原因2解释，第20行的代码可以由原因2与原因3解释，第17行的代码可以由原因5所解释。
还是觉得这个算法好神奇……
接下来的东西就和昨天的考试有关了，先是求每个数的因数个数：
记每个数因数个数为e(i)，显然函数e(i)是积性函数但不是完全积性函数。
首先，当p是素数时e(i)=2。
对于i与prime[j]互素的情况，直接相乘即可。而对于i与prime[j]不互素的情况，就有点复杂了。
我们要再记录一个num数组，num[i]表示i的最小素因子（即每次递推时的prime[j]）的次数。
对于i与prime[j]不互素的情况，e(i×prime[j])=e(i)÷(num[i]+1)×(num[i]+2)。
num数组的维护也是很方便的：对于所有n为素数或n=i×prime[j]（i与prime[j]互素）的情况，num(i)=1，对于另外一种情况num[i*prime[j]]=num[i]+1。
再是求因数和还有因数平方和（两者差不了多少）
设数n的因数和为q(n)，有q(n)=(1+p1+p12+...+p1k1)×(1+p2+p22+...+p2k2)×...×(1+pn+pn2+...+pnkn)。
显然，函数q是积性的，我们只用考虑i与prime[j]不互质的情况，一种办法是用乘法逆元搞，显然这是不科学的，真正的做法如下
q(i×prime[j])=q(i÷prime[j]num[i])×q(prime[j]num[i]+1)
其中的prime[j]num[i]是可以预先存好的。这种做法还是利用了积性函数的性质。
最后附上昨天第二题的代码
1 #include <iostream> 2 #include <cstdio> 3 using namespace std; 4 inline long long imax(long long a,long long b){return a>b?a:b;} 5 const int mo=1000000007; 6 long long q[2100000]; 7 long long prime[11000000],pnum,ynum[11000000],ysum[11000000],c[11000000]; 8 bool isprime[10000001]; 9 void getlist(int listsize) 10 { 11 memset(isprime,1,sizeof(isprime)); 12 ynum[1]=ysum[1]=c[1]=1; 13 isprime[1]=false; 14 for(int i=2;i<=listsize;i++) 15 { 16 if(isprime[i]) 17 { 18 prime[++pnum]=i; 19 c[i]=i;ynum[i]=2; 20 ysum[i]=((long long)i*(long long)i+1)%mo; 21 } 22 for(int j=1;j<=pnum&&i*prime[j]<=listsize;j++) 23 { 24 int num=i*prime[j]; 25 isprime[num]=false; 26 if(i%prime[j]==0) 27 { 28 c[num]=c[i]*prime[j]; 29 if(num!=c[num]) 30 { 31 ynum[num]=(long long)ynum[c[num]]*(long long)ynum[num/c[num]]%mo; 32 ysum[num]=(long long)ysum[c[num]]*(long long)ysum[num/c[num]]%mo; 33 } 34 else 35 { 36 ynum[num]=(ynum[i]+1)%mo; 37 ysum[num]=(ysum[i]+(long long)num*(long long)num)%mo; 38 } 39 break; 40 } 41 else 42 { 43 c[num]=prime[j]; 44 ynum[num]=(long long)ynum[i]*(long long)ynum[prime[j]]%mo; 45 ysum[num]=(long long)ysum[i]*(long long)ysum[prime[j]]%mo; 46 } 47 } 48 } 49 } 50 int main(int argc, char *argv[]) 51 { 52 freopen("math.in","r",stdin); 53 freopen("math.out","w",stdout); 54 int n;scanf("%d",&n); 55 long long a,b,c;scanf("%I64d%I64d%I64d%I64d",&q[1],&a,&b,&c); 56 long long biggest=q[1]; 57 for(int i=2;i<=n;i++) 58 { 59 q[i]=((long long)q[i-1]*a+b)%c+1; 60 biggest=imax(biggest,q[i]); 61 } 62 getlist(biggest); 63 long long ansnum=0,anssum=0; 64 for(int i=1;i<=n;i++) 65 { 66 ansnum=(ansnum+ynum[q[i]])%mo; 67 anssum=(anssum+ysum[q[i]])%mo; 68 if((q[i]&1)&&q[i]!=1) 69 { 70 ansnum=(ansnum+1)%mo; 71 anssum=(anssum+4)%mo; 72 } 73 } 74 cout<<ansnum<<endl<<anssum<<endl; 75 return 0; 76 }
