primegenerate.cpp

RSA加密算法的实现 2^1024大的素数实现

//Generate prime number
#include<iostream>
#include<math.h>
#include"function.h"

using namespace std;
const int TIMES=5;//Times for RabbinMillerTest
/*------------------------------------------------------------------
  The famous Montgomery algorithm to calculate (n^p)%m
  instruction: we know that when we calculate n^p,we can use
  the simple skill that n^p = (n^2)^(p/2), but remember to multiply
  n when p is an odd.
  -------------------------------------------------------------------*/
 /*测试中发现的程序要求:无论是用Montgomery还是QuickMod,在迭代计算的过程中
  因为都mod m(namely the public base),几次后n<m,但是和m的位数很可能一样多,
  或少几位,出现了计算n*n,也就是说为了保证计算不越界,m最大只能取64位,这就
  是说我们最大只用到32位的素数!!何必辛辛苦苦的等一个甚至几个小时去算64位的
  素数?但是在没发现这个问题之前,我计算了一次128位的public_base和private_key
  这样最大的数可能达到256位,没关系,如果非想用这组数的话只要把VSIZE改为256
  就OK~了,我进行了测试,但发现并不是base越大越好,这样虽然在某种程序上增强了
  安全性,但因为数据的增长是十分巨大的,加密时效率变得很低.折中考虑,取base为64位
  -----------------------------------------------------------------------------*/
mathInt Montgomery(mathInt n,mathInt p,mathInt m)
{
    mathInt r = n%m;
	mathInt accul = 1;//the odd part accumulated; 注意:要重载!!!
	mathInt mOne = 1;
	mathInt mTwo = 2;

	while(mOne<p)//thar is p>1;
	{
	    if(p.is_odd())//if p is an odd;
			accul = (accul*r)%m;
		r = (r*r)%m;
		p = p/mTwo;//要重载;+ , - , * , / , % , = , <,==;
	}
	return (r*accul)%m;
}
//Turn a mathInt into a binary num.
vector<int> to_binary(mathInt num)
{
	int len = num.length();
	int unit_len=0;
	int temp = BASE;
	while(temp=temp/2)
		unit_len++;

	vector<int> v = num.get_coeff();
	vector<int> result;
	for(int i=0;i<=len-1;i++)
	{
		for(int j=0;j<=unit_len-1;j++)
		{
			result.push_back(v[i]&1);//that is v[j+8*i]%2;
			v[i]=v[i]>>1;//v[j+8*i] = v[j+8*i]/2;
		}
	}
	return result;
}
/*-----------------------------------------------------------------
  Instruction:when we claculate n^p%m, we could 
  turn p to a binary num-->b0b1b2b3....bn,we only
  need to calculate the nonzero part in b0b1b2b3....bn for p^0 = 1
  -----------------------------------------------------------------*/
mathInt QuickMod(mathInt n, mathInt p,mathInt m)//calculate (n^p)%m
{
     vector<int> pow = to_binary(p);
	 int len = pow.size();
	 mathInt remain=1;
	 n = n%m;
	 for(int i= 0;i<=len-1;i++)
	 {
		 if(pow[i]==1)
            remain = (remain*n)%m;
		 n = (n*n)%m;
	 }

     return remain;
}

/* ----------------------------------------------------------------------------------------
   we know that if p is a prime number then a^(p-1)%p=1 for any a;
   when we random chose a for m time,and each time a^(p-1)%p=1,then the 
   probability that p is't a prime is 1/(4^m); here we test five a to 
   make sure that the probability that p is't a prime is less than 1/1000.

   p-1 must be an even, we need to find such s and d,to make p-1=2^s * d
   we can extend a^(p-1) to (a^((2^s)*d) + 1)(a^((2^(s - 1))*d) + 1).....(a^d + 1)(a^d - 1).
   -----------------------------------------------------------------------------------------*/
bool RabbinMillerTest(mathInt num)
{

	if(!num.is_odd())
		return false;//an even can't be a prime;

	mathInt temp = num-1;
	mathInt mTwo = 2;
	int pass = 0;
	int testTimes = 0;
	int s=0;
	mathInt d;
	int power = 2;
	while(!temp.is_odd())
	{ 
		temp = temp/mTwo;
		s++;
	}
	d = temp;


	for(int i=1;i<=TIMES;i++)
	{
		mathInt a = (rand()%100 + 2);
		bool flag = 1;
		mathInt last = QuickMod(a,d,num);
		if(last==1||last==num-1)
		{
			pass++;
			flag = 0;
		}
		if(flag)
			for(int j=1;j<=s-1;j++)
			{
				if(QuickMod(a,d*power,num)==num-1)
			   {
					pass++;
					break;
			   }
			   power = power*2;
			}
		power = 2;
		testTimes++;//用于记录做了几次测试,如果有一次没通过,则不必再测试.
		if(pass<testTimes)
			return false;	
	}
	return true;

}
mathInt GeneratePrime()
{
	mathInt p;
	int count = 0;
	
	//srand(time(NULL));
	p.set_coeff(1,PRI_LEN);
	while(!RabbinMillerTest(p))
	{
		p.set_coeff(1,PRI_LEN);
		//cout<<"The "<<count<<"th time:"<<endl;
		//p.display();
		//count++;
	}
	return p;
}