问题描述
我目前正在开发一个需要使用任意大精度实数来计算 pi 到非常大的精度范围的应用程序。我发现 MPFR 是一个令人惊叹的库，高度优化且高效，符合我的目的。但是，我对如何实现任意精度实数非常好奇。我非常熟悉任意精度有符号整数的实现，并且我已经在 C 中自己实现了一个，但我不知道了解这些库如何非常有效地操作实数运算。
如果我处理整数，我可以将数字 mod 2**32 作为数组的元素存储在堆上，然后使用传统的教科书数学方法进行基本的加法、减法、乘法、除法等。
我认为为实数开发任意精度的实现可能是一个很好的挑战。我感谢每一个帮助我朝着正确方向前进的帮助。
解决方法
您可以检查库 directly on their github 的实现：
特别是来自 this file（第 910 行）：
#define HAVE_DECIMAL128_IEEE 1
unsigned int t3:32;
unsigned int t2:32;
unsigned int t1:32;
unsigned int t0:14;
unsigned int comb:17;
unsigned int sig:1;
所以库使用 128 位来存储浮点信息，并使用精度、指数、咒语和肢体的组合来计算浮点之间的运算：
#define MPFR_PREC(x) ((x)->_mpfr_prec)
#define MPFR_EXP(x) ((x)->_mpfr_exp)
#define MPFR_MANT(x) ((x)->_mpfr_d)
#define MPFR_GET_PREC(x) MPFR_PREC_IN_RANGE (MPFR_PREC (x))
#define MPFR_LAST_LIMB(x) ((MPFR_GET_PREC (x) - 1) / GMP_NUMB_BITS)
#define MPFR_LIMB_SIZE(x) (MPFR_LAST_LIMB (x) + 1)
并且 here 你有乘法的源文件：
/* now the full product is {h,l,v + w + high(b0*c0),low(b0*c0)},where the lower part contributes to less than 3 ulps to {h,l} */
/* If h has its most significant bit set and the low sh-1 bits of l are not
000...000 nor 111...111 nor 111...110,then we can round correctly;
if h has zero as most significant bit,we have to shift left h and l,thus if the low sh-2 bits are not 000...000 nor 111...111 nor 111...110,then we can round correctly. To avoid an extra test we consider the latter
case (if we can round,we can also round in the former case).
For sh <= 3,we have mask <= 7,thus (mask>>2) <= 1,and the approximation
cannot be enough. */
494 if (MPFR_LIKELY(((l + 2) & (mask >> 2)) > 2))
495 sb = sb2 = 1; /* result cannot be exact in that case */
496 else
497 {
498 umul_ppmm (sb,sb2,bp[0],cp[0]);
499 /* the full product is {h,sb + v + w,sb2} */
500 sb += v;
501 l += (sb < v);
502 h += (l == 0) && (sb < v);
503 sb += w;
504 l += (sb < w);
505 h += (l == 0) && (sb < w);
506 }
507 if (h < MPFR_LIMB_HIGHBIT)
508 {
509 ax --;
510 h = (h << 1) | (l >> (GMP_NUMB_BITS - 1));
511 l = (l << 1) | (sb >> (GMP_NUMB_BITS - 1));
512 sb <<= 1;
513 /* no need to shift sb2 since we only want to know if it is zero or not */
514 }
515 ap[1] = h;
516 rb = l & (MPFR_LIMB_ONE << (sh - 1));
517 sb |= ((l & mask) ^ rb) | sb2;
518 ap[0] = l & ~mask;
519
520 MPFR_SIGN(a) = MPFR_MULT_SIGN (MPFR_SIGN (b),MPFR_SIGN (c));
关于算法here的解释（来源：chtz）：
一切都在 THE MPFR LIBRARY: ALGORITHMS AND PROOFS 和不同的源文件中得到了很好的解释。