算法系列2——RSA

1. RSA介绍

RSA公钥加密算法是1977年由Ron Rivest、Adi
Shamirh和LenAdleman在（美国麻省理工学院）开发的。RSA取名来自开发他们三者的名字。RSA是目前最有影响力的公钥加密算法，它能够抵抗到目前为止已知的所有密码攻击，已被ISO推荐为公钥数据加密标准。RSA算法基于一个十分简单的数论事实：将两个大素数相乘十分容易，但那时想要对其乘积进行因式分解却极其困难，因此可以将乘积公开作为加密密钥。

RSA的安全性依赖于大数的因子分解，但并没有从理论上证明破译RSA的难度与大数分解难度等价。即RSA的重大缺陷是无法从理论上把握它的保密性能如何，而且密码学界多数人士倾向于因子分解不是NPC问题。

RSA的缺点主要有：A）产生密钥很麻烦，受到素数产生技术的限制，因而难以做到一次一密。B）分组长度太大，为保证安全性，n 至少也要600BITS以上，使运算代价很高，尤其是速度较慢，较对称密码算法慢几个数量级；且随着大数分解技术的发展，这个长度还在增加，不利于数据格式的标准化。目前，SET(Secure Electronic Transaction）协议中要求CA采用2048BITS长的密钥，其他实体使用1024比特的密钥。C)RSA密钥长度随着保密级别提高，增加很快。下表列出了对同一安全级别所对应的密钥长度。

该可逆算法是经批准用于加密和生成数字签名的算法。公钥指数的值只允许是3和216＋1。该算法产生的密文及数字签名的长度与模长相等。

认证中心公钥模的长度NCA，发卡行公钥模的长度NI，IC卡公钥模的长度NIC，必须满足NIC≤NI≤NCA和NPE≤NI≤NCA。

注： IC卡中一个记录的长度最长不超过254字节（包括Tag和Length），因而实际IC卡公钥和发卡行公钥长度应小于最大长度248字节。命令数据长度最长为255字节，响应数据最长为256字节，动态签名数据作为IC卡响应数据，也限制了IC卡公钥的最大长度。

如卡片支持DDA和CDA，包含IC卡证书的记录模板的长度，即IC卡公钥证书长度加上证书（’9F46’）和记录模板（’70’）的Tag和Length不超过254字节，则IC卡公钥长度不超过247字节，因而发卡行公钥长度最大长度也不超过247字节。

根据发卡行应用数据长度不同，IC卡公钥最大长度在205到240之间。如果Generate AC命令响应包含其他可选数据，IC卡公钥最大长度还应减去这些数据的长度（包括Tag和Length）。如果卡片应用支持INTERNALAUTHENTICATION格式二，IC卡公钥最大长度还应减去7字节。

在选择公钥模长时，应该考虑到比较密钥的生命周期同预期的因数分解进程。

发卡行公钥指数和IC卡公钥指数的值由发卡行决定。认证中心，发卡行和IC卡公钥指数必须等于3或2¹⁶＋1。

标识本数字签名算法的公钥算法标识必须编码为十六进制‘01’。

2.RSA算法源码

<<RSA.h>>

// LargeInt.h
// 大整数类的声明，包含了质数检查，随机数生成，以及RSA密钥生成及加密解密函数。

//////////////////////////////////////////////////////////////////////////////

//RSA加密里所说的1024 位加密、2048位加密，是不是只要
//CLargeInt::CreateRandom(e,1024);
//这句代码里的参数设置成1024、2048就行了？
//
////不是的，2048bit，是指P*Q的积是2048bit，也就是说，当你P和Q都设置为1024的时候，加密强度是2048bit，事实上P和Q的位数可以不等，并且过分接近的P和Q将导致破解难度极大降低。
//
//
//_data 和 _len 这两个变量到底起什么作用？
//_data 占了4092个字节的大小，_len＝64时，只输出了512个字节。
////_data和_len其实是模拟了一个动态数组，用来容纳长度不定的大整数。_data 占据的大小是预先设置好的，这个数值在LargeInt.h中有定义：
//enum LIMIT{ LENGTH = 0x3FF,SMALLPRIMES = 5000};
//LENGTH 可以根据需要修改。
//_len是以4字节为单位的目前已经占用的数组的长度，_len = 64就是当前大整数的长度是64*4 = 256字节
//输出的字节是以16进制表示的字符串，每个字节需要两个符号来表示，所以长度是256*2 = 512字节。

//必须保证T_DWORD是4字节。
typedef unsigned long T_DWORD;

class CLargeInt
{
public:
	//说明：
	//基本的操作函数之所以没有设置成运算符的重载，主要是因为，
	//大数运算的过程中，经常需要额外的参数，例如错位相加减的offset值。
	//以及某些额外的结果
	//例如，除法运算的商和余数。
	//为了代码尽量通用，都写成了静态成员函数。
	//运算符的重载可以通过对这些静态成员函数进一步封装来实现。

	//复制
	static void Copy(const CLargeInt& source,CLargeInt& target);

	//比较
	static bool LargerEqual(const CLargeInt& source,const CLargeInt& target,T_DWORD offset = 0);
	static bool Equal(const CLargeInt& source,const CLargeInt& target);

	//和T_DWORD类型交互的底层操作函数
	static void Mul(const CLargeInt& faciend,T_DWORD multiplier,CLargeInt& product);
	static void Div(const CLargeInt& dividend,T_DWORD divisor,CLargeInt& quotient,T_DWORD& residual);

	//移位：
	//右移一位，用来取代 /2 操作。
	static void RCR(CLargeInt& n);
	//左移一位，用来取代 *2 操作。
	static void RCL(CLargeInt& n);

