Precisions and rounding modes

In the first part of this mini-series we have gone through the motions of creating bindings from Guile to the floating point C library MPFR. As a running example, we have used IEEE 754 standard 113 bit quadruple precision. But the real goal of MPFR is to handle arbitrary (and in particular, arbitrarily large) precisions. Additionally, following the IEEE standard, it provides for different rounding modes, while so far we have limited ourselves to (admittedly the most common) rounding to nearest. This instalment treats these two aspects in a way that is idiomatic to Guile, provides flexibility and leaves programmers using guile-mpfr in control.

Choices of precisions

Precision handling can be achieved by setting a global default precision, as we did so far; after modifying the global parameter, all subsequent variables that are initialised by a call to mpfr_init will have this new precision. In a functional programming context such a recourse to state is not very idiomatic; more down to earth, it also poses questions in multithreaded code. In fact each MPFR variable carries its own precision, and the function mpfr_init2 takes as additional argument the precision at which a variable is initialised. This makes the precision explicit and should generally be preferred (the MPC library, written after MPFR, made the decision to completely drop the concept of a global precision and to only provide the function mpc_init2).

In MPFR, the variable that is supposed to hold the result of a computation is passed as an argument to the C function; in this way, also the desired target precision is passed along. In the functional context (but also when operators are overloaded, as in C++), results are instead returned as new objects, and choosing an appropriate precision becomes an issue with several possible solutions: One could use the global default precision regardless of the precision of the input values. One could use the maximum precision among all inputs with the hope of minimising precision losses. One could use the minimum precision among all inputs; implicitly this acknowledges that floating point inputs are already tainted by numerical imprecisions, that the last bit, say, may already be wrong, and that this error spills over into the result. One could explicitly consider all inputs as intervals with a given error (in the last bit, say), and use a precision that optimally approximates a result with the same semantics.

The choice made in guile-mpfr is to refuse any choice and to let the programmer remain responsible of precision handling as in the C library.

Handling precisions in the C wrapper library

In a first step, we modify guile_mpfr.c such that each C function returning an MPFR value takes as additional parameter the target precision. In C this would have the mpfr_prec_t type (usually a long); since we want to integrate the functions into the Guile REPL, the argument needs to be of SCM type, and we have to add an explicit conversion. So we add a new function scm_to_prec and modify the initialisation function mpfr_malloc to take an additional precision parameter, which is used to call mpfr_init2 instead of mpfr_init:

static mpfr_prec_t scm_to_prec (SCM prec)
{
   return (mpfr_prec_t) scm_to_long (prec);
}

static SCM mpfr_malloc (SCM prec)
{
   mpfr_ptr x;

   x = (mpfr_ptr) malloc (sizeof (__mpfr_struct));
   mpfr_init2 (x, scm_to_prec (prec));

   return scm_from_mpfrptr (x);
}

After that the conversion functions returning an MPFR value, as well as the macros wrapping computational functions, handle an additional parameter for the desired precision of the return value, which is passed through to mpfr_malloc; for instance:

SCM integer_to_mpfr (SCM val, SCM prec)
{
   SCM res;
   mpfr_ptr x;
   mpz_t z;

   res = mpfr_malloc (prec);
   x = scm_to_mpfrptr (res);
   mpz_init (z);
   scm_to_mpz (val, z);
   mpfr_set_z (x, z, MPFR_RNDN);
   mpz_clear (z);

   return res;
}

#define GUILE_MPFR_F_FF(name) \
SCM scm_mpfr_ ## name (SCM val1, SCM val2, SCM prec) \
{ \
   SCM res; \
   mpfr_ptr x1, x2, x; \
   res = mpfr_malloc (prec); \
   x = scm_to_mpfrptr (res); \
   x1 = scm_to_mpfrptr (val1); \
   x2 = scm_to_mpfrptr (val2); \
   mpfr_ ## name (x, x1, x2, MPFR_RNDN); \
   return res; \
}

For all these functions, we modify the library initialisation routine init_guile_mpfr by incrementing the number of arguments recorded through calls to scm_c_define_gsubr. In preparation of the following section, we also mark them as internal by prepending their names with an underscore:

