大数除法之迭代法

落叶

上一篇文章讲到了估商法的原理，有了一个不错的效率，但在要求精度较大时，速度和迭代法相比差距很大。
      除法：u/y=u*(1/y)；
      先讲一下倒数迭代式：x1=(2-y*x0)*x0,x0是y的倒数的近似值，它必须要小于y的倒数。另外迭代式中的乘法子程序要选用快速乘法（如FFT算法的乘法子程序）。
      否则迭代法的除法速度是很慢的，远小于估商法。
     以求9的倒数为例演示迭代法的使用，求9的32精度的倒数：
     1/9=1.1111111111111111111111111111111e-1;这里计算的是9的32位精度的倒数
    (2-9*0.11)*0.11=1.111e-1;                  现在我们用初值0.11迭代计算32位精度的倒数
    (2-9*0.1111)*0.1111=1.1111111e-1;
    (2-9*0.11111111)*0.11111111=1.111111111111111e-1;
    (2-9*1.111111111111111e-1)*1.111111111111111e-1 = 1.1111111111111111111111111111111e-1;
    对于如何确定y的初值，网上找不到相关的文献，我最早是用估商法算y的前N位的初值，然后迭代算出精度更高的倒数。
     但是这种方法在综合运算时，估商法存在大量的进制转换，和指数对位，效率不高。
    前一段时间才突破用计算器的cpu的除法功能计算除数倒数的初值，这个问题原理很简单：（就是取除数的前几位有效数，用cpu算倒数，然后取倒数的前几位有效位，转换成高精表现形式，用这个高精形式的倒数作为迭代式的初值，可以迭代出你所需的更高精度的倒数）
   但实现起来还是相当复杂的，三分的算法原理，七份的调试。
   
   利用cpu 算初值，由于cpu除法的局限性，过小的或过大的数它都计算不了，所以这些数你必须调整规范为0.xxx...*10^n的形式，其中xxx必须为有效数字，1234.5678可以表示为0.12345678*10^4，但不能为0.012345678*10^5，这都是调试才能理解的技术，不调试你很难理解，
  选取除数的前几个有效数字用cpu算倒数，这里采用1/1234，算出的得数是8.1037277147487844408427876823339e-4，现在我们取其中的8.1037e-4，具体应该取多少位有效数字，可以调试得出，然后把8.1037e-4识别并换算成高精度数的表示形式（这里还要考虑转换后的指数位，这个也需要慢慢的调试），然后用转换好的数进行后面的迭代运算，后面的迭代运算需要有理数的加，减，乘支持。
  纯迭代法实现的除法代码更简洁，效率更高，特别是作为子程序调用时可以减少一些中间转换，在综合运算中效率要高于估商配迭代法的算法，是个很实用的算法。

kastin

“不调试你很难理解”，你如果学过《计算机组成原理》，就不会这么说。

落叶

不好意思，是我用词不当，我当时是不断调试，不断改进才把这个问题解决的，所以我解说时也只是还原了当时的情形，用cpu算初值，我很早就有这个想法，一直难以实现，这次也是不断调试测试才得也解决，所以这次用词上对调试说的重了一些，《计算机组成原理》这本书我曾学过，不过学的很浅，只懂一点大致的原理。

zeroieme

借楼讨论下迭代求倒的优化。
牛顿迭代是精度倍翻的，怎么利用这个特点减少乘法计算量呢。

一，我实现过。 \(y*x_i \approx 1\)，将\(x_1=(2-y*x_0)*x_0\)替换成\(x_1=(1-y*x_0)*x_0+x_0\)，可以把每步迭代中第二次乘法长度减少一半。
二，在这里(数学研发论坛)见过：把y实行逐次增加精度：\(y_0\)、\(y_1\)、\(y_2\)……。这样把每步迭代中第一次乘法长度减少。

三，结合以上优化。产生了一个设想，不知道是否实现的实例。因为原除数y是逐次增加精度进入运算，除数倒数x也是逐次增加精度。就是\(y_{i+1}=y_i + \Delta y_{i+1}\)  ，\(x_{i+1}=x_i + \Delta x_{i+1}\)  有什么办法充分利用之前步骤的计算结果。


谢谢站长帮助

liangbch

我们知道，使用牛顿迭代法求倒数属于浮点运算，或者为非完全精度运算。浮点运算和整数运算相对，整数运算属于完全精度运算，意为每一步运算都取全部精度，不丢弃任何数字，结果也一定是精确的。浮点运算的结果是非精确的，只需要达到指定的精度即可。

和整数运算相比，浮点运算最麻烦的在于精度控制，如果精度控制做的太保守（即在运算过程中保留的位数太多）则会导致运算时间过长。反之，则会导致最终结果达不到要求的精度。所以，在运算过程中需要保留多少位数字，才能恰到好处，需要精心设计和验证。

