数学研发论坛

 找回密码
 欢迎注册
楼主: 落叶

[原创] 大数除法之迭代法

[复制链接]
发表于 2018-1-1 19:56:26 | 显示全部楼层
误差的原因是进位产生的,进位的不确定性也决定了冗余精度的不确定性


这样的描述不准确。总体来说,物理中的误差是由测量引起的。数学中的误差是因为,绝大多数有理数的小数表示是无限循环小数的形式,所有的无理数都是无限不循环小数,故用小数的形式(不管是什么进制)来表示数,一定会有误差。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-1-1 19:59:24 | 显示全部楼层
误差分为绝对误差和相对误差。比如,误差小于\( 2^{-100} \)是指绝对误差。又如,精确到4位有效数字,这是指相对误差。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2018-1-1 20:23:00 | 显示全部楼层
确实不准确,当时想了半天,没找到合适的用语。应该说该例中出现的问题和进位与否有一定关系。编程中的精度控制是误差形成的一个主因。不同的精度控制会形成不同的误差特征,有些特征更利于编程开发。
大数运算中的精度控制是指相对误差。这种说法才是正确的。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-1-1 20:34:04 | 显示全部楼层
冗余精度的做法广泛存在于 计算过程,和数的存储过程中。

先说数的存储中的冗余精度,如果进制为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.

点评

没想到舍入方式还有这么多这么细,我的舍入太粗糙了!  发表于 2018-1-1 21:09
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-1-1 20:54:51 | 显示全部楼层
一般的,我们只能尽可能使最终结果达到指定的精度,但完全做到很难。四则运算中,加,乘, 除这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 应该可以做到这一点。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-1-8 13:53:55 | 显示全部楼层
@只是呼吸
我确实是那个qsd573。不光是qsd573,qsd572和qsd574也都是我。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-1-13 00:56:52 | 显示全部楼层
本帖最后由 只是呼吸 于 2018-1-13 01:38 编辑

  以下是我自己写的牛顿迭代除法函数。
       除法函数的入口条件:
                   要添加<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.         len=LENdividend-LENdivR+1;//计算商的有效位数len=LL-LR       
  6.         N=(int)(floor(log(len)/log(2)));  //计算迭代的次数 N=floor(LOG(Len)/LOG(2)); ;        

  7.         Work[0]=(UINT32)(h64/divR[LENdivR]);        //求初商,估商s0=100000000 00000000/divR[LR];
  8.         Work[0]/=10;
  9.         Work[0]--;

  10.               LenWork=0;
  11.                   len=1;
  12.                   
  13.                   for(i=0;i<=N;i++)
  14.                   {
  15.                           ALU_mul(qut,divR,Work,LENdivR+1,LenWork+1);// 1 乘法函数:求a*s0  
  16.                           
  17.                           s=LENdivR+LenWork+1;
  18.                           s=Newton_LeftZero(qut,s);// 寻找 qut[] 的非零最高位的下标 函数
  19.   
  20.                           Newton_sub(qut,s); // 2 减法函数,  2-a*s0

  21.                           ALU_mul(temp,qut,Work,s+2,LenWork+1);   ///   乘法(2-a*s0)*s0
  22.                          
  23.                           LEN2=s+LenWork+2;
  24.                           LEN2=Newton_LeftZero(temp,LEN2);// 寻找 temp[] 的非零最高位的下标 函数
  25.                           Newton_clear(Work,LenWork+2);  //清零函数  把 Work[] 清零
  26.                           len<<=1;                                 // 精度加倍
  27.                           LenWork=len;// 计算迭代后的精度,也正好是Work 的长度
  28.        
  29.                           Newton_set(Work,temp,LEN2,LenWork);///赋值,运算关系Work<--temp, 赋的长度是LenWork, 裁剪掉temp 中的多余的右边的数(低位)。
  30.                           
  31.                           Newton_Zero(qut,s+4);   //  清零函数  把 qut[] 清零
  32.                           Newton_Zero(temp,LEN2+4);  //清零函数  把 temp[] 清零

  33.                   }///迭代的边界

  34.        
  35.         len=LENdividend-LENdivR+1;
  36.         for(i=0;i<=len-1;i++)
  37.         {
  38.                 Work[i]=Work[LenWork+1+i-len];//取Work左面的len位,赋给Work.
  39.         }
  40.                 for(i=len;i<=LenWork;i++)  ///将多余的数清零
  41.                         Work[i]=0;
  42.        
  43.       ALU_mul(qut,Cdiv,Work,LENdividend+1,len);  ///  乘法,用倒数乘 被除数得到商。qut[] 存储商

  44.         s=Newton_LeftZero(qut,LENdividend+len+1);//  寻找 qut[] 的非零最高位的下标 函数
  45.     len=quotient_Num(Cdiv,divR,LENdividend,LENdivR,qut[s]);///  求商的精确长度的函数,(不是把商求出来,只是求出商的精确长度)


  46.         for(i=0;i<=len-1;i++)
  47.         {
  48.                 qut[i]=qut[s+i-len+1];//取Work左面的s位,赋给Work.
  49.         }
  50.        
  51.         for(i=len;i<=s;i++)  ///将多余的数清零
  52.                 qut[i]=0;
  53.       
  54.    // return 0;
  55. }//*///
复制代码



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

  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乘法试试。


点评

抱歉:程序的第29行写错了,应是 Newton_Zero(); 我自己在网站的编辑框中改了清零函数的名称,忘了把这一行也改了。  发表于 2018-1-13 12:13
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2018-1-13 14:02:31 | 显示全部楼层
本帖最后由 落叶 于 2018-1-13 14:07 编辑
只是呼吸 发表于 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[] 的非零最高位的下标 函数
把两个正整数的乘法转换成了一个数乘以它的倒数,这里面的细节是怎么通过调试实现的?你没有专门设变量存储小数位,难道是用乘法产生的前导零或尾零来标记?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-1-13 15:23:31 | 显示全部楼层
本帖最后由 只是呼吸 于 2018-1-13 15:36 编辑
那么是不是这段代码
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\) 行的那个。

点评

80999999 运行到程序的第25行时,精度加倍,其16位数是 8100000073710000,把它作为下一次迭代的初值  发表于 2018-1-13 20:31
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2018-1-13 17:03:06 | 显示全部楼层
你的迭代法确实很不错,可是说绝大部分都是调试找规律的结晶,凭想象凭原理都是很难实现的,另外也很佩服你在没有理论(宝宝关于估商法的理论)的支持下,加长除数长度增加估商正确率的估商法,我也是在宝宝的理论和你的程序提示下才对估商法有了新的认识。
我的迭代法也是规约了除数的,和你的原理差不多,但没有规约被除数,我在计算倒数的过程中,有些关键环节也是在作整数运算,但我增加了逻辑上的约束,就是始终把运算结果当作小数看待,我有专门的指数位监视中间的小数点变化。进入最终运算环节的也是小数。
你的迭代法通过规约除数,被除数,把后续的运算控制在易于编程调试的范围。
我在编写高精表达式计算器中大量的用到字符串规范化,和高精数的规范化,和高精数范围缩小化,好处就是逻辑清晰(这个很重要,逻辑清晰在整个程序设计中都是非常重要的!),程序设计简单易调试。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2019-3-26 18:41 , Processed in 0.048919 second(s), 15 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表