prime(1).cpp

网友提供的大数计算的程序

// primenumber.cpp : Defines the entry point for the console application.
//

#include <stdafx.h>
#include <stdio.h> 
#include <stdlib.h>
#include <conio.h>
#include <time.h>
#include <math.h>

int add(int m[101],int n[101],int sum[102])			//加函数，100位加100位

{
	int i;
	sum[101]=0;
	for (i=1;i<101;i++)
		sum[i]=m[i]+n[i];
	for (i=1;i<101;i++)
		if (sum[i]>9) 
		{sum[i+1]+=1;
		sum[i]-=10;
		}
	return 0;
}


int minus(int m[101],int n[101],int difference[101])	//减函数，100位减100位

{
	//因暂无必要，所以先忽略判断被减数与减数的大小关系
	int i;
	for (i=1;i<101;i++)
		difference[i]=m[i]-n[i];
	for (i=1;i<101;i++)
		if (difference[i]<0)
		{difference[i+1]-=1;
		difference[i]+=10;
		}
	return 0;
}


int multiply(int m[101],int n[101],int product[201])	//乘函数，100位乘100位

{
	int sum1[102],sum2[102],i,j,k;
	for (i=1; i<=101;i++)
	{
		sum2[i]=0;
	}
	for (j=1;j<101;j++)
	{
		for(i=1;i<101;i++)sum1[i]=0;
		for (i=1;i<=n[j];i++) 
		{
		 add(m,sum1,sum2);
		 for (k=1;k<=101;  k++)
			 sum1[k]=sum2[k];
		}
		for (k=j;k<=100+j;k++) 
		    product[k]=sum2[k-j+1]+product[k];
	}
	for (i=1;i<201;i++)
	{
	 j=product[i]/10;
	 product[i+1]+=j;
	 product[i]-=j*10;
	}
    return 0;
}



int cmp(int a[301], int ab, int b[201], int bb)  //比较两个数的大小
{
	int i, result;
	result = 0;		  // a == b;
	for(i = 300; i > ab; i--)
		if(a[i] != 0) return 1;
	for(i = 0; i < bb; i++)
	{
		if(a[ab - i] > b[bb - i])
		{ 
			result = 1;	   // a > b;
			break;
		}
		if(a[ab - i] < b[bb - i])
		{ 
			result = -1;	   // a < b;
			break;
		}
	}
	return result;
}


void division(int a[301], int b[201], int c[301], int d[201])	//除函数，暂时300位除200位
//c为商，d为余数
{
	int al, bl, at[301];     // al = length of a, so is bl
	int i, j, al2;

	for(i = 0; i < 301; i++)
	{
		c[i] = 0;
		if(i < 202) d[i] = 0;
		at[i] = a[i];
	}
	for(i = 300; i > 0; i--)
		if(a[i] != 0)
		{
			al = i;
			break;
		}

	for(i = 200; i > 0; i--)
		if(b[i] != 0)
		{
			bl = i;
			break;
		}
	al2 = al;
	for(i = 0; i < al - bl + 1; i++)
	{
		while(cmp(at, al2, b, bl)!=-1)
		{
			for(j = 0; j < bl; j++)
				at[al2 - j] -= b[bl - j];
			for(j = 1; j < 301; j++)
				if(at[j] < 0)
				{ 
					at[j] += 10;
					at[j + 1]--;
				}
			c[al - bl + 1 - i]++;
		}
		al2--;
	}
	for(i = 0; i< 201; i++) d[i] = at[i];
}




int random(int p[101])  //随机产生一个100位的数

{
	p[1]=1;            //低位1为保证该素数为奇数
	int i;
	
    for (i=2;i<=100;i++) 
		p[i]=rand()*10/RAND_MAX;
	while (p[100]==0)           //高位不能为0，保证素数达到要求的长度
		p[100]=rand()*10/RAND_MAX;
//	for (i=100;i>=1;i--)
//	{
//		printf ("%d",p[i]);
	
//	}
//    	printf ("\n");
	return 0;
}


int take( int n[55])     //取小于256的所有素数

{

	int i,j,k=3;
	n[1]=2;
	n[2]=3;
	for (i=4;i<254;i++)
	{
		for (j=2;j<=sqrt(i)+1;j++)
		    if (i%j==0) break;
		 if (i%j!=0)
		 {
			 n[k]=i;
		     k++;
		 }
	}
	return 0;
}




int test(int p[101])   //测试大数,用小于256的所有素数分别去除大数
                       //如果通过，可以排除80%的合数  
{
	int num256[55], m[101];
	take (num256);
	int i,j,sign=1;
    int quo[101],arith[101];
	for (i=1;i<=100;i++)
		m[i]=0;
	while (sign)
	{
		for (i=2;i<55;i++)
		{
			if (num256[i]<10) m[1]=num256[i];
			else
			{ 
				if (num256[i]<100) 
				{
					m[1]=num256[i]%10;
					m[2]=num256[i]/10;
				}
				else 
				{
					m[1]=num256[i]%10;
					m[2]=(num256[i]%100)/10;
					m[3]=num256[i]/100;
				}
			}
			division(p,m,quo,arith);
			sign = 1;
			for (j=1;j<=100;j++)
			{
				if(arith[j] != 0)
				{
					sign=0;
					break;
				}
			}
			if (sign==1) break;	
		}

		
		if (sign==1)
		{
			p[1]+=2;
			for (i=1;i<101;i++)
				if (p[i]>9) 
				{
					p[i+1]+=1;
					p[i]-=10;
				}
		}
	}
	return 0;
}