	//CLargeInt类型底层操作函数，
	static void Add(const CLargeInt& augend, const CLargeInt& addend, CLargeInt& sum,T_DWORD offset = 0);
	static void Sub(const CLargeInt& minuend,const CLargeInt& subtrahend, CLargeInt& difference,T_DWORD offset = 0);
	static void Mul(const CLargeInt& faciend, const CLargeInt& multiplier, CLargeInt& product);
	static void Div(const CLargeInt& dividend, const CLargeInt& divisor, CLargeInt& quotient,CLargeInt& residual);
	static void ExpMod(const CLargeInt& source, const CLargeInt& exponent, const CLargeInt& modulo,CLargeInt& result);

	//高阶的测试
	static bool RabinMillerTest(const CLargeInt& source, const CLargeInt& base);
	static bool Coprime(const CLargeInt &source, const CLargeInt &target);

	//RSA算法相关函数：
	static bool IsPrime(const CLargeInt &n);
	static void CreatePrime(CLargeInt &n);
	static bool RSACreate( const CLargeInt &p,const CLargeInt & q,const CLargeInt &e,CLargeInt &d,CLargeInt &n);
	static void RSAEncode( const CLargeInt &n,const CLargeInt &d,const CLargeInt &m,CLargeInt &c);
	static void RSADecode( const CLargeInt &n,const CLargeInt &e,const CLargeInt &c,CLargeInt &m);
	static string RSADecode(const char* N, const char* E, const char* I);
	//建立随机数。
	static void CreateRandom(CLargeInt &n,T_DWORD bitcount);
public:
	CLargeInt(T_DWORD value = 0);
	CLargeInt(const char* str);
	CLargeInt(const CLargeInt& other);

	~CLargeInt();

private:
	enum LIMIT{ LENGTH = 0x3FF, SMALLPRIMES = 5000};
	T_DWORD _len;
	mutable T_DWORD _data[LENGTH]; //之所以用 mutable 是因为需要在计算中获取指针，并且有时还需要在不改变数值的情况下修改某些数组元素。

	//辅助用的小质数表和相关的函数。
	static std::vector<T_DWORD> _smallPrimes;
	static void InitSmallPrimes(T_DWORD count = SMALLPRIMES);
	static bool SmallPrimeTest(const CLargeInt& n);
	string GetHexStr();
};

// LARGEINT_HEADER

<<RSA.CPP>>

std::vector<T_DWORD> CLargeInt::_smallPrimes;
void CLargeInt::InitSmallPrimes(T_DWORD	count)
{
	if (!_smallPrimes.size()) //尚未初始化
	{
		_smallPrimes.reserve(count);
		_smallPrimes.push_back(3);
		_smallPrimes.push_back(5);
		_smallPrimes.push_back(7);
		_smallPrimes.push_back(11);
		_smallPrimes.push_back(13);
		T_DWORD	i = 0;
		for( i = 17; _smallPrimes.size() < count; i+= 2)
		{
			bool bePrime = true;
			for( T_DWORD j = 0; j <	_smallPrimes.size(); j++ )
			{
				if( ( _smallPrimes[j] *	_smallPrimes[j]	) > i )	//已经不可能了。
					break;
				if( (i % _smallPrimes[j]) == 0)	//可以整除
					bePrime	= false;
			}
			if( bePrime ) //是素数
			{
				_smallPrimes.push_back(i);
			}
		}
	}
}

bool CLargeInt::SmallPrimeTest(const CLargeInt&	n)
{
	if (n._data[0] == 2 && n._len == 1)
	{
		return true; //2是质数。
	}

	if (n._len == 1	&& n._data[0] <= _smallPrimes.back()) //整数比表中的数据小。
	{
		if (std::binary_search(_smallPrimes.begin(),_smallPrimes.end(),n._data[0])) //找到，n存在于质数表中。
		{
			return true; //通过测试，数n是质数。
		}

		return false; //没有通过，n是合数。
	}
	else
	{
		CLargeInt temp;
		T_DWORD	r;
		for( T_DWORD i = 0; i <	_smallPrimes.size(); i++)
		{
			Div(n,_smallPrimes[i],temp,r); //依次检查是否能整除。
			if( r == 0) return false; //能整除，说明 n 有小于它的质因子。n 是合数。
		}

		return true; //	n 没有表中已有的质因子，n有可能是质数。
	}
}

bool CLargeInt::IsPrime(const CLargeInt	&n)
{
	InitSmallPrimes();
	if (n._data[0]&1) //奇数
	{
		if (n._len == 1	&& n._data[0] <= _smallPrimes.back())
		{
			return SmallPrimeTest(n);  //是小质数表中的一个质数。
		}
		else
		{
			if( SmallPrimeTest(n) )	//通过小质数测试
			{
				for( int i = 0;	i < 5; i++)
				{
					CLargeInt a(_smallPrimes[ _smallPrimes.size() -	i]);
					if( !RabinMillerTest(n,a) ) return false;
				}

				return true; //99.9%的可能是质数。
			}
			else
			{
				return false;
			}
		}
	}

	return false;
}

void CLargeInt::CreatePrime(CLargeInt &n)
{
	n._data[0] |= 0x01;
	while (!IsPrime(n))
	{
		Add(n, 2, n); //步进2，寻找最接近的下一个质数。这个算法并不是寻找质数最好的算法，但是它可以保证寻找到相邻的质数。
	}
}

CLargeInt::CLargeInt(T_DWORD value) : _len(1)
{
	_data[0] = value;
}

CLargeInt::CLargeInt(const CLargeInt& other)
{
	Copy(other,*this);
}

