GMPFRXX: An integrated C++ Implementation of the GNU MP and MPFR Multiple-precision Libraries

Many important problems in scientific computing are extremely ill-conditioned, and require intermediate calculations to be done with many more digits of precision than the desired accuracy of the answer. Examples include resolving clusters of eigenvalues, designing numerical quadrature schemes for arbitrary families of functions, and computing asymptotic expansions of solutions of elliptic partial differential equations near corners (see the figure)).

In such problems, it is desirable to write a single code in which the floating point type can be easily switched between double precision for speed, or arbitrary precision for accuracy. MPFR (developed at INRIA) is an arbitrary precision floating point library in which arithmetic, elementary functions, and many special functions are provided with guaranteed correct rounding. We found the native C implementation of this library to be difficult to use as even simple algebraic expressions such as x=a*b+c*d have to be expanded out into separate function calls. Our code uses operator overloading and templates to hide these function calls from the user, dramatically simplifying the interface to the MPFR library. Our code also allows mixing of built-in types (int, double), arbitrary precision integer and rational types from the GMP library, and arbitrary precision floating point types from MPFR in a single algebraic expression, where each type is promoted as necessary to the more general type.

GMPFRXX is one of the official C++ interfaces of the MPFR project, and the only one that also supports GMP. Since its release in November, 2008, the webpage for GMPFRXX has been viewed from 2286 unique IP addresses and the code has been downloaded from 660 unique IP addresses, excluding web crawlers.

The computation of singular solutions of the equations of elasticity near a corner or interface junction involves locating zeros of the characteristic determinant of an analytic matrix function. Such zeros often come in clusters corresponding to nearly degenerate basis functions that must be separated before they can be used in a finite element computation. This procedure involves many steps in which double-precision arithmetic can lead to unacceptable loss of accuracy.