大数除法之迭代法

liangbch

误差的原因是进位产生的，进位的不确定性也决定了冗余精度的不确定性

这样的描述不准确。总体来说，物理中的误差是由测量引起的。数学中的误差是因为，绝大多数有理数的小数表示是无限循环小数的形式，所有的无理数都是无限不循环小数，故用小数的形式（不管是什么进制）来表示数，一定会有误差。

liangbch

误差分为绝对误差和相对误差。比如，误差小于\( 2^{-100} \)是指绝对误差。又如，精确到4位有效数字，这是指相对误差。

落叶

确实不准确，当时想了半天，没找到合适的用语。应该说该例中出现的问题和进位与否有一定关系。编程中的精度控制是误差形成的一个主因。不同的精度控制会形成不同的误差特征，有些特征更利于编程开发。
大数运算中的精度控制是指相对误差。这种说法才是正确的。

liangbch

冗余精度的做法广泛存在于 计算过程，和数的存储过程中。

先说数的存储中的冗余精度，如果进制为2^b, 即一个limb存储b比特，则对于一个x比特的数来说， 最多需要 (x + b-1)/b个limb就够了。但你检查GMP的代码，__GMPF_BITS_TO_PREC的定义如下
#define __GMPF_BITS_TO_PREC(n)                                                \
  ((mp_size_t) ((__GMP_MAX (53, n) + 2 * GMP_NUMB_BITS - 1) / GMP_NUMB_BITS))
显然，mpf_t 类型的变量的存储是增加了冗余的。

相应的，如果进制为10^k, 则对于一个x位的数来说， 最多需要 (x +k-1)/k个limb就够了。例如
1,2,3,4,5,6,7,8 位10进制数，采用万进制存储，分别需要1,1,1,1,2,2,2,2个limb

再说计算过程，很多浮点函数的计算涉及到级数求和，在计算过程中有累加误差。所以可能在计算过程中使用比最终结果的表示更高的精度。
关于舍入方式，mpfr的每个函数都要求给出舍入方法，
舍入方式包括以下5种
MPFR_RNDN: 舍入到最接近的数
MPFR_RNDZ: 向0方向舍入
MPFR_RNDU: 向正无穷舍入
MPFR_RNDD: 向负无穷舍入
MPFR_RNDA: 远离0

如果在计算过程中和最终存储是都保留了额外的冗余。那么，在计算过程中无需考虑四舍五入，毕竟直接丢弃低位更方便，而四舍五入需要考虑进位。如果数的存储没有冗余，则四舍五入更准确些。GMP的输出函数，即数转化为字符串函数 mpf_get_str 在输出是总是四舍五入，故内部存储为24.99 和25.01，在输出时，如指定2位数字，则总是会输出25.

liangbch

一般的，我们只能尽可能使最终结果达到指定的精度，但完全做到很难。四则运算中，加，乘, 除这3种运算， 一般的，结果的精度大体等于源操作数的精度，或者说，源操作数中，精度较低的那个数的精度。但减法例外，如果2个数很接近， 如a=355/113=3.1415929203539823008849557522124, b=3.1415926535897932384626433832795, a,b到保持10位精度，但c=a-b就只有3位有效数字了。尽管mpfr提供了几种舍入模式，即保证最低bit是正确的。但mpfr依然不能追踪结果的精度。原文"MPFR does not keep track of the accuracy of a computation", MPFI 应该可以做到这一点。

Ickiverar

@只是呼吸
我确实是那个qsd573。不光是qsd573，qsd572和qsd574也都是我。

