51nod 1186 质数检测 V2

#include <bits/stdc++.h>
using namespace std;

const int m10=1e9;
const int z10=9;
const int MAXL=4;

struct bnum
{
	int data[MAXL];
	
	void read()
	{
		memset(data,0,sizeof(data));
		char buf[32];
		scanf(" %s",buf);
		int i,j,cnt,len=strlen(buf);
		for(i=len-1,j=1,cnt=0;i>=0;i--,j*=10)
		{
			if(j==m10)
			{
				j=1;
				cnt++;
			}
			data[cnt]+=j*(buf[i]-'0');
		}
	}
	bool operator ==(const bnum &x)
	{
		for(int i=0;i<MAXL;i++)
		{
			if(data[i]!=x.data[i])
				return 0;
		}
		return 1;
	}
	bnum &operator =(const int x)
	{
		memset(data,0,sizeof(data));
		data[0]=x;
		return *this;
	}
	bnum operator +(const bnum &x)
	{
		bnum ret;
		int i,laz=0;
		for(i=0;i<MAXL;i++)
		{
			ret.data[i]=data[i]+x.data[i]+laz;
			laz=ret.data[i]/m10;
			ret.data[i]%=m10;
		}
		return ret;
	}
	bnum operator -(const bnum &x)
	{
		int i,laz=0;
		bnum ret;
		for(i=0;i<MAXL;i++)
		{
			ret.data[i]=data[i]-x.data[i]-laz;
			if(ret.data[i]<0)
			{
				ret.data[i]+=m10;
				laz=1;
			}
			else
				laz=0;
		}
		return ret;
	}
	// assume *this<x*2
	bnum operator %(const bnum &x)
	{
		int i;
		for(i=MAXL-1;i>=0;i--)
		{
			if(data[i]<x.data[i])
			{
				return *this;
			}
			else if(data[i]>x.data[i])
			{
				break;
			}
		}
		return ((*this)-x);
	}
	bnum &div2()
	{
		int i,laz=0,tmp;
		for(i=MAXL-1;i>=0;i--)
		{
			tmp=data[i]&1;
			data[i]=(data[i]+laz)>>1;
			laz=tmp*m10;
		}
		return *this;
	}
	bool is_odd()
	{
		return data[0]&1;
	}
	bool is_zero()
	{
		for(int i=0;i<MAXL;i++)
		{
			if(data[i])
				return 0;
		}
		return 1;
	}
};

bnum multi(bnum &a0,bnum &b0,bnum &p)
{
	bnum tmp=a0,b=b0,ret;
	ret=0;
	while(!b.is_zero())
	{
		if(b.is_odd())
			ret=(ret+tmp)%p;
		tmp=(tmp+tmp)%p;
		b.div2();
	}
	return ret;
}

bnum powmod(bnum &a0,bnum &b0,bnum &p)
{
	bnum tmp=a0,b=b0,ret;
	ret=1;
	while(!b.is_zero())
	{
		if(b.is_odd())
			ret=multi(ret,tmp,p);
		tmp=multi(tmp,tmp,p);
		b.div2();
	}
	return ret;
}

bool miller_rabin_test(bnum &p,int iter)//检测p是否满足费马小定理，检测inter次 
{
	int i,j,small=0,d=0;
	for(i=1;i<MAXL;i++)
	{
		if(p.data[i])
			break;
	}
	if(i==MAXL)
	{
		if(p.data[0]<2)
			return 0;
		if(p.data[0]==2)
			return 1;
		small=1;
	}
	if(!p.is_odd())
		return 0;
	bnum a,s,m,one,pdl;
	one=1;
	s=pdl=p-one;
	while(!s.is_odd())//p-1=2^d * s，计算d和奇数s
	{
		d++;
		s.div2();
	}
	for(i=0;i<iter;i++)
	{
		//建立随机数 a，2≤a ≤p-2
		a=rand();
		if(small)
		{
			a.data[0]=a.data[0]%(p.data[0]-1)+1;
		}
		else
		{
			a.data[1]=a.data[0]/m10;
			a.data[0]=a.data[0]%m10;
		}
		if(a==one)
			continue;
		//m=a^s %p
		m=powmod(a,s,p);
		//若 a^(s*(2^j))%p==1，则p为合数
		for(j=0;j<d&&!(m==one)&&!(m==pdl);j++)
		{
			m=multi(m,m,p);
		}
		if(!(m==pdl)&&j>0)
			return 0;
	}
	return 1;
}

int main()
{
	bnum x;
	x.read();
	if(miller_rabin_test(x,5))
		printf("Yes\n");
	else
		printf("No\n");
}