CLargeInt::CLargeInt(const char	* str):	_len(1)
{
	_data[0] = 0;

	//if(str && strlen(str) >	2 && ( strncmp(str,"0x",2) == 0	|| strncmp(str,"0X",2) == 0 ) )
	if (str)
	{
		CLargeInt temp;
		char ch;
		T_DWORD	n = 0;
		const char * p = str;
		while (*p)
		{
			ch	= *p;
			if ((ch >= '0') && (ch <= '9'))
			{
				n = ch - '0';
			}
			else if ((ch >= 'a') && (ch <= 'f'))
			{
				n = ch - 'a' + 10;       //a->f分别代表10->15   a的AscII为97，所以这里应该减去87
			}
			else if ((ch >= 'A') && (ch <= 'F'))
			{
				n = ch - 'A' + 10;       //同上
			}
			p++;
			Mul(*this, 16, temp);
			Add(temp, n, *this);
		}
	}
	//else
	//{//这里应该用来处理除了0x开头的16进制字符串的情况。
	//	__asm int 3;
	//}
}

CLargeInt::~CLargeInt()
{

}

void CLargeInt::Copy(const CLargeInt& source, CLargeInt& target)
{
	if (&source != &target) //避免自身拷贝。
	{
		target._len = source._len;

		//相信库函数能提供最高效的实现。
		memcpy(target._data, source._data, source._len * 4);
	}
}

bool CLargeInt::Equal(const CLargeInt& source, const CLargeInt& target)
{
	if (source._len	== target._len)
	{
		//相信库函数能提供最高效的实现。
		return !(memcmp(source._data, target._data,source._len * 4));
	}
	else
	{
		return false;
	}
}

bool CLargeInt::LargerEqual(const CLargeInt& source,const CLargeInt& target,T_DWORD offset )
{
	if( source._len	> (target._len + offset) ) return true;
	else if( source._len < (target._len + offset) )	return false;
	else
	{
		int end	= offset;
		for( int i = source._len - 1; i	>= end ; i--)
		{
			T_DWORD	ii = i - offset;
			if( source._data[i] > target._data[ii] ) return	true; //大于
			else if( source._data[i] < target._data[ii] )	return false;
		}
		return true; //相等。
	}

}

void CLargeInt::Add(const CLargeInt &augend, const CLargeInt &addend, CLargeInt	&sum,T_DWORD offset)
{
	T_DWORD	len,len1,len2; //len绝对不能<1，否则会出错。
	T_DWORD	*p1 = 0;
	T_DWORD	*p2 = 0;
	T_DWORD	*p3 = 0;

	if( augend._len	>= (addend._len+offset)	)
	{
		len = augend._len - offset;
		len1 = addend._len;
		len2 =	len - len1;
		p1 = augend._data+offset;
		p2 = addend._data;
		p3 = sum._data+offset;
	}
	else
	{
		len = addend._len;
		if( augend._len	<= offset )
		{
			augend._data[offset] = 0;
			len1 = 1;
			len2 = addend._len - len1;
		}
		else
		{
			len1 = augend._len - offset;
			len2 = addend._len - len1;
		}
		p1 = addend._data;
		p2 = augend._data + offset;
		p3 = sum._data + offset;
	}

	__asm
	{
		pushad;
		//载入计算地址
		mov esi,p1;
		mov ebx,p2;
		mov edi,p3;

		mov eax,len1;
		sal eax,2;

		//计算末地址
		add esi,eax;
		add ebx,eax;
		add edi,eax;

		mov ecx,eax;
		add ecx,4*7;
		and ecx,0xFFFFFFE0; //ecx - ecx%32,这里的ecx就变成了偏移址，同时也是累加器。
		neg ecx; //取负，用作累加。

		sar eax,2; //回复到原来的长度

		neg eax; //取负
		and eax,dword ptr 7;
		jz labeladd0; //如果低位是0，说明长度至少有8(因为长度为0是非法的)

		mov edx,dword ptr 0xC;
		mul edx; //eax变为相对偏移址。

		lea edx,labeladd1; //载入起始地址
		sub edx,dword ptr 0xC;
		add eax,edx; //计算出跳转地址

		clc; //清除标志位。
		jmp eax; //跳转
labeladdloop:
		popf;
labeladd0:
		mov eax,[esi+ecx+0*4];
		adc eax,[ebx+ecx+0*4];
		mov [edi+ecx+0*4],eax;
labeladd1:
		mov eax,[esi+ecx+1*4];
		adc eax,[ebx+ecx+1*4];
		mov [edi+ecx+1*4],eax;
		//labeladd2:
		mov eax,[esi+ecx+2*4];
		adc eax,[ebx+ecx+2*4];
		mov [edi+ecx+2*4],eax;
		//labeladd3:
		mov eax,[esi+ecx+3*4];
		adc eax,[ebx+ecx+3*4];
		mov [edi+ecx+3*4],eax;
		//labeladd4:
		mov eax,[esi+ecx+4*4];
		adc eax,[ebx+ecx+4*4];
		mov [edi+ecx+4*4],eax;
		//labeladd5:
		mov eax,[esi+ecx+5*4];
		adc eax,[ebx+ecx+5*4];
		mov [edi+ecx+5*4],eax;
		//labeladd6:
		mov eax,[esi+ecx+6*4];
		adc eax,[ebx+ecx+6*4];
		mov [edi+ecx+6*4],eax;
		//labeladd7:
		mov eax,[esi+ecx+7*4];
		adc eax,[ebx+ecx+7*4];
		mov [edi+ecx+7*4],eax;
		//labeladd8:
		pushf;
		add ecx,32;
		jnz labeladdloop;

		//运行到这里了，说明两数重合部分已经处理完了。开始处理两数相差的部分。

		mov ecx,len2; //载入长度

		xor ebx,ebx; //因为ebx指向的缓冲区废弃不用了，所以使用已经空闲的ebx寄存器。
		cmp ebx,ecx; //判断ecx是否等于ebx(0)
		jz labelcarry; //如果相等，说明已经处理完了，跳转到处理最后一次进位的地方。

labelfix:
		popf; //弹出上次的标志位。

		mov eax,[esi+ebx*4];
		adc eax,0;
		mov [edi+ebx*4],eax;

		//如果这里没有进位，剩下的部分可以直接拷贝。拷贝的长度就是ecx-1
		jnc labelcopy;

		pushf; //保存标志位。

		inc ebx; //步进
		dec ecx;
		jnz labelfix; //循环。

labelcarry:
		popf; //弹出标志位
		jnc labelend; //没有进位就直接结束
		//没有跳转说明有进位：
		mov dword ptr [edi+ebx*4],dword	ptr 0x00000001;	//直接置1
		lea eax,len;
		inc [eax]; //总长度+1
		inc ecx; //补一次。
labelcopy:
		dec ecx;
		sal ecx,2;
		cmp edi,esi; //比较源地址和目标地址，如果相同就跳过。
		jz labelend;

		sal ebx,2;
		add ebx,4;

		add edi,ebx;
		add esi,ebx;

		rep movs dword ptr [edi],dword ptr [esi];

labelend:
		popad;
	}
	sum._len = len + offset;
}

