线性筛（欧拉筛）

昨天的考试跪的一塌糊涂：第一题水过，第二题带WA的朴素，最后题忘了特判左端点全跪，分数比起预计得分整整打了个对折啊！

步入正题：线性筛（欧拉筛）

一般的筛法（PPT里叫埃拉托斯特尼筛法，名字异常高贵）的效率是O(NlglgN)（其实很接近O(n)啊！），对于一些例如N=10000000的残暴数据会跪，于是，线性筛登场了…

#include <cstring>
using namespace std;
int prime[1100000],primesize,phi[11000000];
bool isprime[11000000];
void getlist(int listsize)
{
    memset(isprime,1,sizeof(isprime));
    isprime[1]=false;
    for(int i=2;i<=listsize;i++)
    {
        if(isprime[i])prime[++primesize]=i;
         for(int j=1;j<=primesize&&i*prime[j]<=listsize;j++)
         {
            isprime[i*prime[j]]=false;
            if(i%prime[j]==0)break;
        }
    }
}

 

以上是线性筛代码。

就我的理解，线性筛有两个地方与一般筛不同：

1.两层循环的顺序不同（一般筛是第一维prime[i] 第二维j，欧拉筛是第一维i 第二位prime[j]）

2.一行神奇的代码：

            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函数（反演不怎么懂）：

void getlist(int listsize)
{
    memset(isprime,1,sizeof(isprime));
    isprime[1]=false;
    for(int i=2;i<=listsize;i++)
    {
        if(isprime[i])
        {
             prime[++primesize]=i;
             phi[i]=i-1;
         }
         for(int j=1;j<=primesize&&i*prime[j]<=listsize;j++)
         {
            isprime[i*prime[j]]=false;
            if(i%prime[j]==0)
             {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
}

至于为什么这样的递推是成立的，有如下的几个原因：

　　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(p^k)=(p-1)×p^k-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+p₁+p₁²+...+p₁^k₁)×(1+p₂+p₂²+...+p₂^k₂)×...×(1+p_n+p_n²+...+p_n^k_n)。

显然，函数q是积性的，我们只用考虑i与prime[j]不互质的情况，一种办法是用乘法逆元搞，显然这是不科学的，真正的做法如下

q(i×prime[j])=q(i÷prime[j]^num[i])×q(prime[j]^num[i]+1)

其中的prime[j]^num[i]是可以预先存好的。这种做法还是利用了积性函数的性质。

最后附上昨天第二题的代码

#include <iostream>
#include <cstdio>
using namespace std;
inline long long imax(long long a,long long b){return a>b?a:b;}
const int mo=1000000007;
long long q[2100000];
long long prime[11000000],pnum,ynum[11000000],ysum[11000000],c[11000000];
bool isprime[10000001];
void getlist(int listsize)
{
    memset(isprime,1,sizeof(isprime));
    ynum[1]=ysum[1]=c[1]=1;
    isprime[1]=false;
    for(int i=2;i<=listsize;i++)
    {
        if(isprime[i])
        {
            prime[++pnum]=i;
            c[i]=i;ynum[i]=2;
            ysum[i]=((long long)i*(long long)i+1)%mo;
        }
        for(int j=1;j<=pnum&&i*prime[j]<=listsize;j++)
        {
            int num=i*prime[j];
            isprime[num]=false;
            if(i%prime[j]==0)
            {
                c[num]=c[i]*prime[j];
                if(num!=c[num])
                {
                    ynum[num]=(long long)ynum[c[num]]*(long long)ynum[num/c[num]]%mo;
                    ysum[num]=(long long)ysum[c[num]]*(long long)ysum[num/c[num]]%mo;
                }
                else 
                {
                    ynum[num]=(ynum[i]+1)%mo;
                    ysum[num]=(ysum[i]+(long long)num*(long long)num)%mo;
                }
                break;
            }
            else 
            {
                c[num]=prime[j];
                ynum[num]=(long long)ynum[i]*(long long)ynum[prime[j]]%mo;
                ysum[num]=(long long)ysum[i]*(long long)ysum[prime[j]]%mo;
            }
        }
    }
}
int main(int argc, char *argv[])
{
    freopen("math.in","r",stdin);
    freopen("math.out","w",stdout);
    int n;scanf("%d",&n);
    long long a,b,c;scanf("%I64d%I64d%I64d%I64d",&q[1],&a,&b,&c);
    long long biggest=q[1];
    for(int i=2;i<=n;i++)
    {
        q[i]=((long long)q[i-1]*a+b)%c+1;
        biggest=imax(biggest,q[i]);
    }
    getlist(biggest);
    long long ansnum=0,anssum=0;
    for(int i=1;i<=n;i++)
    {
        ansnum=(ansnum+ynum[q[i]])%mo;
        anssum=(anssum+ysum[q[i]])%mo;
        if((q[i]&1)&&q[i]!=1)
        {
            ansnum=(ansnum+1)%mo;
            anssum=(anssum+4)%mo;
        }
    }
    cout<<ansnum<<endl<<anssum<<endl;
    return 0;
}