RSA大数运算实现（1024位n）（3）

  在(1)的基础上，采用Knuth提供的估商法来实现除法，会使得程序运行速度大幅加快，实际上整个程序的运行时间主要取决于除法的质量，使用Knuth大神的方法是绝佳选择。大神不愧是大神，方法tql！
  因为公式编辑不太方便，所以引用《计算机程序设计艺术·第2卷》中的一些图片。
  后面实现了另一种比较快的求逆算法，以及求贝祖等式和蒙哥马利算法之后再次更新。

  首先是对于被除数和除数的说明:


  重点是，用qhat(有帽子符号的q)去估计q，通过下面Knuth提出的定理，可以总结出来，|qhat-q|≤2，误差只会是2，而且满足条件的时候，一定是qhat-2≤q≤qhat。

实现算法

  在《密码学C/C++实现》中给出了一种可行的算法（算法部分中文版中有印刷错误，原理部分英文版和中文版均有错误，所以给的是Knuth原书中的部分）

实现细节

① uint32_t 的数移动32位的结果是不动不是0，如果要移位32位，哪怕是uint32_t的也要先换成uint64_t进行操作。
② 超过32位的两个数相乘结果会溢出，超过64位。
③ 做检验时，括号中不会超过64位，但是括号中结果乘以B就不知道了，可以做检验，检验括号中的数是否比B(232)大，如果大于B，不等式右边就比BB大；不等式左边bn-2q’小于BB，不用做不等式比较，跳过了溢出这种可能性；重复做的时候也是需要这样考虑。
④ 第5步判断，可以先执行r-qb，如果最后有借位，说明商估计大了，需要把一个b加回去，q-1是真实的商。

源代码

