大数除法之迭代法

liangbch

进入第一次循环，此时，参入运算的y的精度是40，x0的精度为20，第一次乘法是40位乘以20位，....
第二次乘法也是40位乘以20位，...

按照我的流程，第二次乘法，计算s2= s1 * x0，因为s1的高半部分为0，所以20位乘以20位,比你的方法节省一半的运算量。

落叶

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精度？

liangbch

我在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。

Ickiverar

我半年前对着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。

Ickiverar

这是$v_3+2^64-2^128/x$的散点图。其中横坐标是$x$，从$2^63$取到$2^64-1$。很明显，所有散点被$-9/8$和$0$这两个界限制住了，这表明我们的定理(4)是正确的。
事实上我证明了一个更紧的下界，就是那根紫色曲线，但那个曲线表达式太复杂，写出来要一屏幕，我就只给了个弱点的简单下界。
由于$x$可能的取值太多，不可能举尽，图上只大概取了十万个随机点描述一下基本图像。

落叶

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

因为我没有专攻数学，我对数学的应用是学懂一点原理，直接吸取别人的成果用于开发 ，我对程序开发或技术研究的理念是团队合作，科学管理，随着知识量的大爆炸，在科学研究上是没有所谓的超人的，我也没能力学到大而全。
其实科学上的很多原理，几大数学猜想，它们都没有得到充分的证明，但实际应用中或多或少把它们作为理论依据，包括我们初中所学的物理，化学的定理，最后发现有问题，但它们在一定范围可以应用于生活，或作为科学研究的参考指导。

落叶

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

落叶

四舍五入本身就是一种不精确操作，这里主要是探讨能否精确推算冗余精度，而且大数运算四舍五入会使误差更大。

zeroieme

落叶 发表于 2017-12-31 16:02
四舍五入本身就是一种不精确操作，这里主要是探讨能否精确推算冗余精度，而且大数运算四舍五入会使误差更大 ...

四舍五入与直接舍去都是不精确操作。四舍五入相当于引入-0.5~0.5的误差；直接舍去相当于引入0~1的误差。
用牛顿法会修补这一点误差。当然尽量减少误差是最好的。

落叶

四舍五入的综合误差是要远大于直接舍去，当你进行复杂运算时，对中间结果四舍五入简直没法算，当初我也是采取的四舍五入，我估计大多人碰到这种情况都会四舍五入，我也是。
经过大量的测试，我决定采用直接舍去。
当然这里面还有其它一些技术细节或理论，我就不多说了。我保留一点秘密。这可能就是理论和实践的区别，理论上分析四舍五入和直接舍入误差差不多，但事实却相反！

补充内容 (2018-1-1 08:07):
(四舍五入的综合误差是要远大于直接舍去)改为：四舍五入的综合误差在实际应用中，因为各种客观限制和操作关联误差要远大于直接舍去