How does libraries like MPFR store and manipulate real numbers of arbitrary precision?

Question

I am currently working on an application that requires, the usage of arbitrary large precision real numbers to calculate pi to a very large precision range. I've found MPFR an amazing library highly optimized and efficient that suits my purpose. But, I'm extremely curious as to how arbitrary precision real numbers are implemented. I am extremely familiar with the implementation of arbitrary precision signed integers, and I've already implemented one by myself in C , but I've no idea as to how these libraries manipulate real number operations very efficiently. If I was dealing with integers, I could store the digits mod 2**32 as elements of an array on the heap, and then do traditional school book mathematical methods for basic addition, subtraction, multiplication, division, etc. I think developing an arbitrary precision implementation for real numbers, could be a good challenge. I appreciate every help that will nudge me in the right direction.

It is a huge challenge to develop such thing, see the number of lines of code they have in their library, arbitrary precision for signed integers is a joke in comparison, nevertheless it is indeed an interesting project! — Antonin GAVREL, Mar 25 '21 at 18:44
This is quite well documented here: https://www.mpfr.org/algorithms.pdf. I think for stackoverflow the question is too broad. — chtz, Mar 25 '21 at 18:55
You do not need floating point for computing Pi ... not even bigints are needed see [BBP](https://stackoverflow.com/a/56035284/2521214) also this might be interesting [Baking-Pi Challenge - Understanding & Improving](https://stackoverflow.com/a/22295383/2521214) for you... For arbitrary precision numbers you first need to decide if you want floating point or fixed point. However the core stuff is still just bigint — Spektre, Mar 29 '21 at 07:39

Antonin GAVREL · Answer 1 · 2021-03-25T19:51:00.107

You can check implementation of the library directly on their github:

Especially from this file (line 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;

So the library use 128 bits to store the floating point informatiom and use a combination of precision, exponent, mantra and limb to compute operations between floating points:

#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)

And here you have the source file for multiplications:

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

And explanation about the algorithm here (credit: chtz):

Everything is nicely explained in the both the THE MPFR LIBRARY: ALGORITHMS AND PROOFS and in the different source files.

How does libraries like MPFR store and manipulate real numbers of arbitrary precision?

1 Answers1