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

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

[复制链接]
发表于 2017-12-27 22:57:33 | 显示全部楼层
进入第一次循环,此时,参入运算的y的精度是40,x0的精度为20,第一次乘法是40位乘以20位,....
第二次乘法也是40位乘以20位,...

按照我的流程,第二次乘法,计算s2= s1 * x0,因为s1的高半部分为0,所以20位乘以20位,比你的方法节省一半的运算量。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-12-28 14:09:27 | 显示全部楼层
liangbch 发表于 2017-12-27 22:57
按照我的流程,第二次乘法,计算s2= s1 * x0,因为s1的高半部分为0,所以20位乘以20位,比你的方法节省一 ...

根据你在一篇文章中说到,要想得到N位的精度结果,参入运算的操作数精度要大于等于N位,在迭代循环中,以我的程序为例来分析,在第一个循环中,y的精度是40,x0的精度是20,根据上面的理论分析,是会出现计算结果达不到精度40的情况,我在实际的测试中也验证了这种情况,我的解决方法是增加冗余精度,在第一个循环中我的设定精度是大于40的,相应的y精度也大于40,x0的精度是大于20的,冗余精度的长度我用的土方法,一点点试出来的。不知道有什么方法可以推出这个冗余精度?我的循环中40精度乘以20精度数会出现计算结果精度达不到40的情况,不知你程序中(按照我的流程,第二次乘法,计算s2= s1 * x0,因为s1的高半部分为0,所以20位乘以20位,比你的方法节省一半的运算量),这段代码所在循环是求20的精度还是求40精度?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-12-28 14:41:57 | 显示全部楼层
我在5楼说过
如果你在每次迭代时,源操作精度仅仅保持目标精度的一半,你会发现,最终精度是达不到要求的。

请看5楼的代码,我增加了2个额外的比特。这是靠广泛的测试得到的最佳值。测试时,要保证足够的覆盖率,足够的测试次数。
t_prc.b += EXTRA_BITS;                // enlarge target precision to guarantee to get enough precision

步骤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. 在步骤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


以x0的精度是20位为例,再详细解释下上面的步骤。
1. a取高40位,x0取20位,结果大约是60位,丢弃低20位,仅保留40位。
2. 一般情况下,s1的高20位为0,所以s1仅包含20位有效数字。
3. 做s2= s1 * x0,结果为40位,丢弃低20位,仅仅保留高20位。
4. 将x0作精度扩展,扩展至40位,低20位为0.做x1=x0+s2或者x1=x0-s2。注意,低位对齐。换言之,s2的高20位为0。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-12-30 03:49:55 | 显示全部楼层
本帖最后由 Ickiverar 于 2017-12-30 04:14 编辑

我半年前对着GMP的源码把
1.查表法初步求逆(9bit 初始精度)
2.牛顿迭代法求逆(64bit,单个字长数)
3.牛顿迭代法求高精逆(多个字长的高精度数)
4.除法的估商法
5.除法的乘逆求商法
这些东西的误差全部从数学上严格地证明了一遍,可以做到精度控制到最后1bit或者最后一个字。
然后我觉得,做这种东西,如果想要代码很优雅,就必须从数学上证明它的正确性,而不是仅凭经验给出一个大概的精度。
比如我举个例子(虽然例子用$2^64$做基,但原理是可以适用于任意基的):
对一个64bit-整数$2^63<=x<2^64$,如何计算$v=\floor{\frac{2^128-1}x}-2^64$?
根据$v$的定义很容易发现一定有$0<v<2^64$,也就是$v$可以用一个uint64无损表示,这也是为什么我们在$v$的定义式后面减一个$2^64$,抗溢出啊。
这个其实就是计算除数的最高64位$x$的逆$v$。这种定义使得$x(v+2^64)$尽可能接近但严格小于$2^128$,这个性质在后面各种求高精逆、求商的误差分析里还有很大用处。
然后有如下算法:
(1).根据$x$的最高9bit(即$\floor{x/2^55}$)查表,这个表是事先计算的,有$2^9=512$个元素,可以直接查出$v_0=\round{2^19/(1/2+\floor{x/2^55})}$,round是四舍五入。
(2).计算$v_1=v_0*2^10-\floor{(v_0^2\floor{x/2^24})/2^41}$
(3).计算$v_2=v_1*2^15+\floor{(v_1*(2^59-v_1*\floor{x/2^24}))/2^44}$
(4).计算$v_3=v_2*2^30+\floor{(v_2*(2^96-\floor{(v_2*x)/4}-1))/2^66}-2^64$
$v_3$就是对真正的$v$的一个非常好的估计,可以严格证明(留作习题,较难):$-9/8<=v_3+2^64-2^128/x<=0$
(5).由(4)中结论可以证明(留作习题,较易):$v=(v_3-\floor{((v_3+2^64+1)x)/2^64}) mod 2^64$.
算法里所有的除法的除数都是2的幂,而且紧跟着一个floor,这其实就是位运算的右移,实际上根本不用做除法。
以上只用了整数乘法和位运算就计算出了$x$的64bit精确逆$v$,汇编很容易实现。经严格证明,结果一丝不差。

