误差的原因是进位产生的，进位的不确定性也决定了冗余精度的不确定性
这样的描述不准确。总体来说，物理中的误差是由测量引起的。数学中的误差是因为，绝大多数有理数的小数表示是无限循环小数的形式，所有的无理数都是无限不循环小数，故用小数的形式（不管是什么进制）来表示数，一定会有误差。
误差分为绝对误差和相对误差。比如，误差小于\( 2^{-100} \)是指绝对误差。又如，精确到4位有效数字，这是指相对误差。
确实不准确，当时想了半天，没找到合适的用语。应该说该例中出现的问题和进位与否有一定关系。编程中的精度控制是误差形成的一个主因。不同的精度控制会形成不同的误差特征，有些特征更利于编程开发。
大数运算中的精度控制是指相对误差。这种说法才是正确的。
冗余精度的做法广泛存在于 计算过程，和数的存储过程中。
先说数的存储中的冗余精度，如果进制为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.
一般的，我们只能尽可能使最终结果达到指定的精度，但完全做到很难。四则运算中，加，乘, 除这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 应该可以做到这一点。
@只是呼吸
我确实是那个qsd573。不光是qsd573，qsd572和qsd574也都是我。:lol
本帖最后由 只是呼吸 于 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位十进制数
当除数是一个单精度的 unsignedlong 时，调用其它函数处理，不进入 Newton_div()。
void Newton_div(UINT32 * const qut,UINT32 * const temp,UINT32 * const Work,UINT32 * const Cdiv,int LENdividend, UINT32 * const divR,int LENdivR)
{
UINT64 h64=100000000000000000;//17个 0
int len,i,s,N,LEN2,LenWork;
len=LENdividend-LENdivR+1;//计算商的有效位数len=LL-LR
N=(int)(floor(log(len)/log(2)));//计算迭代的次数 N=floor(LOG(Len)/LOG(2)); ;
Work=(UINT32)(h64/divR); //求初商，估商s0=100000000 00000000/divR;
Work/=10;
Work--;
LenWork=0;
len=1;
for(i=0;i<=N;i++)
{
ALU_mul(qut,divR,Work,LENdivR+1,LenWork+1);// 1 乘法函数：求a*s0
s=LENdivR+LenWork+1;
s=Newton_LeftZero(qut,s);// 寻找 qut[] 的非零最高位的下标 函数
Newton_sub(qut,s); // 2 减法函数，2-a*s0
ALU_mul(temp,qut,Work,s+2,LenWork+1); /// 乘法(2-a*s0)*s0
LEN2=s+LenWork+2;
LEN2=Newton_LeftZero(temp,LEN2);// 寻找 temp[] 的非零最高位的下标 函数
Newton_clear(Work,LenWork+2);//清零函数把 Work[] 清零
len<<=1; // 精度加倍
LenWork=len;// 计算迭代后的精度，也正好是Work 的长度
Newton_set(Work,temp,LEN2,LenWork);///赋值，运算关系Work<--temp, 赋的长度是LenWork, 裁剪掉temp 中的多余的右边的数(低位）。
Newton_Zero(qut,s+4); //清零函数把 qut[] 清零
Newton_Zero(temp,LEN2+4);//清零函数把 temp[] 清零
}///迭代的边界
len=LENdividend-LENdivR+1;
for(i=0;i<=len-1;i++)
{
Work=Work;//取Work左面的len位，赋给Work.
}
for(i=len;i<=LenWork;i++)///将多余的数清零
Work=0;
ALU_mul(qut,Cdiv,Work,LENdividend+1,len);///乘法，用倒数乘 被除数得到商。qut[] 存储商
s=Newton_LeftZero(qut,LENdividend+len+1);//寻找 qut[] 的非零最高位的下标 函数
len=quotient_Num(Cdiv,divR,LENdividend,LENdivR,qut);///求商的精确长度的函数，（不是把商求出来，只是求出商的精确长度）
for(i=0;i<=len-1;i++)
{
qut=qut;//取Work左面的s位，赋给Work.
}
for(i=len;i<=s;i++)///将多余的数清零
qut=0;
// return 0;
}//*///
乘法函数就是我写的那个“大整数乘法代码”，把它做了优化后暂时用着。
减法函数是专门针对牛顿除法的.代码如下：
static int Newton_sub(UINT32 *dst,int size2)///
{
int i,s=0;
while(dst==0)
{
s++;
}
dst=100000000-dst;
for(i=s+1;i<=size2;i++)
{
dst=99999999-dst;
}
dst=1;//(///
return 0;
}////////
其余的函数都没有什么特殊的，就不写出来了。
这个牛顿除法终于完工，花了很多时间，一开始的时候，一个数字都不对。后来才慢慢正确。
实际运算：35800位十进制数除以 11000位十进制数，我写的牛顿除法用时35毫秒。我写的试商除法用时 110毫秒。
可能是乘法不好，这么大的数字用硬乘法也许勉强了。等下一步写一个karatsuba乘法试试。
本帖最后由 落叶 于 2018-1-13 14:07 编辑
只是呼吸 发表于 2018-1-13 00:56
以下是我自己写的牛顿迭代除法函数。
除法函数的入口条件：
要添加 头文件 ...
请教一下！
如果把除数123456789代入
Work=(UINT32)(h64/divR); //求初商，估商s0=100000000 00000000/divR;
Work/=10;
Work--;
的具体数为：
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: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\) 行的那个。
你的迭代法确实很不错，可是说绝大部分都是调试找规律的结晶，凭想象凭原理都是很难实现的，另外也很佩服你在没有理论（宝宝关于估商法的理论）的支持下，加长除数长度增加估商正确率的估商法，我也是在宝宝的理论和你的程序提示下才对估商法有了新的认识。
我的迭代法也是规约了除数的，和你的原理差不多，但没有规约被除数，我在计算倒数的过程中，有些关键环节也是在作整数运算，但我增加了逻辑上的约束，就是始终把运算结果当作小数看待，我有专门的指数位监视中间的小数点变化。进入最终运算环节的也是小数。
你的迭代法通过规约除数，被除数，把后续的运算控制在易于编程调试的范围。
我在编写高精表达式计算器中大量的用到字符串规范化，和高精数的规范化，和高精数范围缩小化，好处就是逻辑清晰（这个很重要，逻辑清晰在整个程序设计中都是非常重要的！），程序设计简单易调试。