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;
}