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

#include "stdafx.h"

#define TRUE 1
#define FALSE 0
#define VLONG    unsigned long
//參考書：《計算機密碼學》（第二版）盧開澄 清華大學出版社
//模冪算法，快速計算：x^r(mod p)，p76
VLONG MONmod(VLONG x,VLONG r,VLONG p)
{
	VLONG a,b,c;
	a=x;
	b=r;
	c=1;
	while (b!=0)
	{
		if ((b&0x01)==0)
		{
			b=b/2;
			a=(a*a)%p;
		}
		else
		{
			b=b-1;
			c=(a*c)%p;
		}
	}
	return c;
}
/*
 * Rabin-Miller
 * 這是個很容易且廣泛使用的簡單算法，它基于Gary Miller的部分象法，
 * 有Michael Rabin發展。事實上，這是在NIST的DSS建議中推薦的算法的一個簡化版。 
 * 首先選擇一個代測的隨機數p，計算b，b是2整除p-1的次數。然后計算m，使得n=1+(2^b)m。
 * （1） 選擇一個小于p的隨機數a。
 * （2） 設j=0且z=a^m mod p
 * （3） 如果z=1或z=p-1，那麼p通過測試，可能使素數
 * （4） 如果j>0且z=1, 那麼p不是素數
 * （5） 設j=j+1。如果j<b且z<>p-1,設z=z^2 mod p，然后回到(4)。
 *       如果z=p-1，那麼p通過測試，可能為素數。
 * （6） 如果j=b 且z<>p-1，不是素數
 * 這個測試較前一個速度快。數a被當成證據的概率為75%。這意味著當迭代次數為t時，它產生一個假的素數所花
 * 費的時間不超過1/4^t。實際上，對大多數隨機數，幾乎99.99%肯定a是證據。

 * 實際考慮：
 * 在實際算法，產生素數是很快的。
 * （1） 產生一個n-位的隨機數p
 * （2） 設高位和低位為1（設高位是為了保證位數，設低位是為了保證位奇數）
 * （3） 檢查以確保p不能被任何小素數整除：如3，5，7，11等等。有效的方法
 *       是測試小于2000的素數。使用字輪方法更快
 * （4） 對某隨機數a運行Rabin-Miller檢測，如果p通過，則另外產生一個隨機
 *		 數a,在測試。選取較小的a值，以保證速度。做5次 Rabin-Miller測試如
 *		 果p在其中失敗，從新產生p，再測試。
 */
//參考書：《計算機密碼學》（第二版）盧開澄 清華大學出版社
//模冪算法，素數的米勒-勒賓(MIller-Rabin)概率測試法，p151
//輸入:被 b要檢測的數字,n素數產生的范圍滿足:(2<=b<n-1)
bool TestPrimeNum(VLONG b,VLONG n)
{
	VLONG m=1;
	VLONG l=0;
	VLONG V;
	VLONG i;
	while (m<(n-1)){m=m*2;l++;}
	if ((b<2)||(b>=(n-1))) return (false);
	V=MONmod(b,m,n);
	if (V==1) return (true);
	i=1;
	do
	{
		if (V==(n-1)) return (true);
		if (i==l) return (false);
		V=MONmod(V,2,n);
		i++;
	}
	while (1);
	return(false);
}
//參考書：《計算機密碼學》（第二版）盧開澄 清華大學出版社
//求gcd{n1,n2}的快速算法，快速計算：gcd{n1,n2}=an1+bn2，p132
VLONG gcd(VLONG /*mod*/ n1,VLONG n2)
{
	VLONG c;
	VLONG t;
	if ((n1<0) || (n2<=0)) return 0;
	//s1
	c=0;
	//s2
	while ((n1%2==0)&&(n2%2==0))
	{
		n1=n1/2;
		n2=n2/2;
		c++;
	}
	//s3
	if  (n2%2==0)
	{
		t=n1;
		n1=n2;
		n2=t;
	}

	do
	{		
		//s4
		while (n1%2==0)
		{
			n1=n1/2;
		}
		//s5
		if ((long)(n1-n2)<0)
		{
			t=n1;
			n1=n2;
			n2=t;
		}
		//s6
		n1=(n1-n2)/2;
	}
	while (n1!=0);
	
	return (n2<<c);
}

