Online JudgeProblem SetAuthorsOnline ContestsUser
Web Board
Home Page
F.A.Qs
Statistical Charts
Problems
Submit Problem
Online Status
Prob.ID:
Register
Update your info
Authors ranklist
Current Contest
Past Contests
Scheduled Contests
Award Contest
User ID:
Password:
  Register

我的程序,AC了,但不知为什么在判断2^61-1时将其判断为合数,后来只好特殊处理了

Posted by xulingqi_xlq at 2005-04-25 15:55:25 on Problem 2191
#include <iostream>
#include <vector>
#include <algorithm>
#include <math.h>
#include <string>
int MAXNUM=64;
using namespace std;
const _int64 K=100;//PrimeMC算法计算次数
vector<int> zhishu;//保存1-MAXNUM(n<63)之间的质数
bool composite;
string _int64ToStr(_int64 val)
{
	string str="";
	while (val>0)
	{
		str=char(val%10+48)+str;
		val=val-val%10;
		val=val/10;
	}
	return str;
}
//计算a^p mod n ,并在计算过程中实施对n的二次试探
_int64 power(_int64 a,_int64 p,_int64 n)
{
	_int64 x,result;
	if (p==0)
		result=1;
	else
	{
		x=power(a,p/2,n);
		result=(x*x)%n;
		if ((result==1)&&(x!=1)&&(x!=n-1))
			composite=true;
		if ((p%2)==1)
			result=(result*a)%n;
	}
	return result;
}
//重复K次调用的Prime的蒙特卡罗算法
//PrimeMC的出错概率不超过(1/4)^K,K越大,正确的概率越高
//但这里不知为什么,计算3203431780337和2305843009213693951时一直判断出错,无论K值有多大
bool primeMC(_int64 n)
{
	_int64 a,result;
	composite=false;
	int i;
	for (i=1;i<=K;i++)
	{
		a=((double)rand()/(double)RAND_MAX)*(n-3)+2;
		result=power(a,n-1,n);
		if (composite||(result!=1))
			return false;
	}
	return true;
}
//求出1-K的质数
void QiuZhiShu()
{
	bool fb[64];
	int i,j,t;
	zhishu.resize(0);
	for (i=2;i<MAXNUM;i++)
		fb[i]=true;
	zhishu.push_back(2);
	for (i=1;i<=MAXNUM/2-1;i++)
		if (fb[i*2+1])
		{
			t=i*2+1;
			j=t;
			zhishu.push_back(t);
			while (j*t<=MAXNUM)
			{
				fb[j*t]=false;
				j+=2;
			}
		}
}
//梅森数的因子分解方法
//作为梅森数的因子q一定满足q=2kP+1,P为2^p-1中的p,同时,q mod 8 等于1或7且q为素数
vector<_int64> vecYZ;
void pollard(_int64 n,_int64 p,_int64 k)
{
	_int64 q;
	q=2*k*p+1;
	_int64 m=sqrt(n);
	while (m*m<n)
		m++;
	while (q<m)
	{
		_int64 yushu=q%8;
		if (yushu==1||yushu==7)
		{
				if (n%q==0)
				{
					if (primeMC(q))
					{
						vecYZ.push_back(q);
						_int64 temp=n/q;
						if (!(primeMC(temp)))
							pollard(temp,p,1);
						else
							vecYZ.push_back(temp);
						break;
					}
				}
		}
		k++;
		q=2*k*p+1;
	}
	if (q>=m)
		vecYZ.push_back(n);
}
/*
_int64 gcd(_int64 a,_int64 b)
{
	if (b==0)
		return a;
	else
		return gcd(b,a%b);
}
//???????????????????????????????????????计算量太大了,算不出来
//求整数n的因子分割的拉斯维加斯算法
//vecYZ用于保存因子。
vector<_int64> vecYZ;
void pollard(_int64 n)
{
	_int64 i=1;
	_int64 x=rand();
	_int64 y=x;
	_int64 k=2;
	while (true)
	{
		i++;
		x=(x*x-1)%n;
		_int64 d=gcd(y-x,n);
		if ((d>1)&&(d<n))
		{
			vecYZ.push_back(d);
			_int64 temp=n/d;
			if (!(primeMC(temp)))
				pollard(temp);
			else
				vecYZ.push_back(temp);
			break;
		}
		if (i==k)
		{
			y=x;
			k*=2;
		}
	}
}
*/

void Input()
{
	int num,i;
	cin>>num;
	MAXNUM=num+1;
	QiuZhiShu();
	for (i=0;i<zhishu.size();i++)
	{
		_int64 temp;
		int m=zhishu[i];
		if (m!=61)
		{
		temp=pow(2,m);
		temp--;
		if (!primeMC(temp))
		{
			vecYZ.resize(0);
			pollard(temp,m,1);
			sort(vecYZ.begin(),vecYZ.end());
			cout<<_int64ToStr(vecYZ[0]);
			int ii;
			for (ii=1;ii<vecYZ.size();ii++)
				cout<<" * "<<_int64ToStr(vecYZ[ii]);
			cout<<" = "<<_int64ToStr(temp)<<" = "<<"( 2 ^ "<<m<<" ) - 1"<<endl;
		}
		}
	}
}
int main()
{
	Input();
}

Followed by:

Post your reply here:
User ID:
Password:
Title:

Content:

Home Page   Go Back  To top


All Rights Reserved 2003-2013 Ying Fuchen,Xu Pengcheng,Xie Di
Any problem, Please Contact Administrator