void init_guile_mpfr (void)
{
   init_mpfr_type ();
   scm_c_define_gsubr ("_integer->mpfr", 2, 0, 0, integer_to_mpfr);
   scm_c_define_gsubr ("mpfr->double", 1, 0, 0, mpfr_to_double);
   scm_c_define_gsubr ("_mpfr-sin", 2, 0, 0, scm_mpfr_sin);
   scm_c_define_gsubr ("_mpfr-mul", 3, 0, 0, scm_mpfr_mul);
   scm_c_define_gsubr ("_mpfr-div", 3, 0, 0, scm_mpfr_div);
   scm_c_define_gsubr ("_mpfr-pi", 1, 0, 0, scm_mpfr_pi);
   scm_c_define_gsubr ("mpfr-print", 1, 0, 0, scm_mpfr_print);
}

The names of the exported functions at the Guile side are adapted accordingly. Again to prepare the following section, we add an indirection on the module level and record the changes in a new module mpfr-internal in the file mpfr-internal.scm:

(define-module (mpfr-internal)
  #:export (_integer->mpfr
            mpfr->double
            mpfr-print
            _mpfr-sin
            _mpfr-mul
            _mpfr-div
            _mpfr-pi))

(load-extension "guile_mpfr" "init_guile_mpfr")

The reason for the indirection is that load-extension should be the last line in a module. The loaded procedures can be used after that, but since they are loaded dynamically, the module compiler does not know their names and emits a warning about potentially unbound variables. Once the symbols are exported, however, they can be used without raising a warning. Since we do not want to export the symbols with an underscore to the outside world, we need to export them from an internal module that is supposed to not be imported into user code.

We can now do the same computations as in the previous post by starting a Guile REPL with

make
export GUILE_EXTENSIONS_PATH=.:$GUILE_EXTENSIONS_PATH
guile -L .

and executing

(use-modules (mpfr-internal))
(mpfr-print (_mpfr-mul (_integer->mpfr 2 53) (_mpfr-sin (_mpfr-div (_mpfr-pi 53) (_integer->mpfr 3 53) 53) 53) 53))
(mpfr-print (_mpfr-mul (_integer->mpfr 2 113) (_mpfr-sin (_mpfr-div (_mpfr-pi 113) (_integer->mpfr 3 113) 113) 113) 113))

to obtain the square root of 3 once in double and then in quadruple precision; the only difference to the previous code is that to each procedure call we add the precision at which it is carried out. The current state of the code is recorded in commit cf6da40 of the guile-mpfr project.

Parameterising the precision

Flexibility we have achieved, but convenience we have definitely lost. In a typical computation one may wish to change the precision from time to time (for instance when realising that the current precision is too low to obtain a satisfactory result), but usually, as in the example at the end of the previous section, whole swathes of computation are done at the same precision, and having to repeat the precision for each operation is cumbersome.

To solve this problem, we do something that at first sight looks like going back to the previously frowned upon solution of introducing a global variable; a variable, however, that will allow for local modification, namely a Guile parameter. On the surface a parameter is just a variable; more precisely it is a procedure that serves as a container for a value: If called without an argument, it returns the current value; if called with an argument, it stores the argument as the new current value. We add such a parameter to the (mpfr) module:

(define prec (make-parameter 53))

setting it to the default MPFR default precision. Then we can redefine all the procedures returning an MPFR value to take the current value of prec into account; for instance:

(define (mpfr-mul x y)
  (_mpfr-mul x y (prec)))

so that the (exported) procedure mpfr-mul calls the (hidden inside the mpfr-internal module) procedure _mpfr-mul with its own arguments, to which the current value of the precision parameter is added. This results in the following file mpfr.scm:

(define-module (mpfr)
  #:use-module (mpfr-internal)
  #:re-export (mpfr->double
               mpfr-print)
  #:export (integer->mpfr
            mpfr-sin
            mpfr-mul
            mpfr-div
            mpfr-pi
            prec))

(define prec (make-parameter 53))

(define (integer->mpfr x)
  (_integer->mpfr x (prec)))

(define (mpfr-sin x)
  (_mpfr-sin x (prec)))

(define (mpfr-mul x y)
  (_mpfr-mul x y (prec)))

(define (mpfr-div x y)
  (_mpfr-div x y (prec)))

(define (mpfr-pi)
  (_mpfr-pi (prec)))

Notice how the internal procedures are first imported from (mpfr-internal); the procedures not requiring a precision are then passed through by re-exporting them, while all others are exported under their final name. With this code, which is commit fef0e0c, the previous example becomes:

(use-modules (mpfr))
(mpfr-print (mpfr-mul (integer->mpfr 2) (mpfr-sin (mpfr-div (mpfr-pi) (integer->mpfr 3)))))
(prec 113)
(mpfr-print (mpfr-mul (integer->mpfr 2) (mpfr-sin (mpfr-div (mpfr-pi) (integer->mpfr 3)))))