void CLargeInt::Sub(const CLargeInt& minuend,const CLargeInt& subtrahend,CLargeInt& difference,T_DWORD offset)
{
	T_DWORD	len,len1,len2; //len绝对不能<1，否则会出错。
	T_DWORD	*p1 = 0;
	T_DWORD	*p2 = 0;
	T_DWORD	*p3 = 0;

	if( minuend._len >= (subtrahend._len+offset) )
	{
		len = minuend._len- offset;
		len1 = subtrahend._len;
		len2 =	len - len1;
		p1 = minuend._data + offset;
		p2 = subtrahend._data;
		p3 = difference._data +	offset;
	}
	else
	{
		__asm int 3; //出错。
	}

	__asm
	{
		pushad;
		//载入计算地址
		mov esi,p1;
		mov ebx,p2;
		mov edi,p3;

		mov eax,len1;
		sal eax,2;

		//计算末地址
		add esi,eax;
		add ebx,eax;
		add edi,eax;

		mov ecx,eax;
		add ecx,4*7;
		and ecx,0xFFFFFFE0; //ecx - ecx%32,这里的ecx就变成了偏移址，同时也是累加器。
		neg ecx; //取负，用作累加。

		sar eax,2; //回复到原来的长度

		neg eax; //取负
		and eax,dword ptr 7;
		jz labeladd0; //如果低位是0，说明长度至少有8(因为长度为0是非法的)

		mov edx,dword ptr 0xC;
		mul edx; //eax变为相对偏移址。

		lea edx,labeladd1; //载入起始地址
		sub edx,dword ptr 0xC;
		add eax,edx; //计算出跳转地址

		clc; //清除标志位。
		jmp eax; //跳转
labeladdloop:
		popf;
labeladd0:
		mov eax,[esi+ecx+0*4];
		sbb eax,[ebx+ecx+0*4];
		mov [edi+ecx+0*4],eax;
labeladd1:
		mov eax,[esi+ecx+1*4];
		sbb eax,[ebx+ecx+1*4];
		mov [edi+ecx+1*4],eax;
		//labeladd2:
		mov eax,[esi+ecx+2*4];
		sbb eax,[ebx+ecx+2*4];
		mov [edi+ecx+2*4],eax;
		//labeladd3:
		mov eax,[esi+ecx+3*4];
		sbb eax,[ebx+ecx+3*4];
		mov [edi+ecx+3*4],eax;
		//labeladd4:
		mov eax,[esi+ecx+4*4];
		sbb eax,[ebx+ecx+4*4];
		mov [edi+ecx+4*4],eax;
		//labeladd5:
		mov eax,[esi+ecx+5*4];
		sbb eax,[ebx+ecx+5*4];
		mov [edi+ecx+5*4],eax;
		//labeladd6:
		mov eax,[esi+ecx+6*4];
		sbb eax,[ebx+ecx+6*4];
		mov [edi+ecx+6*4],eax;
		//labeladd7:
		mov eax,[esi+ecx+7*4];
		sbb eax,[ebx+ecx+7*4];
		mov [edi+ecx+7*4],eax;
		//labeladd8:
		pushf;
		add ecx,32;
		jnz labeladdloop;

		//运行到这里了，说明两数重合部分已经处理完了。开始处理两数相差的部分。

		mov ecx,len2; //载入长度,
		//只处理p1的buf

		xor ebx,ebx; //因为ebx指向的缓冲区废弃不用了，所以使用已经空闲的ebx寄存器。
		cmp ebx,ecx; //判断ecx是否等于ebx
		jz labelcarry; //如果相等，说明已经处理完了，跳转到处理最后一次进位的地方。

labelfix:
		popf; //弹出上次的标志位。

		mov eax,[esi+ebx*4];
		sbb eax,0;
		mov [edi+ebx*4],eax;

		//如果这里没有进位，剩下的部分可以直接拷贝。拷贝的长度就是ecx-1
		jnc labelcopy;

		pushf; //保存标志位。

		inc ebx; //步进
		dec ecx;
		jnz labelfix; //循环。

labelcarry:
		popf; //弹出标志位
		jnc labelend; //没有进位就直接结束
		//没有跳转说明有进位：
		int 3; //说明减数不够，结果为负
		mov dword ptr [edi+ebx*4],dword	ptr 0x00000001;	//直接置1
		lea eax,len;
		dec [eax]; //总长度-1
		inc ecx; //补一次。
labelcopy:
		dec ecx;
		sal ecx,2;
		cmp edi,esi; //比较源地址和目标地址，如果相同就跳过。
		jz labelend;

		sal ebx,2;
		add ebx,4;

		add edi,ebx;
		add esi,ebx;

		rep movs dword ptr [edi],dword ptr [esi];

labelend:
		popad;
	}
	difference._len	= len+offset;
	for( T_DWORD i = difference._len-1; i >	0; i--)
	{
		if( difference._data[i]	== 0 ) //末尾是0
		{
			difference._len--;
		}
		else break;
	}
}

