| ||||||||||
| Online Judge | Problem Set | Authors | Online Contests | User | ||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Web Board Home Page F.A.Qs Statistical Charts | Current Contest Past Contests Scheduled Contests Award Contest | |||||||||
取前100个质数去做Rabin-Miller的成功率比随机数更高一些,况且一般都用不着100个,前11、12个质数就应该很够了In Reply To:我的程序,AC了,但不知为什么在判断2^61-1时将其判断为合数,后来只好特殊处理了 Posted by:xulingqi_xlq at 2005-04-25 15:55:25 > #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: |
All Rights Reserved 2003-2013 Ying Fuchen,Xu Pengcheng,Xie Di
Any problem, Please Contact Administrator