Excursion: Guile macros

Again the transformation of underscored internal functions to exported functions follows the same pattern over and over again and depends only on the arity of the function, which calls for the use of macros. There are different solutions for this problem, and I thank Sergio Pastor Pérez and Ludovic Courtès for discussions and suggestions.

Coming from the C preprocessor, the following solution may feel natural. As in the C case with ##, we wish to programmatically concatenate tokens, respectively symbols in the Guile context, and fixed strings such as "_". This can be done conveniently by going back and forth between symbols and strings and doing the transformation on the strings. For instance, using define-macro we can write a macro for functions with two MPFR inputs and use it to define mpfr-mul and mpfr-div as follows:

(define-macro (private->public-2 f)
  (let ((private-f (string->symbol (string-append "_" (symbol->string f)))))
    `(define-public (,f x y)
      (,private-f x y (prec)))))

(private->public-2 mpfr-mul)
(private->public-2 mpfr-div)

The quasi-quoted define-public expression creates a nested list of symbols; it starts with 'define-public, then contains a list of three symbols, the macro argument f as a symbol such as 'mpfr_mul (here unquoting is needed since otherwise we would end up with the symbol 'f), 'x and 'y, and so on. Otherwise said, we end up with a quoted S-expression (or sexp) that holds the Scheme code we would otherwise write by hand. The calls to this macro are then evaluated in the current environment and create new public procedures mpfr-mul and mpfr-div. Since we have used define-public instead of define the procedures are also automatically exported from the module and can be dropped from the #:export clause of the module definition. Adding macros for the functions of other arities, we end up with the following complete module in mpfr.scm:

(define-module (mpfr)
  #:use-module (mpfr-internal)
  #:re-export (mpfr->double
               mpfr-print)
  #:export (prec))

(define prec (make-parameter 53))

(define-macro (private->public-2 f)
  (let ((private-f (string->symbol (string-append "_" (symbol->string f)))))
    `(define-public (,f x y)
      (,private-f x y (prec)))))

(define-macro (private->public-1 f)
  (let ((private-f (string->symbol (string-append "_" (symbol->string f)))))
    `(define-public (,f x)
      (,private-f x (prec)))))

(define-macro (private->public-0 f)
  (let ((private-f (string->symbol (string-append "_" (symbol->string f)))))
    `(define-public (,f)
      (,private-f (prec)))))

(private->public-1 integer->mpfr)
(private->public-2 mpfr-mul)
(private->public-2 mpfr-div)
(private->public-1 mpfr-sin)
(private->public-0 mpfr-pi)

Except for quasi-quoting and unquoting, which allows for a bit more flexibility (the computation with symbols part), this is essentially syntactic source code rewriting like for the C preprocessor. Said otherwise, this is very close to text file handling using a script in an editor, combinations of sed and awk on the command line, or with the Autotools.

Scheme purists, however, will worry about hygiene, although I think that the previous macros are a reasonable approach in our well circumscribed situation: They are not exported outside the module, in which they are used with very little context, so that referential transparency should not be an issue. If one wants to be in style, one should instead use define-syntax together with syntax-rules, or the define-syntax-rules shortcut. The previous approach of creating one symbol from another is not hygienically possible. So one needs to pass both the names of the existing hidden and the new public procedure as arguments, to obtain the following macro for a procedure with two arguments used to define multiplication:

(define-syntax-rule (private->public-2 private-f public-f)
  (define-public
    (public-f x y)
      (private-f x y (prec))))

(private->public-2 _mpfr-mul mpfr-mul)

Indeed the macro definition looks cleaner and less twisted, but there is a bit more boiler plate typing for using it. The resulting file mpfr.scm remains as compact, and the module can be used to run the same computations as previously:

(define-module (mpfr)
  #:use-module (mpfr-internal)
  #:re-export (mpfr->double
               mpfr-print)
  #:export (prec))

(define prec (make-parameter 53))

(define-syntax-rule (private->public-2 private-f public-f)
  (define-public
    (public-f x y)
      (private-f x y (prec))))

(define-syntax-rule (private->public-1 private-f public-f)
  (define-public
    (public-f x)
      (private-f x (prec))))

(define-syntax-rule (private->public-0 private-f public-f)
  (define-public
    (public-f)
      (private-f (prec))))

