注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

OI之路,漫漫人生

只为梦想,没有理由

 
 
 

日志

 
 

欧拉筛素数,线性筛欧拉函数  

2014-05-07 22:28:27|  分类: 算法杂类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
1.O(N)欧拉筛素数:
scanf("%d",&n);
 for (i=2;i<=n-1;i++)
            {
            if (!p[i])p[i]=prime[++tot]=i;
            for (j=1;(j<=tot)&&(i*prime[j]<=n-1);j++)
                           {
                           p[i*prime[j]]=prime[j];
                           if (!(i%prime[j]))break;
                           }
            }
其中 if (!(i%prime[j]))break;保证了每个数只被它的最小质因数筛一次。
其中p保存i的最小质因数,用来求接下phi,即欧拉函数。先看下容易理解的:
phi【i】是求在1----i-1中与i互质的对数,两个数互质指:GCD(A,B)=1;
1.首先对于任意正整数i,都能被它的质因数拆分。如12=2^2*3^1;
2.假设N=a^k,那么容易知道和N不能组成互质的数的个数为a^k-1,如16,:有4、8、12、16共4^(2-1)个,即a^k/a=a^k-1个。
这是phi的大写和小写:欧拉筛素数,线性筛欧拉函数 - 山哥 - OI之路,漫漫人生

然后是比较难的:

pq的欧拉函数

假设p,q是两个互质的正整数,则pq的欧拉函数为

                Φ(pq) = Φ(p)Φ(q)           gcd(p,q)=1
证明:
∵M= pq,gcd(p,q) =1
∴根据中国余数定理,有
ZM和Zp×Zq之间存在一一映射
所以M的完全余数集中元素的个数等于集合Zp×Zq元素的个数
而后者的元素个数为Φ(p)Φ(q),所以有
Φ(pq) = Φ(p)Φ(q)//
即欧拉函数为积性函数。

任意正整数的欧拉函数

根据前面积性函数和1的结论,很容易得出它的欧拉函数为:

 Φ(n) = ∏piki-1(pi-1)                   pi表示所有n的质因数,ki表示某个质因数个数。
if n为质数再乘N-1;
接下来是例题1:http://www.lydsy.com/JudgeOnline/problem.php?id=2190
通过以上的基本知识,我们可以知道:phi[I]=i*∏(1-1/pi)。1-1/pi即与i互质的概率,那么所有是I的质因数相乘便是1到I-1中与i互质的概率,再乘个i便是个数。
然后我们分析如何线性实现筛欧拉函数:
我们可以发现如下递推式:
if (p[i]==p[i/p[i]])phi[i]=phi[i/p[i]]*p[i];
           else phi[i]=phi[i/p[i]]*(p[i]-1);
当p[I]与P[I/P[I]]相等时,显然答案就是再乘上P[I]这段多出来的质因数,不需要再乘1-1/P[I];
当不相等是需要乘P[I]*(1-1/P[i])=p[i]-1;
贴下代码:
#include<cstdio>
using namespace std;
int i,j,n,ans,tot;
int a[40010],phi[40010],prime[40010],p[40010];
int main()
{
 scanf("%d",&n);
 for (i=2;i<=n-1;i++)
            {
            if (!p[i])p[i]=prime[++tot]=i;
            for (j=1;(j<=tot)&&(i*prime[j]<=n-1);j++)
                           {
                           p[i*prime[j]]=prime[j];
                           if (!(i%prime[j]))break;
                           }
            }
phi[1]=1;
for (i=2;i<=n-1;i++)
           {
           if (p[i]==p[i/p[i]])phi[i]=phi[i/p[i]]*p[i];
           else phi[i]=phi[i/p[i]]*(p[i]-1);
           ans+=phi[i];
           }
ans=ans*2+3;
if (n==1)ans=0;
printf("%d\n",ans);
return 0;
}

例题2
http://www.lydsy.com/JudgeOnline/problem.php?id=2705
此题可以转换一下,求GCD(I,N)=求GCD(X,Y)=1的个数,当然如果我们枚举一个GCD(I,N)的最大公因数k,那么GCD(I/K,N/K)=1的个数再乘个i就是答案了,而求与N/K互质的数就是欧拉函数嘛。于是问题效率转换为O(SQRT(N)),可以通过。

然而我们不可能一开始线性筛好欧拉函数phi,于是运用积性函数的定理和公式:
Φ(n) = ∏piki-1(pi-1)                   
即枚举n的质因数练成起来即可。
直接贴下代码:
#include<cstdio>
using namespace std;
typedef long long ll;
ll n,ans,i;
ll phi(ll p)
 {ll j,a,ok=1,k;
 for (j=2;j*j<=p;j++){
                a=0;
      while (!(p%j)){a++;p/=j;}
      if (a){
      for (k=1;k<=a-1;k++)ok*=j;
       ok*=(j-1);
             }
                          }
 if (p>1)ok*=(p-1);
 return ok;
 }
int main()
{
 scanf("%lld",&n);ans=0;
 for (i=1;i*i<=n;i++)
   if (n%i==0){
      ans+=i*phi(n/i);
      if (i*i!=n)ans+=n/i*phi(i);
              }
printf("%lld\n",ans);
return 0;
}

例题3:
http://www.lydsy.com/JudgeOnline/problem.php?id=2818
与上题类似的,我们只要同时除以素数x,即可转换为求GCD(I,J)=1的形式,比如6,枚举6的素数2,求出6/2=3的phi为1、3和2、3,那么都乘上2就是2、6和4、6。当然乘上2满足条件的不止phi3,phi2成立的也成立。所以我们统计答案的时候只要类似前缀和一搞即可,注意相反*2和1、1的特殊情况。
贴下代码:
#include<cstdio>
#define nn 10000001
using namespace std;
long long ans,n,tot,i,j,prime[nn],p[nn],phi[nn];
int main()
{
 scanf("%lld",&n);tot=0;
 for (i=2;i<=n;i++)
           {
      if (!p[i])p[i]=prime[++tot]=i;
      for (j=1;(j<=tot)&&(i*prime[j]<=n);j++)
                   {
                   p[i*prime[j]]=prime[j];
                   if (!(i%prime[j]))break;
                   }
           }
phi[1]=1;
for (i=2;i<=n;i++)
               if (p[i]==p[i/p[i]])phi[i]=phi[i/p[i]]*p[i];
               else phi[i]=phi[i/p[i]]*(p[i]-1);
for (i=2;i<=n;i++)phi[i]=phi[i]*2;
for (i=2;i<=n;i++)phi[i]+=phi[i-1];
ans=0;
for (i=1;i<=tot;i++)ans+=phi[n/prime[i]];
printf("%lld\n",ans);
return 0;
}
  评论这张
 
阅读(130896)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018