这只是1.2.的过程,楼主应该对3.的过程更感兴趣,不过分析的思路都是一样的,只是精度的单位从bit变成了64bit-word。

点评

写错了,是“数论变换吧”  发表于 2018-1-6 20:48
想问 Ickiverar 一个问题( 如涉及个人隐私,可以不回答):你是不是“数论变化吧”的那个 qsd573?  发表于 2018-1-6 15:56
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-12-30 04:14:05 | 显示全部楼层

这是$v_3+2^64-2^128/x$的散点图。其中横坐标是$x$,从$2^63$取到$2^64-1$。很明显,所有散点被$-9/8$和$0$这两个界限制住了,这表明我们的定理(4)是正确的。
事实上我证明了一个更紧的下界,就是那根紫色曲线,但那个曲线表达式太复杂,写出来要一屏幕,我就只给了个弱点的简单下界。
由于$x$可能的取值太多,不可能举尽,图上只大概取了十万个随机点描述一下基本图像。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-12-30 19:53:54 | 显示全部楼层
本帖最后由 落叶 于 2017-12-30 20:11 编辑

冗余精度的推算,我数学太差,没有能力从数学上精确推算,只能根据精度原理(宝宝解析的国外文献),和一些代表性数据测试进行模糊推理,然后根据结果指导程序的后期开发。我认同你的观念,也希望做到更好。

因为我没有专攻数学,我对数学的应用是学懂一点原理,直接吸取别人的成果用于开发 ,我对程序开发或技术研究的理念是团队合作,科学管理,随着知识量的大爆炸,在科学研究上是没有所谓的超人的,我也没能力学到大而全。
其实科学上的很多原理,几大数学猜想,它们都没有得到充分的证明,但实际应用中或多或少把它们作为理论依据,包括我们初中所学的物理,化学的定理,最后发现有问题,但它们在一定范围可以应用于生活,或作为科学研究的参考指导。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-12-31 04:45:52 | 显示全部楼层
本帖最后由 落叶 于 2017-12-31 04:48 编辑

我又重新思考了一下冗余精度,在做加法时,所设精度之外保留多少个基本运算单位,才能保证计算的完全正确呢?我采用的是一个基本运算单位,1234.4444444444........5+1234.5555555555.........5应等于2469,但当我们设定精度为有效数四位时,冗余精度设更长(1或大于1),最终的结果都是2468,误差的原因是进位产生的,进位的不确定性也决定了冗余精度的不确定性,但实际运算采用完全精度是不可能或不实用的,所以只能折中处理了。

点评

断尾应当四舍五入。  发表于 2017-12-31 15:41
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-12-31 16:02:12 | 显示全部楼层
四舍五入本身就是一种不精确操作,这里主要是探讨能否精确推算冗余精度,而且大数运算四舍五入会使误差更大。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-12-31 16:12:58 | 显示全部楼层
落叶 发表于 2017-12-31 16:02
四舍五入本身就是一种不精确操作,这里主要是探讨能否精确推算冗余精度,而且大数运算四舍五入会使误差更大 ...

四舍五入与直接舍去都是不精确操作。四舍五入相当于引入-0.5~0.5的误差;直接舍去相当于引入0~1的误差。
用牛顿法会修补这一点误差。当然尽量减少误差是最好的。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-12-31 17:13:46 | 显示全部楼层
四舍五入的综合误差是要远大于直接舍去,当你进行复杂运算时,对中间结果四舍五入简直没法算,当初我也是采取的四舍五入,我估计大多人碰到这种情况都会四舍五入,我也是。
经过大量的测试,我决定采用直接舍去。
当然这里面还有其它一些技术细节或理论,我就不多说了。我保留一点秘密。这可能就是理论和实践的区别,理论上分析四舍五入和直接舍入误差差不多,但事实却相反!

补充内容 (2018-1-1 08:07):
(四舍五入的综合误差是要远大于直接舍去)改为:四舍五入的综合误差在实际应用中,因为各种客观限制和操作关联误差要远大于直接舍去
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-3-29 01:00 , Processed in 0.045236 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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