牛顿迭代法求倒数，理论上，每迭代一次，精度加倍。如果你在每次迭代时，源操作精度仅仅保持目标精度的一半，你会发现，最终精度是达不到要求的。

例如，最终精度需要10000比特，依此倒推，最后一次迭代的目标精度为10000，倒数第2次迭代的目标精度为5000。
故前几次迭代的源精度和目标精度依次为（假定x0 可通过其他指令求得，且精度可达到30比特）
30，  40
40，  79
79，  157
157， 313
313， 625
625， 1250
1250，2500
2500，5000
5000，10000

如果你真的这样做，你会发现，最终精度是达不到要求的。事实上，为了保证精度，需要的目标精度必须稍微大于理论精度

我的做法是：
首先，通过最终目标精度，推导出每次迭代过程中需要的精度，并将其存入数组。
精度用2个成员(i,b)来表示，i表示多少个limbs,b表示多少个比特。（如果采用10^k进制，i表示多少个limbs,b表示多少位）
例如采用\( 2^{30} \)进制，精度为50比特,则i=50/30=1, b=50 mod 30=20

下面是我的代码

1. #define DO_FORMAT(res)        \
   
2. do {                                        \
   
3.         if ( res.b>=BCMP_NUMB_BITS)                        \
   
4.         {        res.i++; res.b-=BCMP_NUMB_BITS; }        \
   
5. } while (0)
   
6. 
   
7. 
   
8. typedef struct _precision_st
   
9. {
   
10.         int i;
    
11.         int b;
    
12. }PRC_ST;
    
13. 
    
14. #define EXTRA_BITS         2
    
15. #define BCMP_NUMB_BITS  30
    
16. 
    
17. PRC_ST tn_prc_arr[32];                // the precision for each iteration
    
18. PRC_ST t_prc;                        // the target precision
    
19. PRC_ST s_prc;                        // the source precision
    

复制代码

首先，将 tgt_prec赋值为目标精度+EXTRA_BITS 比特，然后执行下面的代码

1. 
   
2. t_prc = tgt_prec;
   
3. i=0; tn_prc_arr[0]=t_prc;                // tgt_prec have included EXTRA_BITS and need not add EXTRA_BITS
   
4. while (1)
   
5. {
   
6.         m= 0-(t_prc.i & 1);                // if t_prc.i is odd number, then m= 0xffffffff, else m=0
   
7.         t_prc.i /= 2;
   
8.         t_prc.b= (t_prc.b+1 + (m & BCMP_NUMB_BITS))/2;
   
9.         t_prc.b += EXTRA_BITS;                // enlarge target precision to guarantee to get enough precision
   
10.         DO_FORMAT(t_prc);
    
11.         if ( t_prc.i==0 )
    
12.         break;
    
13.         tn_prc_arr[++i]=t_prc;
    
14. }

复制代码

我的使用牛顿迭代方求x=1/a的大致流程
步骤1: s1=a * x0,
步骤2: 如果(s1>1), 那么s1=s1-1, doAdd=true, 否则s1=1-s1, doAdd=false
步骤3: s2= s1 * x0
步骤4: 如果 doAdd, 那么 x1= x0 + s2, 否则x1=x0-s2
步骤5: 将x1赋值给x0,此时x0的精度较原始值加倍

下面是迭代过程的精度控制

1. 
   
2. for (;i>=0;i--)
   
3. {
   
4.         bcmp_size_t SN;                // the needed source precison in limb
   
5.         bcmp_size_t TN;                // the needed target precison in limb
   
6.         bcmp_size_t TN1;        // TN1=TN+1
   
7.         bcmp_size_t zCnt;        // the lead zero count
   
8.         BOOL doAdd;
   
9. 
   
10.         t_prc = tn_prc_arr[i];        // t_prc, target precision
    
11.                
    
12.         m= 0-(t_prc.i & 1);        // if t_prc.i is odd, m=0xffffffff, else m=0
    
13.         s_prc.i = t_prc.i/2;
    
14.         s_prc.b = ( (m & BCMP_NUMB_BITS) + t_prc.b+1)/2;       
    
15.         DO_FORMAT(s_prc);
    
16. 
    
17.         SN = s_prc.i + (s_prc.b>0);               
    
18.         TN = t_prc.i + (t_prc.b>0);
    
19.         TN1 = TN+1;
    
20.         .....        此处略
    
21. }

复制代码

下面是具体的精度控制方法

