c - MPFR 之类的库如何存储和操作任意精度的实数？

Question

我目前正在开发一个应用程序，该应用程序需要使用任意大精度实数来计算pi到非常大的精度范围。我发现 MPFR 是一个了不起的库，高度优化和高效，适合我的目的。但是，我对如何实现任意精度实数非常好奇。我非常熟悉任意精度有符号整数的实现，并且我已经自己实现了一个C，但我不知道这些库如何非常有效地操作实数运算。如果我在处理整数，我可以存储数字mod 2**32作为堆上数组的元素，然后对基本的加法、减法、乘法、除法等进行传统的教科书数学方法。我认为为实数开发任意精度的实现可能是一个很好的挑战。我感谢每一个帮助我朝着正确方向前进的帮助。

标签： cfloating-pointimplementationarbitrary-precision

Answer 1

您可以直接在他们的 github 上检查库的实现：

特别是从此文件（第 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)

---

在这里你有乘法的源文件：

/* 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, l, 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));

以及这里关于算法的解释（来源：chtz）：

在MPFR 库：算法和证明以及不同的源文件中，所有内容都得到了很好的解释。