class="markdown_views prism-atom-one-light">   在(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 { 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 ", rh); //printf("rh*B = %llX ",rh < //printf("rh*B+ri-2 = %llX ", (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; }