poj3641大素数判断（拉宾米勒）和快速模平方（蒙哥马利）

/*题意：给出p和a，如果p不是素数且a^p === a（mod p）的话，输出yes，否者输出no
解题：蒙哥马利模平方计算+拉宾米勒素数判定方法
0msACc++代码*/
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include<algorithm>
using namespace std;
typedef __int64 int_type;
int_type Montgomery(int_type n,int_type p,int_type m)
{ //蒙哥马利法快速计算(n^p)%m的值
int_type k = 1;
n%=m;
while(p!=1)
{
if(0!=(p&1))k=(k*n)%m;
n=(n*n)%m;
p>>=1;
}
return(n*k)%m;
}
int_type g_aPrimeList[] = {2, 3, 5, 7, 11,13,17,19,23};
bool RabbinMillerTest( int_type n )
{
if (n<2)
{
// 小于2的数即不是合数也不是素数
throw 0;
}
const int_type nPrimeListSize=sizeof(g_aPrimeList)/sizeof(int_type);//求素数表元素个数
for(int i=0;i<nPrimeListSize;++i)
{
// 按照素数表中的数对当前素数进行判断
if (n/2+1<=g_aPrimeList[i])
{
// 如果已经小于当前素数表的数，则一定是素数
return true;
}
if (0==n%g_aPrimeList[i])
{
// 余数为0则说明一定不是素数
return false;
}
}
// 找到r和m，使得n = 2^r * m + 1;
int r = 0, m = n - 1; // ( n - 1 ) 一定是合数
while ( 0 == ( m & 1 ) )
{
m >>= 1; // 右移一位
r++; // 统计右移的次数
}
const int_type nTestCnt = 8; // 表示进行测试的次数,错判概率为小于10^-5
for ( int_type i = 0; i < nTestCnt; ++i )
{
// 利用随机数进行测试，
int_type a = g_aPrimeList[ rand() % nPrimeListSize ];
if ( 1 != Montgomery( a, m, n ) )
{
int_type j = 0;
int_type e = 1;
for ( ; j < r; ++j )
{
if ( n - 1 == Montgomery( a, m * e, n ) )
{
break;
}
e <<= 1;
}
if (j == r)
{
return false;
}
}
}
return true;
}
int main()
{
int_type p,a;
while(scanf("%I64d %I64d",&p,&a),p||a)
{
if(((a%p) == Montgomery(a,p,p)) && (false == RabbinMillerTest(p)))
printf("yes\n");
else
printf("no\n");
}
return 0;
}

本文出自 “Qero” 博客，请务必保留此出处http://8590696.blog.51cto.com/8580696/1358923

poj3641大素数判断（拉宾米勒）和快速模平方（蒙哥马利）

原文：http://8590696.blog.51cto.com/8580696/1358923