(private->public-1 _integer->mpfr integer->mpfr)
(private->public-2 _mpfr-mul mpfr-mul)
(private->public-2 _mpfr-div mpfr-div)
(private->public-1 _mpfr-sin mpfr-sin)
(private->public-0 _mpfr-pi mpfr-pi)

Think local

The parameterize construction of Guile brings us back to the functional idiom. Much like let, which creates local bindings, it takes a list of parameters and values to which the parameters are dynamically bound, so that these values are used in all computations internal to the envrionment. Since we already took care to export the prec parameter from the (mpfr) module, we can try it out immediately in the Guile REPL:

(use-modules (mpfr))
(parameterize ((prec 113))
  (mpfr-print (mpfr-mul (integer->mpfr 2) (mpfr-sin (mpfr-div (mpfr-pi) (integer->mpfr 3))))))
(mpfr-print (mpfr-mul (integer->mpfr 2) (mpfr-sin (mpfr-div (mpfr-pi) (integer->mpfr 3)))))

Notice how the first series of computations is carried out at local quadruple precision, whereas the value of prec remains globally unchanged, so that the second series is carried out with the default double precision.

To make things even more comfortable, we may hide the internal prec parameter and renounce at exporting it from the module (after all, its exact name does not matter), and instead export a bit of syntax:

(define-syntax-rule (with-prec p exp ...)
   (parameterize ((prec p))
     exp ...))

Then altogether, the previous example can be rewritten as

(use-modules (mpfr))
(define (sqrt3) (mpfr-mul (integer->mpfr 2) (mpfr-sin (mpfr-div (mpfr-pi) (integer->mpfr 3)))))
(with-prec 113
  (mpfr-print (sqrt3)))
(mpfr-print (sqrt3))

to compute an expression at two different precisions. (I have the vague impression that with-prec creates a monad; it is left as an exercise to the more knowledgeable reader to prove or disprove this statement.) The resulting code can be found in commit bae58c8.

Rounding modes

MPFR currently provides six rounding modes. When a result is exactly representable in the target precision, then is is returned as such. Otherwise the exact result falls in between two representable values, and it is rounded to one of them depending on the rounding mode. While in the MPFR C library the precision of a result is implicit and determined by the precision with which the variable holding it has been initialised, the rounding mode is chosen explicitly at each function invocation by passing one of six constants as an additional argument. Paradoxically, the easiest to understand rounding modes are the directed ones, which always round into the same direction, potentially depending on the sign of the number: MPFR_RNDU for rounding up (or towards +∞); MPFR_RNDD for rounding down (or towards -∞); MPFR_RNDZ for rounding towards 0 (or down for positive, up for negative values); MPFR_RNDA for rounding away towards infinity (or up for positive, down for negative numbers). Clearly these two come in pairs: If MPFR_RNDU returns one of the possibilities of a non-representable number, then MPFR_RNDD returns the other one, and the same holds for MPFR_RNDZ and MPFR_RNDA. The most conventional rounding mode, MPFR_RNDN, rounds to nearest, but suffers from an ambiguity when a number is exactly in the middle between two representable numbers, which are then both at the same distance. The strategy chosen by MPFR (and in fact the IEEE 754 standard) is to resolve ties to even, that is, to round towards the binary number with a 0 in the last position. Unlike its name suggests, MPFR_RNDF, or faithful rounding, is rather fickle: It non-deterministically returns one of the two possibilities, but without guaranteeing which one. So it is potentially faster at the price of losing reproducibility.

In guile-mpfr we export all these values to Guile as numerical constants. Rounding modes themselves can then be handled in an analogous way to the precision, by adding a hidden parameter (defaulting to MPFR_RNDN) and additional syntax with-rnd, which locally binds the rounding mode, and with-prec-rnd, which locally binds the precision and the rounding mode. Then the following code computes the square root of 3 at different precisions and rounding modes (it returns successively smaller values).

(use-modules (mpfr))
(define (sqrt3) (mpfr-mul (integer->mpfr 2) (mpfr-sin (mpfr-div (mpfr-pi) (integer->mpfr 3)))))
(mpfr-print (with-prec 113 (sqrt3)))
(mpfr-print (with-prec-rnd 113 MPFR_RNDD (sqrt3)))
(mpfr-print (with-rnd MPFR_RNDD (sqrt3)))

This finishes the mini-series on guile-mpfr; the current state of the code is given by commit b5cea0f. The scaffolding is in place, what is missing now is completing the functionality by wrapping all of MPFR!