只是呼吸

  以下是我自己写的牛顿迭代除法函数。
       除法函数的入口条件：
                   要添加<math.h> 头文件。
                   qut[]-------存储商。商为\(10^8\)进制。
                   temp[]-----临时变量(存储中间结果）
                   Work[]-----临时变量（存储倒数）
                   Cdiv[]------规约化了的被除数。采用\(10^8\)进制。LENdividend 是 Cdiv[]的长度减1.  
                   divR[]------规约化了的除数。采用\(10^8\)进制。LENdivR 是 divR[]的长度减1.其中divR[]的最高项一定要有9位十进制数
                 当除数是一个单精度的 unsigned  long 时，调用其它函数处理，不进入 Newton_div()。
                  

1. void Newton_div(UINT32 * const qut,UINT32 * const temp,UINT32 * const Work,UINT32 * const Cdiv,int LENdividend, UINT32 * const divR,int LENdivR)
   
2. {
   
3.         UINT64 h64=100000000000000000;//17个 0
   
4.         int len,i,s,N,LEN2,LenWork;
   
5. 
   
6.         len=LENdividend-LENdivR+1;//计算商的有效位数len=LL-LR       
   
7.         N=(int)(floor(log(len)/log(2)));  //计算迭代的次数 N=floor(LOG(Len)/LOG(2)); ;        
   
8. 
   
9.         Work[0]=(UINT32)(h64/divR[LENdivR]);        //求初商，估商s0=100000000 00000000/divR[LR];
   
10.         Work[0]/=10;
    
11.         Work[0]--;
    
12. 
    
13.               LenWork=0;
    
14.                   len=1;
    
15.                   
    
16.                   for(i=0;i<=N;i++)
    
17.                   {
    
18.                           ALU_mul(qut,divR,Work,LENdivR+1,LenWork+1);// 1 乘法函数：求a*s0  
    
19.                           
    
20.                           s=LENdivR+LenWork+1;
    
21.                           s=Newton_LeftZero(qut,s);// 寻找 qut[] 的非零最高位的下标 函数
    
22.   
    
23.                           Newton_sub(qut,s); // 2 减法函数，  2-a*s0
    
24. 
    
25.                           ALU_mul(temp,qut,Work,s+2,LenWork+1);   ///   乘法(2-a*s0)*s0
    
26.                          
    
27.                           LEN2=s+LenWork+2;
    
28.                           LEN2=Newton_LeftZero(temp,LEN2);// 寻找 temp[] 的非零最高位的下标 函数
    
29.                           Newton_clear(Work,LenWork+2);  //清零函数  把 Work[] 清零
    
30.                           len<<=1;                                 // 精度加倍
    
31.                           LenWork=len;// 计算迭代后的精度，也正好是Work 的长度
    
32.        
    
33.                           Newton_set(Work,temp,LEN2,LenWork);///赋值，运算关系Work<--temp, 赋的长度是LenWork, 裁剪掉temp 中的多余的右边的数(低位）。
    
34.                           
    
35.                           Newton_Zero(qut,s+4);   //  清零函数  把 qut[] 清零
    
36.                           Newton_Zero(temp,LEN2+4);  //清零函数  把 temp[] 清零
    
37. 
    
38.                   }///迭代的边界
    
39. 
    
40.        
    
41.         len=LENdividend-LENdivR+1;
    
42.         for(i=0;i<=len-1;i++)
    
43.         {
    
44.                 Work[i]=Work[LenWork+1+i-len];//取Work左面的len位，赋给Work.
    
45.         }
    
46.                 for(i=len;i<=LenWork;i++)  ///将多余的数清零
    
47.                         Work[i]=0;
    
48.        
    
49.       ALU_mul(qut,Cdiv,Work,LENdividend+1,len);  ///  乘法，用倒数乘 被除数得到商。qut[] 存储商
    
50. 
    
51.         s=Newton_LeftZero(qut,LENdividend+len+1);//  寻找 qut[] 的非零最高位的下标 函数
    
52.     len=quotient_Num(Cdiv,divR,LENdividend,LENdivR,qut[s]);///  求商的精确长度的函数，（不是把商求出来，只是求出商的精确长度）
    
53. 
    
54. 
    
55.         for(i=0;i<=len-1;i++)
    
56.         {
    
57.                 qut[i]=qut[s+i-len+1];//取Work左面的s位，赋给Work.
    
58.         }
    
59.        
    
60.         for(i=len;i<=s;i++)  ///将多余的数清零
    
61.                 qut[i]=0;
    
62.       
    
63.    // return 0;
    
64. }//*///

复制代码

乘法函数就是我写的那个“大整数乘法代码”，把它做了优化后暂时用着。
减法函数是专门针对牛顿除法的.代码如下：

1. static int Newton_sub(UINT32 *dst,int size2)///
   
2. {
   
3.         int i,s=0;
   
4.        
   
5.         while(dst[s]==0)
   
6.         {
   
7.                 s++;
   
8.         }
   
9.         dst[s]=100000000-dst[s];
   
10.        
    
11.                 for(i=s+1;i<=size2;i++)
    
12.                 {
    
13.                         dst[i]=99999999-dst[i];
    
14.                 }
    
15.        
    
16.         dst[size2+1]=1;//(///
    
17.         return 0;
    
18. }////////
    

复制代码

其余的函数都没有什么特殊的，就不写出来了。

这个牛顿除法终于完工，花了很多时间，一开始的时候，一个数字都不对。后来才慢慢正确。

实际运算：35800位十进制数  除以   11000位十进制数，我写的牛顿除法用时  35毫秒。我写的试商除法用时   110毫秒。

可能是乘法不好，这么大的数字用硬乘法也许勉强了。等下一步写一个karatsuba乘法试试。

落叶

只是呼吸 发表于 2018-1-13 00:56
以下是我自己写的牛顿迭代除法函数。
       除法函数的入口条件：
                   要添加 头文件 ...

请教一下！
如果把除数123456789代入
Work[0]=(UINT32)(h64/divR[LENdivR]);        //求初商，估商s0=100000000 00000000/divR[LR];
        Work[0]/=10;
        Work[0]--;
的具体数为：
10000000000000000/123456789/10-1 = 8099999;
到这段代码时 ALU_mul(qut,divR,Work,LENdivR+1,LenWork+1);// 1 乘法函数：求a*s0  
Work中存的是 8099999，它是个正整数
那么是不是这段代码
s=LENdivR+LenWork+1;
s=Newton_LeftZero(qut,s);// 寻找 qut[] 的非零最高位的下标 函数
把两个正整数的乘法转换成了一个数乘以它的倒数，这里面的细节是怎么通过调试实现的？你没有专门设变量存储小数位，难道是用乘法产生的前导零或尾零来标记？

只是呼吸

那么是不是这段代码
s=LENdivR+LenWork+1;
s=Newton_LeftZero(qut,s);// 寻找 qut[] 的非零最高位的下标 函数
把两个正整数的乘法转换成了一个数乘以它的倒数，这里面的细节是怎么通过调试实现的？你没有专门设变量存储小数位，难道是用乘法产生的前导零或尾零来标记？

这个牛顿除法从一开始就没有用小数。也就是说从一开始就是用整数来模拟小数迭代的。
求 \(a\) 的倒数的牛顿迭代公式是：\((2-a*s_{i})s_{i}\) 其中 \(s_{i}\) 是逐渐精确的 \(a\) 的倒数。你例子中的\(8099999\)  就是\(s_{0}\), 这个数运行到 第\(31\) 行时。精度就加倍。我们只取加倍后的\(16\)个数字,比如说是\(8099999     77777777\)，作为 \(s{1}\),进入下一次运算（其他的数字全部抛弃),实现这个目的的函数就是第 \(33\) 行的那个。

落叶

你的迭代法确实很不错，可是说绝大部分都是调试找规律的结晶，凭想象凭原理都是很难实现的，另外也很佩服你在没有理论（宝宝关于估商法的理论）的支持下，加长除数长度增加估商正确率的估商法，我也是在宝宝的理论和你的程序提示下才对估商法有了新的认识。
我的迭代法也是规约了除数的，和你的原理差不多，但没有规约被除数，我在计算倒数的过程中，有些关键环节也是在作整数运算，但我增加了逻辑上的约束，就是始终把运算结果当作小数看待，我有专门的指数位监视中间的小数点变化。进入最终运算环节的也是小数。
你的迭代法通过规约除数，被除数，把后续的运算控制在易于编程调试的范围。
我在编写高精表达式计算器中大量的用到字符串规范化，和高精数的规范化，和高精数范围缩小化，好处就是逻辑清晰（这个很重要，逻辑清晰在整个程序设计中都是非常重要的！），程序设计简单易调试。