线性筛法及扩展

---

以下代码中的定义：
mindiv[i]:i的最小质因子
phi[i]:欧拉函数i的值
mindivq[i]:i的最小质因子的个数
d[i]:i的约数个数
sumd[i]:i的约数和
miu[i]:莫比乌斯函数i的值
inv[i]:i在mod n意义下的乘法逆元

---

标准筛法

欧拉筛法，可以保证每个数只被自己最小的质因子筛去，时间复杂度O(n)
两种等价实现
易于理解版：

void solve(int n){
    for(int i=2;i<=n;i++){
        if(!mindiv[i])
            mindiv[i]=pris[++cnt]=i;
        for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
            mindiv[k]=pris[j];
            if(i%pris[j]==0)break;
        }
    }
}

又短又快版（避免了取模运算）：

void solve(int n){
    for(int i=2;i<=n;i++){
        if(!mindiv[i])
            mindiv[i]=pris[++cnt]=i;
        for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
            if(i%pris[j]==0)break;
            mindiv[k]=pris[j];
        }
    }
}

扩展：求积性函数

积性函数定义：对于一个数论函数f(x)，满足：若gcd(a,b)=1，则f(a*b)=f(a)*f(b)的函数f
由于这个性质，我们可以在筛法的同时求出积性函数的值
先看几个简单的积性函数：

id函数

定义id(i)=i，显然满足积性。
这就不用筛了，直接得到值

e函数

又称单位函数，定义e(1)=1,e(other)=0

phi函数

欧拉函数，很有用，可以用来求逆元
定义phi(i)=sigma{e(gcd(j,i))|1<=j<=i}，即1..n中，与n互质的数的数量
易得对于质数，phi(p)=p-1，
对于合数有以下公式：
n=p1^q1*p2^q2*…*pr^qr
phi(n)=n*(1-1/p1)*(1-1/p2)…
将n带入可得phi(n)=π{(pi-1)*pi^(qi-1)}
这就很容易发现它的积性性质
用筛法怎么求呢？
稍微改一下就行了：

void solve(int n){
    for(int i=2;i<=n;i++){
        if(!mindiv[i]){
            mindiv[i]=pris[++cnt]=i;
            phi[i]=i-1;
        }
        for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
            mindiv[k]=pris[j];
            if(i%pris[j]==0){
                phi[k]=phi[i]*pris[j];
                break;
            }
            phi[k]=phi[i]*(pris[j]-1);
        }
    }
}

为什么呢？
1.对于i%pris[j]!=0的项，它的最小质因子一定是pris[j]，则由积性函数性质可得phi[k]=phi[i]*phi[pris[j]]，其中phi[pris[j]]=pris[j]-1
2.对于i%pris[j]==0的项，则说明pris[j]已经在i中已经出现了，而且是最小的，观察phi函数的计算式，可知此时应乘pris[j]

最小质因数的指数

这个就很简单了：

void solve(int n){
    for(int i=2;i<=n;i++){
        if(!mindiv[i]){
            mindiv[i]=pris[++cnt]=i;
            mindivq[i]=1;
        }
        for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
            mindiv[k]=pris[j];
            if(i%pris[j]==0){
                mindivq[k]=mindivq[i]+1;
                break;
            }
            mindivq[k]=1;
        }
    }
}

约数个数

首先约数个数如何求？
分解质因数，得到n=p1^q1*p2^q2*…*pr^qr
则其约数个数=π(qi+1)
这个结论由乘法原理易得，每个质因子有qi+1种选法
积性就不证了，比较显然。
怎么求呢？观察筛法的过程，i是质数时或i%pris[j]!=0时都比较容易，主要是i%pris[j]==0时需要考虑好。
下面来分析一下：
关键在于如何由d(i)转移到d(i*pris[j])
首先筛法保证了pris[j]一定是最小质因子，那么由于i%pris[j]==0，则意味着最小质因子的指数会+1，这一点在上面的mindivq求解过程中也体现了。
+1会导致什么？
很自然会想到d(k)=d(i)/(mindivq[i]+1)*((mindivq[k]=mindivq[i]+1)+1)
就是这样。

void solve(int n){
    for(int i=2;i<=n;i++){
        if(!mindiv[i]){
            mindiv[i]=pris[++cnt]=i;
            mindivq[i]=1;
            d[i]=2;
        }
        for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
            mindiv[k]=pris[j];
            if(i%pris[j]==0){
                mindivq[k]=mindivq[i]+1;
                d[k]=d[i]/(mindivq[i]+1)*(mindivq[k]+1);
                break;
            }
            mindivq[k]=1;
            d[k]=d[i]*d[pris[j]];
        }
    }
}

约数和

约数和又如何求？
分解质因数，得到n=p1^q1*p2^q2*…*pr^qr
则其约数和=π i ( sigma j {pi^j|0<=j<=qi} )
看起来也是比较显然，我们直接讨论如何利用筛法来求：
依然是利用和上面的类似思路：
只讨论i%pris[j]==0的情况：
需要除以sigma(pris[j]^(0..mindivq[i]))
再乘以sigma(pris[j]^(0..mindivq[k])
其中mindivq[k]=mindiv[i]+1
这样开两个辅助数组记录
t1[i]=sigma(mindiv[i]^(0..mindivq[i]))
t2[i]=mindiv[i]^mindivq
就可以做了：

void solve(int n){
    for(int i=2;i<=n;i++){
        if(!mindiv[i]){
            mindiv[i]=pris[++cnt]=i;
            sumd[i]=1+i;
            t1[i]=1+i;
            t2[i]=i;
        }
        for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
            mindiv[k]=pris[j];
            if(i%pris[j]==0){
                t2[k]=t2[i]*pris[j];
                t1[k]=t1[i]+t2[k];
                sumd[k]=sumd[i]/t1[i]*t1[k];
                break;
            }
            sumd[k]=sumd[i]*sumd[pris[j]];
            t1[k]=1+pris[j];
            t2[k]=pris[j];
        }
    }
}

miu函数

莫比乌斯函数，定义miu(1)=1，miu(i|i is square free number)=0，miu(other)=(-1)^(质因数个数)
用筛法求也很简单：

void solve(int n){
    for(int i=2;i<=n;i++){
        if(!mindiv[i]){
            mindiv[i]=pris[++cnt]=i;
            miu[i]=-1;
        }
        for(int j=1,k;j<=cnt&&(k=i*pris[j])<=n;j++){
            mindiv[k]=pris[j];
            if(i%pris[j]==0){
                miu[k]=0;
                break;
            }
            miu[k]=-miu[i];
        }
    }
}

乘法逆元

由求逆元的欧拉定理方法易知乘法逆元也是积性的，而且有一个很优美的性质：它是完全积性的。
所以求它变得很简单：
inv[p]=p^(phi(n)-1)%n
inv[k]=inv[i]*inv[pris[j]]，gcd(k,n)=1

然而乘法逆元有一个递推的方法，更加简单：
首先逆元满足以下等式（定义）：
inv[P mod i]*(P mod i)=1
将P mod i展开：
inv[P mod i]*(P-P/i*i)=1
进一步展开：
inv[P mod i]i(-P/i)=1
会发现一个神奇的等式：
i*(-P/i)*(inv[P mod i])=1
此时(-P/i)*(inv[P mod i])就是i的逆元！

总结

线性欧拉筛法简洁，积性函数性质优美，充分结合才能发挥更大的效果。同时又要注意到不同函数的性质，用筛法求积性函数的本质是根据积性函数的定义和运算式，利用筛法的特点，构造合适的转移方法。