void CLargeInt::Mul(const CLargeInt& faciend,T_DWORD multiplier,CLargeInt& product)
{
	T_DWORD	len; //len绝对不能<1，否则会出错。
	T_DWORD	*p1 = 0;
	T_DWORD	*p2 = 0;
	T_DWORD	*p3 = 0;

	len = faciend._len;
	p1 = faciend._data;
	p2 = 0;
	p3 = product._data;

	memset(p3,0,(len+1)*4);
	__asm
	{
		pushad;
		//载入计算地址
		mov esi,p1;
		mov edi,p3;
		mov ebx,multiplier

			mov eax,len;
		sal eax,2;

		//计算末地址
		add esi,eax;
		add edi,eax;

		mov ecx,eax;
		add ecx,4*7;
		and ecx,0xFFFFFFE0; //ecx - ecx%32,这里的ecx就变成了偏移址，同时也是累加器。
		neg ecx; //取负，用作累加。

		sar eax,2; //回复到原来的长度

		neg eax; //取负
		and eax,dword ptr 7;
		jz labeladd0; //如果低位是0，说明长度至少有8(因为长度为0是非法的)

		mov edx,dword ptr 0xe; //0xe是跳转地址的系数。
		mul edx; //eax变为相对偏移址。

		lea edx,labeladd1; //载入起始地址
		sub edx,dword ptr 0xe;
		add eax,edx; //计算出跳转地址

		clc; //清除标志位。
		jmp eax; //跳转
labeladdloop:
labeladd0:
		mov eax,[esi+ecx+0*4];
		mul ebx;
		add [edi+ecx+0*4],eax;
		adc [edi+ecx+1*4],edx;
labeladd1:
		mov eax,[esi+ecx+1*4];
		mul ebx;
		add [edi+ecx+1*4],eax;
		adc [edi+ecx+2*4],edx;
		//labeladd2:
		mov eax,[esi+ecx+2*4];
		mul ebx;
		add [edi+ecx+2*4],eax;
		adc [edi+ecx+3*4],edx;
		//labeladd3:
		mov eax,[esi+ecx+3*4];
		mul ebx;
		add [edi+ecx+3*4],eax;
		adc [edi+ecx+4*4],edx;
		//labeladd4:
		mov eax,[esi+ecx+4*4];
		mul ebx;
		add [edi+ecx+4*4],eax;
		adc [edi+ecx+5*4],edx;
		//labeladd5:
		mov eax,[esi+ecx+5*4];
		mul ebx;
		add [edi+ecx+5*4],eax;
		adc [edi+ecx+6*4],edx;
		//labeladd6:
		mov eax,[esi+ecx+6*4];
		mul ebx;
		add [edi+ecx+6*4],eax;
		adc [edi+ecx+7*4],edx;
		//labeladd7:
		mov eax,[esi+ecx+7*4];
		mul ebx;
		add [edi+ecx+7*4],eax;
		adc [edi+ecx+8*4],edx;
		//labeladd8:
		add ecx,32;
		jnz labeladdloop;

		//labelend:
		popad;
	}
	product._len = len+1;

	for( T_DWORD i = product._len-1; i > 0;	i--)
	{
		if( product._data[i] ==	0 ) //末尾是0
		{
			product._len--;
		}
		else break;
	}

}

void CLargeInt::RCR(CLargeInt& n)
{
	T_DWORD	len; //len绝对不能<1，否则会出错。
	T_DWORD	*p1 = 0;
	len = n._len;
	p1 = n._data;

	__asm
	{
		pushad;		

		//载入计算地址
		mov esi,p1;

		mov eax,len;
		mov ecx,eax;

		add ecx,dword ptr 15;
		and ecx,dword ptr 0xFFFFFFF0;

		mov eax,ecx;
		sal ecx,dword ptr 2;

		sub eax,len; //求初始偏移。

		mov ebx,dword ptr 4; //偏移倍率
		mul ebx;		

		lea ebx,label0;
		add eax,ebx;
		clc;
		jmp eax;
labelloop:
		popf;
label0:
		rcr dword ptr [esi + ecx - 1*4],1;
		//label1:
		rcr dword ptr [esi + ecx - 2*4],1;
		//label2:
		rcr dword ptr [esi + ecx - 3*4],1;
		//label3:
		rcr dword ptr [esi + ecx - 4*4],1;
		//label4:
		rcr dword ptr [esi + ecx - 5*4],1;
		//label5:
		rcr dword ptr [esi + ecx - 6*4],1;
		//label6:
		rcr dword ptr [esi + ecx - 7*4],1;
		//label7:
		rcr dword ptr [esi + ecx - 8*4],1;
		//label8:
		rcr dword ptr [esi + ecx - 9*4],1;
		//label9:
		rcr dword ptr [esi + ecx - 10*4],1;
		//label10:
		rcr dword ptr [esi + ecx - 11*4],1;
		//label11:
		rcr dword ptr [esi + ecx - 12*4],1;
		//label12:
		rcr dword ptr [esi + ecx - 13*4],1;
		//label13:
		rcr dword ptr [esi + ecx - 14*4],1;
		//label14:
		rcr dword ptr [esi + ecx - 15*4],1;
		//label15:
		rcr dword ptr [esi + ecx - 16*4],1;
		//label16:
		pushf;
		sub ecx,dword ptr 64;
		jnz labelloop;
		popf;
		popad;
	}
	for( T_DWORD i = n._len-1; i > 0; i--)
	{
		if( n._data[i] == 0 ) //末尾是0
		{
			n._len--;
		}
		else break;
	}
}

