jzp线性筛及其简单应用

前言：

很久以前看过了线性筛，没怎么注意原理，但是后来发现线性筛还有很有用的。。

比如上次做的一道题就需要找出每个数的最小质因子，先筛再找就太慢了。。一看线性筛发现就可以直接在筛的过程中处理出来了！

今天又学习了屌炸天的jzp线性筛，可以在o(n)的时间内求出欧拉函数， 莫比乌斯函数等积性函数

 

原理:

首先jzp线性筛并不是一种新的线性筛。。其实就是jzp大牛对线性筛的一些开发应用

先回忆一下积性函数的定义 若a，b互质 则f(ab)=f(a)*f(b)的函数f 定义为积性函数，不要求a，b互质也满足的称为完全积性函数

欧拉函数和莫比乌斯函数都是积性函数但不是完全积性函数

假如我们要求 欧拉函数f(n)和莫比乌斯函数 mb(n)

显然如果n的所有质因数(p1,p2...)的次数都是1，显然p1,p2....是互质的，满足积性函数定义，则f(n)=f(p1)*f(p2).....同理mb(n)

而如果某个质因数p的次数不为1，假设为k，我们可以看(yy)出 f(p^k)=p^k-p^(k-1)=(p-1)*p^k，同时由mobius函数定义知如果某个质因数次数大于1次，则其函数值为0

那么如何在线性筛中找到次数不为1的质因子呢

我们观察 if(i%prime[j]==0) break; 这句代码，此处要筛的数n =i*prime[j]，而当i%prime[j]==0 时 显然n%(prime[j]*prime[j])==0。

因此可以知道此时在n的质因子中 prime[j]的次数已经大于1了，就可以处理相应的欧拉函数和莫比乌斯函数了！

 

简单应用：

hdu1695

题意：
求[1,n]和[1,m]之间有多少个互质的数

做法：
以前是用容斥做的，但是容斥需要找质因数，再二进制枚举，比较慢

莫比乌斯函数其实就是容斥的系数，所以直接枚举可能出现的约数（其实就是1~n）用莫比乌斯函数求和即可

最后的式子(不判重)为sum(i=1 to n , mb(i)*(n/i)*(m/i));

这里还有一个小优化，由于是整数除法，对于i=[a,n/(n/a)]   n/i都是是一样的 ，比如 100/(21,22...25)都等于4，这样可以提前对莫比乌斯函数求前缀和，直接累加即可

具体实现见代码，大神们证明了这个优化可以把复杂度降到sqrt(n)级别。具体实现起来的确是快多了，hdu直接0ms AC了！

代码：

#include <iostream>
#include <stdio.h>
#include<string.h>
#include<algorithm>
#include<string>
#include<ctype.h>
using namespace std;
#define maxn 100000
bool notprime[maxn+10];
int prime[maxn+10];
int mb[maxn+10];
int f[maxn+10];
long long sum[maxn+10];
int np;
long long n,m;
void jzp()
{
    np=0;
    memset(notprime,0,sizeof(notprime));
    mb[1]=1;
    for(int i=2; i<=maxn; i++)
    {
        if(!notprime[i])
        {
            prime[np++]=i;
            mb[i]=-1;
            //f[i]=i-1;
        }
        for(int j=0; j<np&&i*prime[j]<=maxn; j++)
        {
            notprime[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                mb[i*prime[j]]=0;
                //f[i*prime[j]]=f[i]*prime[j];
                break;
            }
            else
            {
                mb[i*prime[j]]=-mb[i];
                //f[i*prime[j]]=f[i]*(prime[j]-1);
            }
        }
    }
}
int main()
{
    int t,cas=1;
    scanf("%d",&t);
    jzp();
    sum[0]=0;
    for(int i=1;i<=maxn;i++)
    {
        sum[i]=sum[i-1]+mb[i];
    }
    while(t--)
    {
        int a,b,c,d,k;
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        if(k==0)
        {
            printf("Case %d: 0\n",cas++);
            continue;
        }
        n=min(b/k,d/k);
        m=max(b/k,d/k);
        long long ans=0;
        for(int i=1;i<=n;i++)
        {
            int j=n/(n/i);
            ans+=(n/i)*(n/i)*(sum[j]-sum[i-1]);
            i=j;
        }
        ans=-(ans/2);
        for(int i=1;i<=n;i++)
        {
            int j=min(m/(m/i),n/(n/i));
            ans+=(n/i)*(m/i)*(sum[j]-sum[i-1]);
            i=j;
        }
        printf("Case %d: %I64d\n",cas++,ans);
    }
    return 0;
}

 

最后贴jzp筛模板

bool notprime[maxn+10];
int prime[maxn+10];
int mb[maxn+10];   //mobius
int f[maxn+10];      //euler
int np;
void jzp()
{
    np=0;
    memset(notprime,0,sizeof(notprime));
    mb[1]=1;
    for(int i=2;i<=maxn;i++)
    {
        if(!notprime[i])
        {
            prime[np++]=i;
            mb[i]=-1;
            f[i]=i-1;
        }
        for(int j=0;j<np&&i*prime[j]<=maxn;j++)
        {
            notprime[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                mb[i*prime[j]]=0;
                f[i*prime[j]]=f[i]*prime[j];
                break;
            }
            else
            {
                mb[i*prime[j]]=-mb[i];
                f[i*prime[j]]=f[i]*(prime[j]-1);
            }
        }
    }
}