線性篩（歐拉篩）

昨天的考試跪的一塌糊塗：第一題水過，第二題帶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;
}