前言:
很久以前看过了线性筛,没怎么注意原理,但是后来发现线性筛还有很有用的。。
比如上次做的一道题就需要找出每个数的最小质因子,先筛再找就太慢了。。一看线性筛发现就可以直接在筛的过程中处理出来了!
今天又学习了屌炸天的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+];
int prime[maxn+];
int mb[maxn+];
int f[maxn+];
long long sum[maxn+];
int np;
long long n,m;
void jzp()
{
np=;
memset(notprime,,sizeof(notprime));
mb[]=;
for(int i=; i<=maxn; i++)
{
if(!notprime[i])
{
prime[np++]=i;
mb[i]=-;
//f[i]=i-1;
}
for(int j=; j<np&&i*prime[j]<=maxn; j++)
{
notprime[i*prime[j]]=;
if(i%prime[j]==)
{
mb[i*prime[j]]=;
//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=;
scanf("%d",&t);
jzp();
sum[]=;
for(int i=;i<=maxn;i++)
{
sum[i]=sum[i-]+mb[i];
}
while(t--)
{
int a,b,c,d,k;
scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
if(k==)
{
printf("Case %d: 0\n",cas++);
continue;
}
n=min(b/k,d/k);
m=max(b/k,d/k);
long long ans=;
for(int i=;i<=n;i++)
{
int j=n/(n/i);
ans+=(n/i)*(n/i)*(sum[j]-sum[i-]);
i=j;
}
ans=-(ans/);
for(int i=;i<=n;i++)
{
int j=min(m/(m/i),n/(n/i));
ans+=(n/i)*(m/i)*(sum[j]-sum[i-]);
i=j;
}
printf("Case %d: %I64d\n",cas++,ans);
}
return ;
}
最后贴jzp筛模板
bool notprime[maxn+];
int prime[maxn+];
int mb[maxn+]; //mobius
int f[maxn+]; //euler
int np;
void jzp()
{
np=;
memset(notprime,,sizeof(notprime));
mb[]=;
for(int i=;i<=maxn;i++)
{
if(!notprime[i])
{
prime[np++]=i;
mb[i]=-;
f[i]=i-;
}
for(int j=;j<np&&i*prime[j]<=maxn;j++)
{
notprime[i*prime[j]]=;
if(i%prime[j]==)
{
mb[i*prime[j]]=;
f[i*prime[j]]=f[i]*prime[j];
break;
}
else
{
mb[i*prime[j]]=-mb[i];
f[i*prime[j]]=f[i]*(prime[j]-);
}
}
}
}