void CLargeInt::RCL(CLargeInt& n)
{
	T_DWORD	len = n._len; //len绝对不能<1，否则会出错。
	T_DWORD	*p1 = 0;

	n._data[len] = 0;
	n._len++;

	len++;
	p1 = n._data;

	__asm
	{
		pushad;		

		//载入计算地址
		mov esi,p1;

		mov eax,len;
		mov ecx,eax;

		sal eax,dword ptr 2;
		add esi,eax;

		add ecx,dword ptr 15;
		and ecx,dword ptr 0xFFFFFFF0;

		mov eax,ecx;
		sal ecx,dword ptr 2;

		neg ecx;

		sub eax,len; //求初始偏移。
		jz label0;

		mov ebx,dword ptr 4; //偏移倍率
		mul ebx;		

		lea ebx,label1;
		sub ebx,4;
		add eax,ebx;
		clc;
		jmp eax;
labelloop:
		popf;
label0:
		rcl dword ptr [esi + ecx + 0*4],1;
label1:
		rcl dword ptr [esi + ecx + 1*4],1;
		//label2:
		rcl dword ptr [esi + ecx + 2*4],1;
		//label3:
		rcl dword ptr [esi + ecx + 3*4],1;
		//label4:
		rcl dword ptr [esi + ecx + 4*4],1;
		//label5:
		rcl dword ptr [esi + ecx + 5*4],1;
		//label6:
		rcl dword ptr [esi + ecx + 6*4],1;
		//label7:
		rcl dword ptr [esi + ecx + 7*4],1;
		//label8:
		rcl dword ptr [esi + ecx + 8*4],1;
		//label9:
		rcl dword ptr [esi + ecx + 9*4],1;
		//label10:
		rcl dword ptr [esi + ecx + 10*4],1;
		//label11:
		rcl dword ptr [esi + ecx + 11*4],1;
		//label12:
		rcl dword ptr [esi + ecx + 12*4],1;
		//label13:
		rcl dword ptr [esi + ecx + 13*4],1;
		//label14:
		rcl dword ptr [esi + ecx + 14*4],1;
		//label15:
		rcl dword ptr [esi + ecx + 15*4],1;
		//label16:
		pushf;
		add ecx,dword ptr 64;
		jnz labelloop;
		popf;
		popad;
	}
	for( T_DWORD i = n._len-1; i > 0; i--)
	{
		if( n._data[i] == 0 ) //末尾是0
		{
			n._len--;
		}
		else break;
	}
}

void CLargeInt::Mul(const CLargeInt& faciend, const CLargeInt& multiplier, CLargeInt& product)
{
	CLargeInt temp;
	product._len = (faciend._len + multiplier._len);
	memset(product._data, 0, product._len * 4);

	for (T_DWORD i = 0; i <	multiplier._len; i++)
	{
		Mul(faciend,multiplier._data[i],temp);
		Add(product,temp,product,i);
	}
	for (T_DWORD i = product._len - 1; i > 0;	i--)
	{
		if( product._data[i] ==	0 ) //末尾是0
		{
			product._len--;
		}
		else
		{
			break;
		}
	}
}

void CLargeInt::Div(const CLargeInt& dividend,T_DWORD divisor,CLargeInt& quotient,T_DWORD& residual)
{
	T_DWORD	len; //len绝对不能<1，否则会出错。
	T_DWORD	*p1 = 0;
	T_DWORD	*p2 = 0;

	T_DWORD	*pr = &residual;

	len = dividend._len;
	quotient._len =	dividend._len;
	p1 = dividend._data;
	p2 = quotient._data;

	__asm
	{
		pushad;		

		//载入计算地址
		mov esi,p1;
		mov edi,p2;

		mov eax,len;
		mov ecx,eax;

		add ecx,dword ptr 7;
		and ecx,dword ptr 0xFFFFFFF8;

		mov eax,ecx;
		sal ecx,dword ptr 2;

		sub eax,len; //求初始偏移。

		mov ebx,dword ptr 0x0A;	//偏移倍率
		mul ebx;		

		lea edx,label0;
		add eax,edx;

		xor edx,edx;

		mov ebx,divisor;

		clc;
		jmp eax;
labelloop:

label0:
		mov eax,[esi + ecx - 1*4];
		div ebx;
		mov [edi + ecx - 1*4],eax;
		//label1:
		mov eax,[esi + ecx - 2*4];
		div ebx;
		mov [edi + ecx - 2*4],eax;
		//label2:
		mov eax,[esi + ecx - 3*4];
		div ebx;
		mov [edi + ecx - 3*4],eax;
		//label3:
		mov eax,[esi + ecx - 4*4];
		div ebx;
		mov [edi + ecx - 4*4],eax;
		//label4:
		mov eax,[esi + ecx - 5*4];
		div ebx;
		mov [edi + ecx - 5*4],eax;
		//label5:
		mov eax,[esi + ecx - 6*4];
		div ebx;
		mov [edi + ecx - 6*4],eax;
		//label6:
		mov eax,[esi + ecx - 7*4];
		div ebx;
		mov [edi + ecx - 7*4],eax;
		//label7:
		mov eax,[esi + ecx - 8*4];
		div ebx;
		mov [edi + ecx - 8*4],eax;
		//label8:
		sub ecx,dword ptr 32;
		jnz labelloop;

		//运行到这里说明已经计算完毕，最后保存余数。
		mov ecx,pr;
		mov [ecx],edx;
		popad;
	}

	for (T_DWORD i = quotient._len-1; i > 0; i--)
	{
		if (quotient._data[i] == 0) //末尾是0
		{
			quotient._len--;
		}
		else
		{
			break;
		}
	}
}

