poj Pseudoprime numbers

快速幂取模与伪素数检测算法

最新推荐文章于 2020-02-05 21:07:55 发布

weixin_30492047

最新推荐文章于 2020-02-05 21:07:55 发布

阅读量65

点赞数

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/one--world--one--dream/archive/2011/10/30/2229197.html

本文介绍了一种快速幂取模算法及其在检测伪素数方面的应用，通过结合快速幂取模和欧拉定理优化，实现高效判断大数是否为伪素数。文中详细解释了算法原理、代码实现以及实例分析。

Pseudoprime numbers

Time Limit: 1000MS		Memory Limit: 65536K
Total Submissions: 4496		Accepted: 1716

Description

Fermat's theorem states that for any prime number p and for any integer a > 1, a^p = a (mod p). That is, if we raise a to the p^th power and divide by p, the remainder is a. Some (but not very many) non-prime values of p, known as base-a pseudoprimes, have this property for some a. (And some, known as Carmichael Numbers, are base-a pseudoprimes for all a.)

Given 2 < p ≤ 1000000000 and 1 < a < p, determine whether or not p is a base-a pseudoprime.

Input

Input contains several test cases followed by a line containing "0 0". Each test case consists of a line containing p and a.

Output

For each test case, output "yes" if p is a base-a pseudoprime; otherwise output "no".

Sample Input

Sample Output

no
no
yes
no
yes
yes

这题其实是简单题，我把它想复杂了，直接快速求幂取模0（logn）的时间复查度就能过；我用到了欧拉定理优化了下，把时间复杂度降到了更低；
（1）：用米勒测试法判断n是否为素数；
（2）：a^phi(p)==1 (mod p) 如果gcd(a,p)==1,p是合数，a也可能是合数，所以要用到欧拉定理优化就要把a进行素因子分解，a=a1^p1*a2^p2*....*ak^pk,a^p=(a1^p1*a2^p2*....*ak^pk)^p
=a1^(p1*p)*a2^(p2*p)....a3^(pk*p)=a1^(p%phi(p))^p1*a2^(p%phi(p))^p2...ak^(p%phi(p))^ak;
我的代码

View Code

  1 #include <iostream>
  2 #include <time.h>
  3 #include <cstdlib>
  4 #include <cstdio>
  5 
  6 #define maxTest 100//开始我把它改成30就悲剧的wa了，后来抱着试着看的心理改成100就ac了，无语了，不然怎么找不出错了
  7 using namespace std;
  8 
  9 __int64 Random(__int64 n)
 10 {
 11     return (__int64)((double)rand()/RAND_MAX*n+0.5);
 12 }
 13 __int64 Modular_Exp(__int64 a, __int64 b, __int64 n) // a^b mod n
 14 {
 15     __int64 ans;
 16     if(b == 0)
 17            return 1;
 18     if(b == 1)
 19     return a%n;
 20     ans = Modular_Exp(a, b/2, n);
 21     ans = ans*ans%n;
 22     if(b%2)
 23     ans = ans*a%n;
 24     return ans;
 25 }
 26 
 27 bool Miller_Rabbin(__int64 n)
 28 {
 29     for(int i=1; i<=maxTest; i++)
 30     {
 31        __int64 a = Random(n-2)+1;
 32        if(Modular_Exp(a, n-1, n) != 1)
 33       return false;
 34     }
 35     return true;
 36 }
 37 
 38 __int64 phi(__int64 a)
 39 {
 40     __int64 temp=a;
 41     for(__int64 i=2;i*i<=a;i++)
 42     {
 43         if(a%i==0)
 44         {
 45             temp=temp/i*(i-1);
 46             while(a%i==0)
 47               a/=i;
 48         }
 49     }
 50     if(a>1)  temp=temp/a*(a-1);
 51     return temp;
 52 }
 53 
 54 __int64 muti(__int64 p,__int64 n,__int64 mod)
 55 {
 56     __int64 temp=p,ans=1;
 57     while(n)
 58     {
 59         if(n&1)  ans*=temp;
 60         temp=temp*temp%mod;
 61         ans%=mod;
 62         n>>=1;
 63     }
 64     return ans;
 65 }
 66 int prime[100],num[100];
 67 int main()
 68 {
 69     __int64 p,a;
 70     while(cin>>p>>a&&!(p==0&&a==0))
 71     {
 72         __int64 key=a;
 73         if(Miller_Rabbin(p))
 74           {
 75               puts("no");continue;
 76           }
 77         __int64 mod=p;
 78        // cout<<phi(p)<<endl;
 79         p=p%phi(p);
 80         int num1=0;
 81         for(int i=2;i*i<=a;i++)
 82         {
 83              if(a%i==0)
 84              {
 85                   prime[num1]=i;
 86                   int t=0;
 87                  while(a%i==0)
 88                   {
 89                       t++,a/=i;
 90                   }
 91                   num[num1++]=t;
 92              }
 93         }
 94         if(a>1) prime[num1]=a,num[num1++]=1;
 95         __int64 ans=1;
 96        for(int i=0;i<num1;i++)
 97        {
 98            __int64 temp=muti((__int64)prime[i],p,mod);
 99            temp=muti(temp,(__int64)num[i],mod);
100            ans=ans*temp%mod;
101        }
102        puts(ans==key?"yes":"no");
103     }
104     return 0;
105 }