int div_b(BN a, BN b, BN q, BN rem)
{

	memset(rem, 0, sizeof(rem));
	memset(q, 0, sizeof(q));

	BND b_t;//临时存放b
	BND q_t = { 0 };//每个商
	uint32_t r_t[2 + (BNMAXDGT << 1)];

	uint32_t *bptr, *mbptr, *rptr, *rcir, *mrptr, *qptr;
	uint32_t d = 0, qh = 0,qh1=0;
	uint64_t kuo, left, right, borrow, carry;
	uint32_t bn_1 = 0, bn_2 = 0;
	uint32_t ri = 0, ri_1 = 0, ri_2 = 0;//三个r'，乘了d以后的
	cpy_b(r_t, a);
	cpy_b(b_t, b);

	if (DIGITS_B(b) == 0)//如果b是0,除0错误
		return FLAG_DIVZERO;
	else if (DIGITS_B(r_t) == 0)//如果a=0
	{
		SETZERO_B(q);
		SETZERO_B(rem);
		return FLAG_OK;
	}
	else if (cmp_b(r_t, b_t) == -1)//如果a<b,返回a就好了
	{
		cpy_b(rem, r_t);
		SETZERO_B(q);//商为0
		return FLAG_OK;
	}
	else if (cmp_b(r_t, b_t) == 0)//如果a=b,返回1就好了
	{
		q[0] = 1; q[1] = 1;//商为1
		SETZERO_B(rem);//余数为0
		return FLAG_OK;
	}
	else if (DIGITS_B(r_t) == 0)//如果a=0，非常好
	{
		SETZERO_B(q);
		SETZERO_B(rem);
		return FLAG_OK;
	}
	//如果除数位数为1，使用移位除法，后面要求除数大于等于两位
	if (DIGITS_B(b_t) == 1)
	{
		mydiv_b(r_t, b_t, q, rem);
		return FLAG_OK;
	}


	mbptr = MSDPTR_B(b_t);
	bn_1 = *mbptr;

	while (bn_1 < BASEDIV2)
	{
		bn_1 = bn_1 << 1;
		d++;
	}
	uint64_t shiftr =(int)(BITPERDIGIT - d);//左移时配合的右移次数

	if (d > 0)
	{
		bn_1 = bn_1 + (uint32_t)((uint64_t)*(mbptr - 1) >> shiftr);

		if (DIGITS_B(b_t) > 2)
			bn_2 = (uint32_t)((uint32_t)((uint64_t)(*(mbptr - 1)) << d) + (uint32_t)((uint64_t)(*(mbptr - 2)) >> shiftr));
		else//等于两位则直接赋值
		{
			bn_2 = (uint32_t)((uint64_t)(*(mbptr - 1)) << d);
		}
	}
	else//没有移动则原封不动
	{
		bn_2 = (uint32_t)(*(mbptr - 1));
	}

	mbptr = MSDPTR_B(b_t);
	mrptr = MSDPTR_B(r_t) + 1;//指向滑动窗口最高位，填了一位，之后可能是0可能不是，由后面移位决定
	rptr = MSDPTR_B(r_t) - DIGITS_B(b_t) + 1;//指向滑动窗口最低位
	qptr = q + DIGITS_B(r_t) - DIGITS_B(b_t) + 1;//最高即将可能为0可能不是，由后面移位决定
	*mrptr = 0;//初始化为0



	while (rptr >= LSDPTR_B(r_t))//开始做事情，下标对应，r(m+n)对应a(m)，r(m+n-1)对应a(m-1)
	{

		ri = (uint32_t)((uint32_t)((uint64_t)(*mrptr) << d) + (uint32_t)((uint64_t)(*(mrptr - 1)) >> shiftr));
		ri_1 = (uint32_t)(((uint32_t)((uint64_t)(*(mrptr - 1)) << d)  + (uint32_t)((uint64_t)(*(mrptr - 2)) >> shiftr)));

		if (mrptr - 3 > r_t)//不止有3位
		{
			ri_2 = (uint32_t)(((uint32_t)((uint64_t)(*(mrptr - 2)) << d) + (uint32_t)((uint64_t)(*(mrptr - 3)) >> shiftr)));
		}
		else//只有2位
		{
			ri_2 = (uint32_t)((uint64_t)(*(mrptr - 2)) << d);
		}


	    //计算q估计，kuo是64位的，取整数部分
		kuo = (uint64_t)((((uint64_t)ri << BITPERDIGIT) + (uint64_t)ri_1)  / bn_1);

		if (kuo < ALLONE)//选择小的那个，最大估计的q不会超过B-1，不会超过B进制
		{
			qh = (uint32_t)kuo;

		}
		else
			qh = ALLONE;

		kuo = ((uint64_t)ri << BITPERDIGIT) + (uint64_t)ri_1 - (uint64_t)bn_1*(uint64_t)qh;
		if (kuo < BASE)//如果括号中的都比B大了，左边肯定b[n-2]*q小于B*B,小于才比较，大于会溢出，也比较不了
		{
			right = (uint64_t)(kuo << BITPERDIGIT) + (uint64_t)ri_2;
			left = (uint64_t)bn_2 *(uint64_t)qh;
			if (left > right)//没有溢出，需要比较
			{
				qh--;
				//重复这一过程
				kuo = ((uint64_t)ri << BITPERDIGIT) + (uint64_t)ri_1 - (uint64_t)bn_1*(uint64_t)qh;
				if (kuo< BASE)//如果括号中的都比B大了，左边肯定b[n-2]*q小于B*B,小于才比较，大于会溢出，也比较不了
				{
					right = (uint64_t)(kuo << BITPERDIGIT) + (uint64_t)ri_2;
					left = (uint64_t)bn_2 *(uint64_t)qh;
					if (left > right)//没有溢出，需要比较
						qh--;
				}
				else
				{
					//printf("rh = %llX\n", rh);
					//printf("rh*B = %llX\n",rh <<BITPERDIGIT);
					//printf("rh*B+ri-2 = %llX\n", (rh << BITPERDIGIT) + (uint64_t)ri_2);
				}

			}
		}

		borrow = BASE;//借位默认
		carry = 0;
		int i = 0;
		for (bptr = LSDPTR_B(b_t), rcir = rptr; bptr <= mbptr; bptr++, rcir++)
		{
			if (borrow>= BASE)// 没有发生借位的情况，最高处总是BASE
			{			
				carry = (uint64_t)qh * (uint64_t)*bptr + (uint64_t)(uint32_t)(carry >> BITPERDIGIT);
				borrow = (uint64_t)*rcir + BASE - (uint64_t)(uint32_t)carry;
				*rcir = (uint32_t)(borrow);
			}
			else//发生了借位,恢复
			{
				carry = (uint64_t)qh * (uint64_t)*bptr + (uint64_t)(uint32_t)(carry >> BITPERDIGIT);
				borrow = (uint64_t)*rcir + (uint64_t)BASE - (uint64_t)(uint32_t)carry - 1ULL;//借位了
				*rcir = (uint32_t)(borrow);
			}
		}

		if (borrow >= BASE) //没有发生借位的情况，最高处总是BASE
		{
			borrow = (uint64_t)*rcir + BASE - (uint64_t)(uint32_t)(carry >> BITPERDIGIT);
			*rcir = (uint32_t)(borrow);
		}
		else//发生了借位
		{
			borrow = (uint64_t)*rcir + BASE - (uint64_t)(uint32_t)(carry >> BITPERDIGIT) - 1ULL;//借位了
			*rcir = (uint32_t)(borrow);
		}
		*qptr = qh;

		if (borrow < BASE)//borrow还被借位了没有还回来,q大了一位，b加会到a中，和加法过程一样
		{
			carry = 0;
			for (bptr = LSDPTR_B(b_t), rcir = rptr; bptr <= mbptr; bptr++, rcir++)
			{
				carry = ((uint64_t)*rcir + (uint64_t)*bptr + (uint64_t)(uint32_t)(carry >> BITPERDIGIT));
				*rcir = (uint32_t)carry;
			}
			*rcir = *rcir + (uint32_t)(carry >> BITPERDIGIT);
			*qptr = *qptr - 1;
		}
		qptr--; mrptr--; rptr--;
	}
	SETDIGITS_B(q, DIGITS_B(r_t) - DIGITS_B(b_t) + 1);
	RMLDZRS_B(q);
	SETDIGITS_B(r_t, DIGITS_B(b_t));
	cpy_b(rem, r_t);
	return FLAG_OK;

}