void CLargeInt::Div(const CLargeInt& dividend,const CLargeInt& divisor,CLargeInt& quotient,CLargeInt& residual)
{
	if( dividend._len < divisor._len )
	{
		quotient._len =	1;
		quotient._data[0] = 0;
		Copy(dividend,residual);
	}
	else
	{
		if( divisor._len == 1)
		{
			Div(dividend,divisor._data[0],quotient,residual._data[0]);
			residual._len =	1;
		}
		else
		{
			CLargeInt d,r;

			//满位算法
			T_DWORD	head = divisor._data[divisor._len-1];
			__asm
			{
				pushad;
				mov eax,head;
				bsr ecx,eax;
				mov edx,dword ptr 31;
				sub edx,ecx;
				mov ecx,edx;
				mov eax,dword ptr 0x01;
				sal eax,cl;
				mov head,eax;
				popad;
			}

			Mul(dividend,head,d);
			Mul(divisor,head,r);

			quotient._len =	d._len - r._len	+ 1;
			quotient._data[quotient._len - 1] = 0;

			T_DWORD	newhead	= r._data[r._len - 1];
			//以上处理的目的是使得被除数最高位的1移动到DWORD的最高位。
			//处理之后可以极大减少商的上下限之间的差距。

			T_DWORD	highpart = 0;
			T_DWORD	lowpart	 = 0;
			T_DWORD	max = 0,min = 0;
			T_DWORD	offset = 0;
			CLargeInt temp;
			for( T_DWORD i = d._len-1; i >=	r._len-1; i-- )
			{
				offset = i - (r._len - 1);
				lowpart	= d._data[i];

				__asm
				{//这段汇编代码的作用是求商的上下限。min & max
					pushad;
					mov edx,highpart;
					mov eax,lowpart;
					mov ecx,newhead;
					div ecx;
					mov max,eax;

					inc ecx;
					jz _mov;
					mov edx,highpart;
					mov eax,lowpart;
					div ecx;
					mov edx,eax;
_mov:
					mov min,edx;
					popad;
				}

				if (max == 0)
				{
					highpart = d._data[i];
					quotient._data[offset] = 0;
				}
				else
				{
					Mul(r,max,temp);
					if ( max	!= min )
					{
						while ( !LargerEqual(d,temp,offset) )
						{//这里使用的其实是试商法，
							//由于max和min的差距已经很小，所以这里不再折半试商。
							max--;
							Sub(temp,r,temp);
						}
					}	

					quotient._data[offset] = max; //保存商。
					Sub(d,temp,d,offset); //求本阶的余数。
					highpart = d._data[i]; //载入余数。
				}
			}

			//清除高位的0。
			for(T_DWORD i = quotient._len-1; i > 0; i--)
			{
				if( quotient._data[i] == 0 ) //末尾是0
				{
					quotient._len--;
				}
				else break;
			}
			//除法完成，计算余数：
			Div(d,head,residual,newhead);

			if ( newhead != 0 )
			{
				//余数不为0，说明出错。
				__asm int 3;
			}
		}

	}
}

void CLargeInt::ExpMod(const CLargeInt&	source,const CLargeInt&	exponent,const CLargeInt& modulo,CLargeInt& result)
{
	CLargeInt n,e,r(1),temp,temp1;
	Copy(source,n);
	Copy(exponent,e);

	while( e._len >	1 || e._data[e._len - 1] > 1 )
	{
		if( e._data[0] &1 )
		{
			Mul(r,n,temp);
			Div(temp,modulo,temp1,r);
		}
		RCR(e);
		Mul(n,n,temp);
		Div(temp,modulo,temp1,n);
	}
	Mul(r,n,temp);
	Div(temp,modulo,temp1,result);
}

bool CLargeInt::RabinMillerTest(const CLargeInt& source,const CLargeInt& base)
{
	CLargeInt m,temp(1);

	Sub(source,temp,m);

	T_DWORD	count =	0;
	while(!(m._data[0] & 0x01) )
	{
		count++;
		RCR(m);
	}
	/*
	这里注释掉的是旧的算法，可以用后面的等效算法取代。
	改变的算法使得计算ExpMod的次数由 count 减少到 1
	虽然理论上计算次数大大减少了，但是实际效果似乎只是减少了50%左右的运算量。

	for( T_DWORD i = 0; i< count;	i++)
	{
	CLargeInt mod;
	ExpMod(base,m,source,mod);
	if( mod._data[0] == 1	&& mod._len == 1 )
	{
	return true;
	}
	else
	{
	Sub(source,mod,temp);
	if(  temp._data[0] ==	1 && temp._len == 1)
	{
	return true;
	}
	}
	RCL(m);
	}
	*/

	CLargeInt mod;
	ExpMod(base,m,source,mod);
	if( mod._data[0] == 1 && mod._len == 1 )
	{
		return true;
	}
	else
	{
		Sub(source,mod,temp);
		if(  temp._data[0] == 1	&& temp._len ==	1)
		{
			return true;
		}
	}
	for( T_DWORD i = 1; i< count; i++)
	{
		CLargeInt next,t;
		Mul(mod,mod,next);
		Div(next,source,t,mod);
		if( mod._data[0] == 1 && mod._len == 1 )
		{
			return true;
		}
		else
		{
			Sub(source,mod,temp);
			if(  temp._data[0] == 1	&& temp._len ==	1)
			{
				return true;
			}
		}
	}

	return false;
}