1. 在步骤1的计算过程中，a的精度取TN1个limb,x0的精度取SN个limb,结果精度值保留TN1个limb，丢弃低位部分
2. 在步骤2中，高位部分是多个连续的0，需要统计连续0的个数zCnt，和剩余部分的长度s1_n，其总长度为TN1
3. 在步骤3中，s1的精度为s1_n, x0的精度取SN，结果精度取TN-zCnt，丢弃低位部分。
4. 在我的代码中x1和x0共用相同的缓冲区，故步骤4是需要将x0的长度扩展至TN1，低位部分补0，然后加上（或者减去）s2

其他的注意事项。
1. a要做预处理，使得1/a的值>=0.5且小于1,这样，在迭代的过程中，参与运算的非0比特位尽可能多。
2. 在步骤4中，如果发现x1>1.0, 要将其重新格式化为0.9999999...., 如果是2进制表示，将所有bit位置1

liangbch

大数除法，至少有3种，除了你说的估商法和牛顿迭代法，还有分治法。

牛顿迭代法也是很有讲究的，为了追求最大的效率，求得的倒数的精度取多少，需要精心的实验，一般取除数长度的1/2+1。
实验表明，只要在除数长度和商的精度很大时，牛顿迭代法才是做好的选择，否则，分治法可能是最好的选择。只有当除数的长度很小时，估商法才是最优的。
编程时，相对于编程者，任何复杂的算法，都需要调试。就某一具体算法，如果编程者对此很熟练，则很少或者不需要调试，而对于另一个不太懂的编程者，则需要调试。复杂是相对的，是相对于编程者的，不是绝对的。

liangbch

但是这种方法在综合运算时，估商法存在大量的进制转换，和指数对位，效率不高。

这不是估商法的问题，是你设计的问题。

最重要的是，对于大数（或者任意精度）运算，只有在最终输出的时候才需要进制转化，才需要转化为字符串表示。大数的内部表示，必须采用更高的进制,比如，\( 10^9,2^{30},2^{32},2^{64} \)

liangbch

对于如何确定y的初值，网上找不到相关的文献

，本站就有啊，请参阅
http://bbs.emath.ac.cn/forum.php?mod=viewthread&tid=5970

落叶

liangbch 发表于 2017-12-27 12:19
这不是估商法的问题，是你设计的问题。

最重要的是，对于大数（或者任意精度）运算，只有在最终输出的 ...

我开始设计程序时，因没有解决万进制下的小数点对位问题，所以各个子程序之间的数据交换采用的是十进制，子程序运行时转换成自己所需的万进制，运算完后再还原为十进制数，所以来回转换浪费很多，最后在疯狂计算器的作者帮助下解决了万进制下的小数对位，并为三角函数运算重写了专用加减乘除，效率有了一倍的提升。
解决万进制下的小数对位采用的方法就是把操作数转换成纯小数形式，我没能力根据纯小数理论解决万进制下的小数对位，所以我是在大量的调试基础上，慢慢修改完善了这个问题。
具体把十进制转换成效率更高的十六进制进行运算，我开始考虑把十进制数转换成纯整数形式再进行转换，但存在极大数和极小数难以操作的问题，假如把十进制数转换成纯小数后，再转换成十六进制数，我并不知道如何操作？

落叶

zeroieme 发表于 2017-12-25 21:13
借楼讨论下迭代求倒的优化。
牛顿迭代是精度倍翻的，怎么利用这个特点减少乘法计算量呢。

在迭代法乘法优化中，我采用了类似的优化，我的方法如下：
因为我是在为我的高精表达式计算器写的子程序，很早我都加入了精度控制功能，剔除掉多余的运算，
我设定了一个全局变量，保存运算中的设定精度，运算子程序每次运算完后，会根据设定精度把运算结果多出的精度除掉（如设定精度小数点后20位，运算结果如果多于20位，会把多余的长度除掉），子程序在开始计算时会检验一次操作数的精度，例：两个数相乘，设定精度20，传入两个30位操作数，程序开始会把操作数修正为21位，运算结果为42位，再次修正为21位，返回运算结果。
在进行迭代循环前，迭代初值前期处理为有效位20位，迭代循环次数根据所需精度和初值精度提前算出，好处是减少一次迭代运算，看你之前的文章你也提到这个问题，
第一次迭代循环前，保存全局精度，全局精度设为初值精度的两倍，
进入第一次循环，此时，参入运算的y的精度是40，x0的精度为20，第一次乘法是40位乘以20位，此时会产生60位的结果，修正为40位，第二次乘法也是40位乘以20位，此时也会产生60位的结果，修正为40位，循环完的结果是40位的数，它实际上是下一次循环的x0,然后设定全局精度为40*2=80，
进入第二次循环，此时，参入运算的y的精度是80，x0的精度为40，第一次乘法是80位乘以40位，此时会产生120位的结果，修正为80位，第二次乘法也是80位乘以40位，此时也会产生120位的结果，修正为80位，循环完的结果是80位的数，然后设定全局精度为80*2=160..........略。