//參考書：《計算機密碼學》（第二版）盧開澄 清華大學出版社
//求u mod m的逆元素算法，快速計算：gcd{u,m}=au+bm=1，p136(替換getD)
VLONG Reciprocal_u(VLONG /*mod*/ m,VLONG u)
{
	VLONG n1=m;
	VLONG n2=u;
	long b1=0;
	long b2=1;
	long t;
	VLONG q;
	VLONG r;
	do
	{
		q=n1/n2;
		r=n1-q*n2;
		if (r==0) break;
		n1=n2;
		n2=r;
		t=b2;
		b2=b1-q*b2;
		b1=t;
	}
	while (1);

	if (n2!=1)
		return (0);   //不存在
	if (b2<0)
		return (b2+m);
	else
		return (b2);
}
//解密
VLONG decrypt(VLONG c,VLONG n, VLONG d)
{
	VLONG i,g,f;
	if (d%2==0) g=1; else g=c;
	for (i=1;i<=d/2;i++)
	{
		f=c*c % n;
		g=f*g % n;
	}
	return(g);
}
VLONG getD( VLONG e, VLONG PHI)
{
	VLONG u[3]={1,0,PHI};
	VLONG v[3]={0,1,e};
	VLONG q,temp1,temp2,temp3;
	while (v[2]!=0)
	{
		q=floor(u[2]/v[2]);
		temp1=u[0]-q*v[0];
		temp2=u[1]-q*v[1];
		temp3=u[2]-q*v[2];
		u[0]=v[0];
		u[1]=v[1];
		u[2]=v[2];
		v[0]=temp1;
		v[1]=temp2;
		v[2]=temp3;
	}
	if (u[1]<0) 
	     return(u[1]+PHI);
	else 
	     return(u[1]);
}
VLONG get_common_denom(VLONG e, VLONG PHI)
{
	VLONG great,temp,a;
	if (e >PHI)
	{
		while (e % PHI != 0)
		{
			temp= e % PHI;
			e =PHI;
			PHI = temp;
		}
		great = PHI;
	} 
	else
	{
		while (PHI % e != 0)
		{
			a = PHI % e;
			PHI = e;
			e = a;
		}
		great = e;
	}
	return(great);
}

VLONG getE( VLONG PHI)
{
	VLONG great=0, e=2;
	while (great!=1)
	{
		e=e+1;
		great = gcd(e,PHI);
	}
	return(e);
}
void get_prime( VLONG *val)
{
	#define NO_PRIMES 39
	VLONG primes[NO_PRIMES]={3,5,7,11,13,
		                     17,19,23,29,31,
							 37,41,43,53,59,
							 61,67,71,101,103,97,
							 107,109,113,127,131,
							 137,139,149,151,157,
							 163,167,173,179,181,
							 191,193,197};
	VLONG prime,i;
	
	do
	{
		prime=FALSE;
		printf("Enter a prime number >> ");
		scanf("%ld",val);
		for (i=0;i<NO_PRIMES;i++)
			if (*val==primes[i]) 
				prime=TRUE;
	} 
	while (prime==FALSE);
}

int main(int argc, char* argv[])
{
	VLONG a,b,n,e,PHI,d,m,c;
	/*
	a=MONmod(123,7,187);
	a=MONmod(183,7,187);
	a=MONmod(72,7,187);
	a=MONmod(30,7,187);
	a=MONmod(1520,13,2537);
	a=get_common_denom(18,12);
	a=gcd(18,12);
	*/
    //bool result=TestPrimeNum(3,100);
	//result=TestPrimeNum(17,100);
	//result=TestPrimeNum(19,100);
	//result=TestPrimeNum(2001,10000);
	get_prime(&a);
	get_prime(&b);
	
	n=a*b;
	PHI=(a-1)*(b-1);
	e=getE(PHI);
	d= Reciprocal_u(PHI,e);
	printf("Enter input value >> "); 
	scanf("%ld",&m);
	printf("a=%ld b=%ld n=%ld PHI=%ld\n",a,b,n,PHI);
	//模冪算法，快速計算：x^r(mod p)，p76
    //VLONG MONmod(VLONG x,VLONG r,VLONG p)
	c=(VLONG)pow(m,e) % n; /* note, this may overflow with large numbers */
	//c=MONmod(5,173,323);
	//c=MONmod(5,173,323);
	
	c=MONmod(m,e,n);
	/* when e is relatively large */
	printf("encrypt Message is %ld \n",c);
	printf("e=%ld d=%ld c=%ld\n",e,d,c);
	m=decrypt(c,n,d); /* this function required as c to */
	/*the power of d causes an overflow */
	printf("decrypt Message is %ld \n",m);
	scanf("");
	return(0);
}