int rabinmiller(int p[101])	   //Rabin-Miller算法，暂时还不对，不可用。

{
	int a[101],pow2b[101],powam[101],p2[101],z[101];
	int product[101],quo[101],arith[101];
	int i,j=0,k,b=1,m,sign=1;
    for (i=1;i<=100;i++)
	{
		a[i]=0;
		pow2b[i]=0;
	}
	for (i=1;i<=100;i++)
		p2[i]=p[i];
	p2[i]-=1;

	while (sign)
	{
		pow2b[1]=1;
		for (i=1;i<=100;i++)
			pow2b[i]*=2;
		for (i=1;i<=100;i++)
			if (pow2b[i]>9)
			{
				pow2b[i+1]+=1;
				pow2b[i]-=10;
			}
		b++;
		division (p2,pow2b,quo,arith);
		for (i=1;i<=100;i++)
		{
			if (arith[i]==0) sign=1;
			else 
			{
				sign=0;
				b--;
				break;
			}
		}
	}
	m=99/pow(2,b);
	random(a);		  //a暂时取四位，选取较小的a，保证较快的运算速度
	
	for (i=5;i<=100;i++)
		a[i]=0;
	for (i=1;i<=100;i++)	
		powam[i]=a[i];
	for (i=1;i<=m;i++)
	{
		multiply (powam,a,product);
		for (k=1;k<=100;k++)
			powam[i]+=product[i];
	}
	division (powam,p,quo,z);
	
	if (z[1]==1) 
	{
		for (i=2;i<=100;i++)
			if (z[i]==0) return 1;	//return 1表示p通过测试
	}
	int pos=100;
	for (i=100;i<=1;i--)
	{
		if (z[i]==0) pos-=1;
		else break;
	}
	if (cmp(z,pos,p2,100)==0) return 1;
	j=1;
	while (j>0&&j<=b)
	{
		pos=100;
		for (i=100;i<=1;i--)
		{
			if (z[i]==0) pos-=1;
			else break;
		}
		if ((j==b)&&(cmp(z,pos,p2,100)!=0)) return 0;

		division(powam,p,quo,z);
		if (z[1]==1)
		{
			for (i=2;i<=100;i++)
				if (z[i]==0) return 0;	//return 0表示p不是素数
		}
		j++;
		pos=100;
		for (i=100;i<=1;i--)
		{
			if (z[i]==0) pos-=1;
			else break;
		}
		if (cmp(z,pos,p2,100)==0) return 1;
		else
		{
			multiply(z,z,product);
			division(product,p,quo,arith);
			for (i=1;i<=100;i++)
				z[i]=arith[i];
		}
	}
}

void ExtBin(int *u,int *v,int *u1,int *u2,int *u3)  //欧几里德算法
{
}



int main(int argc, char* argv[])
{  
	int i;
	int p1[101],p2[101];      //用来产生两个随机数
    time_t t;                 //这两行保证每次产生的随机数不同
	srand( (unsigned) time( &t ) );         
 	random(p1);
    for (i=100;i>=1;i--)
		printf ("%d",p1[i]);
	printf ("\n");
	random(p2);
	for (i=100;i>=1;i--)
		printf ("%d",p2[i]);
	printf("\n");
    printf("\n");
	
	//以下为测试
	int sum[102],difference[101],product[201],quo[101],arith[101];      
	add(p1,p2,sum);                                
	if (sum[101]>=1&&sum[101]<=9) 
		printf ("%d",sum[101]);
	for (i=100;i>=1;i--)
		printf ("%d",sum[i]);
	printf("\n");

	minus(p1,p2,difference);
	for (i=100;i>=1;i--)
		printf ("%d",difference[i]);
	printf("\n");

	for(i=1;i<201;i++)
		product[i]=0;
	multiply(p1,p2,product);
	if (product[200]>=1&&product[200]<=9)
		printf("%d",product[200]);
	for (i=199;i>=1;i--)
		printf ("%d",product[i]);
	printf("\n");    
	
	division(p1,p2,quo,arith);   //因为除法调用的是300位和200位的数，所以暂时不能出结果
	int j=0;
	for (i=100;i>=1;i--)
	{
		if (j==0)
		{
			if (quo[i]==0)	continue;
			else j=1; 
		}
		printf ("%d",quo[i]);
	}
    printf("\n");

	j=0;
	for (i=100;i>=1;i--)
	{
		if (j==0)
		{
			if (arith[i]==0) continue;
			else j=1;
		}
		printf ("%d",arith[i]);
	}
	printf("\n");
	printf("\n");
	//以上为测试
    
	
	int result=0;
	while (result==0)
	{
		test (p1);
		result=rabinmiller(p1);
		if (result==0) p1[1]+=2;
		for (i=1;i<=100;i++)
		{
			if (p1[i]>9) 
			{
				p1[i+1]+=1;
				p1[i]-=10;
			}
		}
	}
	result=0;
	while (result==0)
	{
		test (p2);
		result=rabinmiller(p2);
		if (result==0) p2[1]+=2;
		for (i=1;i<=100;i++)
		{
			if (p2[i]>9) 
			{
				p2[i+1]+=1;
				p2[i]-=10;
			}
		}
	}
    
  for (i=100;i>=1;i--) printf("%d",p1[i]);
   printf("\n");
	
	for (i=100;i>=1;i--) printf("%d",p2[i]);
    printf("\n");
	printf("\n");

	
	
	

	int m[201],n[201];
	multiply(p1,p2,n);
	p1[1]-=1;
	p2[1]-=1;
	multiply(p1,p2,m);
	return 0;
}