bool CLargeInt::RSACreate( const CLargeInt &p,const CLargeInt &	q,const	CLargeInt &e,CLargeInt &d,CLargeInt &n)
{//由已知的P,Q,E计算N,D，完成RSA密钥生成。
	Mul(p,q,n);
	CLargeInt a(e),b(n);
	Sub(b,p,b);
	Sub(b,q,b);
	Add(b,1,b);

	if( !Coprime(b,e) ) return false;
	//求D的算法，通过辗转相除，最终求得需要的D值。
	CLargeInt k1,c1;
	Div(b,a,k1,c1);
	if( c1._len == 1 && c1._data[0]	== 0 ) return false;
	if( c1._len == 1 && c1._data[0]	== 1 )
	{
		CLargeInt temp;
		Sub(e,1,temp);
		Mul(temp,k1,d);
		Add(d,1,d);
		return true;
	}

	CLargeInt k2,c2,ka2;
	Div(a,c1,k2,c2);
	if( c2._len == 1 && c2._data[0]	== 0 ) return false;
	Mul(k1,k2,ka2);
	Add(ka2,1,ka2);
	if( c2._len == 1 && c2._data[0]	== 1 ) { Copy(ka2,d); return true; }

	CLargeInt k3,c3,ka3;
	Div(c1,c2,k3,c3);
	if( c3._len == 1 && c3._data[0]	== 0 ) return false;
	Mul(k3,ka2,ka3);
	Add(ka3,k1,ka3);
	if( c3._len == 1 && c3._data[0]	== 1 )
	{
		Copy(ka3,d);

		CLargeInt temp,temp1,temp2;
		Mul(d,e,temp);
		Div(temp,b,temp1,temp2);
		Add(temp2,1,temp2);
		if( Equal(temp2,b) )
		{
			Mul(d,d,temp1);
			Div(temp1,b,temp2,d);
			Mul(d,e,temp1);
			Div(temp1,b,temp2,d);
		}
		return true;
	}

	CLargeInt kn_2(k2),cn_2(c2),ka_2(ka2);
	CLargeInt kn_1(k3),cn_1(c3),ka_1(ka3);

	CLargeInt kn,cn,kan;
	do{
		Div(cn_2,cn_1,kn,cn);
		Mul(kn,ka_1,kan);
		Add(kan,ka_2,kan);

		if( cn._len == 1 && cn._data[0]	== 0 ) return false;
		if( cn._len == 1 && cn._data[0]	== 1 )
		{
			Copy(kan,d);

			CLargeInt temp,temp1,temp2;
			Mul(d,e,temp);
			Div(temp,b,temp1,temp2);
			Add(temp2,1,temp2);
			if( Equal(temp2,b) )
			{
				Mul(d,d,temp1);
				Div(temp1,b,temp2,d);
				Mul(d,e,temp1);
				Div(temp1,b,temp2,d);
			}
			return true;
		}
		Copy(kn_1,kn_2);
		Copy(cn_1,cn_2);
		Copy(ka_1,ka_2);

		Copy(kn,kn_1);
		Copy(cn,cn_1);
		Copy(kan,ka_1);
	} while (true);
}

void CLargeInt::RSAEncode( const CLargeInt &n,const CLargeInt &d,const CLargeInt &m,CLargeInt &c)
{
	ExpMod(m,d,n,c);
}

void CLargeInt::RSADecode( const CLargeInt &n,const CLargeInt &e,const CLargeInt &c,CLargeInt &m)
{
	ExpMod(c,e,n,m);
}

string CLargeInt::RSADecode(const char* N, const char* E, const char* I)
{
	CLargeInt o;
	CLargeInt n(N);
	CLargeInt i(I);
	CLargeInt e(E);
	RSADecode(n, e, i, o);
	return o.GetHexStr();
}

string CLargeInt::GetHexStr()
{
	char buffer[128];
	memset(buffer, 0x00, sizeof(buffer));
	string strHex;

	for (int i = _len - 1 ; i >= 0; i--)
	{
		char str[] = "%08X";
		str[2] = '0' + sizeof(T_DWORD) * 2;
		if ( i != (int)(_len - 1) )
		{
			sprintf(buffer, str, _data[i]);
		}
		else
		{
			sprintf(buffer, str, _data[i]);
		}

		strHex += buffer;
	}
	return strHex;
}

typedef	unsigned __int64 ULONGLONG;
inline ULONGLONG GetCycleCount()
{
	__asm RDTSC
}

//#include <windows.h>
void CLargeInt::CreateRandom(CLargeInt &n, T_DWORD bitcount)
{
	n._len = (bitcount + 31) / 32;
	for (T_DWORD i = 0; i <	n._len;	i++)
	{
		T_DWORD	byte[4];

		for (int b = 0;	b < 4; b++)
		{
			for (int j = 0;	j < 4;	j++)
			{
				//Sleep(0);
				j = j + 1 - 1;
			}
			byte[b]	= GetCycleCount() & 0xFF;
		}

		n._data[i] = (byte[0] << 24)
			| (byte[1] << 16)
			| (byte[2] << 8)
			| (byte[3] << 0);
	}
	n._data[0] |= 1;
	n._data[n._len - 1] |= 0x80000000;
	if ((bitcount &	31))
	{
		n._data[n._len - 1] >>=	32 - (bitcount & 31);
	}
}

bool CLargeInt::Coprime(const CLargeInt	&source,const CLargeInt	&target)
{
	CLargeInt m(source),n(target),temp;
	while(true)
	{
		Div(m,n,temp,m);
		if( m._len == 1	&& m._data[0] == 1 ) return true;
		else if( m._len	== 1 &&	m._data[0] == 0	) return false;

		Div(n,m,temp,n);
		if( n._len == 1	&& n._data[0] == 1 ) return true;
		else if( n._len	== 1 &&	n._data[0] == 0	) return false;
	}

}

3.RSA工具

算法系列2——RSA,布布扣,bubuko.com