C, Guile, the REPL and everything in-between

The logical next step following the series on Goblins for number theory is to import software that is useful to number theory into Guile. This means creating bindings to C libraries, in particular libraries handling numbers with arbitrarily many digits. Since Guile uses GMP for handling big integers, there are already C functions transforming between both integer types. This would be enough to create the bindings needed to call my primality proving software from Guile. But I would like to start with a more clear-cut example for learning purposes. With the goal of providing an outcome of independent usefulness I have decided to start creating bindings to MPFR. This could also provide an occasion to look into the numerical tower of Guile.

Integrating C libraries into Guile

The first term that comes to mind when handling from within one language software written in another language is foreign function interface. Guile has mechanisms for opening dynamic libraries and calling the functions in it. This is all nice and svelte when these functions operate on scalar C data types; these have counterparts in the Scheme world. Additionally there is limited support for C structs, but the available functions look dangerously close to bean, erh, byte counting instead of providing a higher level interface. Apart from that, being reduced to what is called wrapped pointers, a euphemism for the infamous void * in C, the interface is not really type aware.

The (maybe more modern?) alternative to this opaque handling of wrapped pointers are foreign object types, which make it possible to define separate types for arbitrarily complex C data structures. These can even be equipped with mechanisms to make them compatible with the garbage collector of Guile. The price to pay is that we have to write a separate C library for providing the interface (much like the Java Native Interface, for instance, as used in PariDroid, or the approach taken in PariTwine to make C libraries available to the PARI/GP computer algebra system). Then the functionality from the foreign function interface comes in handy to dynamically load this intermediate library.

Defining MPFR numbers as a foreign object type

Writing the helper library is essentially a matter of following the Guile documentation for creating foreign object types and objects, and copy-pasting and adapting relevant boilerplate. So let us create a file guile_mpfr.c for the library code. We start by putting the following lines into the file:

#include <stdio.h>
#include <mpfr.h>
#include <libguile.h>

static SCM mpfr_type;

The first line will be invaluable for printf debugging; it needs to come first to enable the input and output functions of MPFR. The second one is necessary because eventually we will use bits and pieces of MPFR, the third one because we already use Guile. And the last one creates a variable that is local to the library and will hold our new type mpfr_type as a C variable of type SCM. As a caveat, notice that the word type is used with two different meanings here. On one hand, in C everything Guile is handled by the unique type SCM (for a reader coming from number theory, this will be reminiscent of the GEN generic C type in the PARI library); this is essentially an array of bytes or words. On the other hand, Guile (and PARI) have their own internal type systems; these are handled through the contents of the SCM (or GEN) values (by reserving initial bytes to hold the type information).

Let us start by doing things backwards and write the functions that transform between Scheme values of the still to be defined new type mpfr_type and MPFR values; the latter are given by C pointers of type mpfr_ptr. With this knowledge the transformation functions are essentially boilerplate:

SCM scm_from_mpfrptr (mpfr_ptr x)
{
   return scm_make_foreign_object_1 (mpfr_type, x);
}

mpfr_ptr scm_to_mpfrptr (SCM val)
{
   scm_assert_foreign_object_type (mpfr_type, val);
   return (mpfr_ptr) scm_foreign_object_ref (val, 0);
}

The call to scm_make_foreign_object_1 takes as input our MPFR variable and wraps and tags it as an mpfr_type. Conversely scm_foreign_object_ref unwraps the MPFR variable as a void *, which we convert to the correct type before returning it. The call to scm_assert_foreign_object_type is a safety measure to throw an error whenever we call the function with an SCM value that is not of the (internal) type mpfr_type.

Next let us write the code that creates and destroys MPFR variables as foreign objects. For good C hygiene, we make sure to write the two functions involving malloc and free in lockstep (although, strictly speaking, we currently only need the code for freeing variables to pass it to the Guile garbage collector).

static SCM mpfr_malloc ()
{
   mpfr_ptr x;

   x = (mpfr_ptr) malloc (sizeof (__mpfr_struct));
   mpfr_init (x);

   return scm_from_mpfrptr (x);
}

void mpfr_free (SCM val)
{
   mpfr_ptr x;

   x = scm_to_mpfrptr (val);
   mpfr_clear (x);
   free (x);
}

The mpfr_ptr type is in fact a pointer to an internal struct __mpfr_struct. In the mpfr_malloc function, we malloc such a data structure, then initialise the actual MPFR variable by a call to mpfr_init (which executes an additional call to malloc for an internal field of the struct), then wrap the resulting MPFR variable into a Scheme value. The mpfr_free function exactly reverses these operations.

Finally we are ready to define the MPFR object type and to assign it to the mpfr_type variable.

static void init_mpfr_type (void)
{
   SCM name, slots;
   scm_t_struct_finalize finaliser;

   name = scm_from_utf8_symbol ("mpfr");
   slots = scm_list_1 (scm_from_utf8_symbol ("data"));
   finaliser = (scm_t_struct_finalize) mpfr_free;
   mpfr_type = scm_make_foreign_object_type (name, slots, finaliser);

   mpfr_set_default_prec (113);
}

The heart of the function is the call to scm_make_foreign_object_type, which takes three arguments. The name is an arbitrary symbol (here, the string "mpfr" is transformed into the symbol 'mpfr; if these internals are confusing, you may think of it as an identifyer given by a string). The slots can handle an arbitrary number (as a list) of symbols that determine data slots; here we need only one slot to hold an MPFR value. We attach the symbol 'data to it, which is arbitrary since we actually handle this unique slot by its number in scm_from_mpfrptr and scm_to_mpfrptr. And finaliser is a C function that is called when the Guile garbage collector tries to delete a foreign object; this function needs to free all memory that has been allocated outside the control of Guile when creating the object (if there is no such memory, one may use the NULL pointer instead of a finaliser function).

Remember that MPFR handles real floating point numbers of arbitrary precision. Each variable can hold its own precision, but there is also a library-wide default value. We take the opportunity in init_mpfr_type to set this value to 113 bits, that is, quadruple precision. So the call to mpfr_init in mpfr_malloc will use this default precision (a more functional, stateless approach would be to pass the desired precision at each variable initialisation as a call mpfr_init2 (x, prec), see the second instalment of this miniseries).

Now the helper library containing all the above code in this order can be compiled via a Makefile containing the lines

guile_mpfr.so : guile_mpfr.c
        gcc -shared -o guile_mpfr.so -fPIC `pkg-config --cflags guile-3.0` `pkg-config --libs guile-3.0` -lmpfr -lgmp guile_mpfr.c

(be careful to use the tabulator key and not spaces at the beginning of the second line). Then running make creates the dynamic library file guile_mpfr.so, using flags directly extracted from the pkg-config file of Guile (in particular on Guix, this adds -I and -L command line arguments pointing directly to the correct subdirectory in /gnu/store/).

Conversion functions

Before being able to compute in Guile with MPFR variables, we need to connect the newly created MPFR island to the existing Guile continent by defining bridges, that is, functions that transform values between them. So we add the following to guile_mpfr.c:

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

   res = mpfr_malloc ();
   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;
}

For this to work, we need Guile to be compiled with GMP support. The default guile package in Guix is currently compiled with mini-gmp, an embedded version of GMP, so that the function scm_to_mpz does not get defined. For the full experience, one needs to install the guile-with-gmp package.

The structure of the function is straightforward. We allocate a new Scheme variable res of type mpfr_type by a call to mpfr_malloc, which will hold the result to be returned by the function. To be able to use MPFR functionality, we define the pointer x as a handle for the wrapped MPFR variable. Then we call the Guile function scm_to_mpz that transforms a Guile integer into a GMP integer (abstracting away from the question whether the integer is an immediate, that is, immediately stored in a machine word, or a pointer to an array of words; and hopefully causing an error when called with inappropriate arguments). The MPFR function mpfr_set_z then assigns this GMP integer to x. Every MPFR function takes an additional argument, the rounding mode; here and in the following, we use MPFR_RNDN, which means that all values that cannot be represented exactly at our target precision of 113 bits are rounded to the nearest representable value. With an integer this can happen if it is particularly large in absolute value.

In the converse direction, let us for the moment content ourselves with transforming an MPFR value to a Scheme value containing a double.

SCM mpfr_to_double (SCM val)
{
   double d;

   d = mpfr_get_d (scm_to_mpfrptr (val), MPFR_RNDN);

   return scm_from_double (d);
}

The code should be self-explaining; again we (usually) need to round when reducing from the 113 bit MPFR precision to the 53 bit double precision.

For good measure, we add one more function that prints an MPFR value. It uses the MPFR function mpfr_out_str to directly print to stdout, using the decimal system (10) and as many digits as necessary to uniquely specify the floating point number at 113 bits of precision (indicated by the argument 0). This will need refinement to handle Guile output ports, but for the time being it comes handy for debugging purposes, and to see if and what our upcoming functions compute with MPFR values.

void scm_mpfr_print (SCM x)
{
   mpfr_out_str (stdout, 10, 0, scm_to_mpfrptr (x), MPFR_RNDN);
   printf ("\n");
}

Arithmetic functions

Having conversion functions is enough to freely mix C code that works with Scheme objects and MPFR functions. However, it is our goal to eventually make the MPFR functions accessible from the Guile REPL. To this purpose, we need to write for every MPFR function an equivalent C function operating on foreign objects of type mpfr_type, which is essentially a mechanical exercise: Create a foreign object for the result, access the mpfr_ptr pointers behind each foreign object and call the MPFR function on the pointers. The code is very similar to what we have seen above for conversions. For example, let us add multiplication and division:

SCM scm_mpfr_mul (SCM val1, SCM val2)
{
   SCM res;
   mpfr_ptr x1, x2, x;
   res = mpfr_malloc ();
   x = scm_to_mpfrptr (res);
   x1 = scm_to_mpfrptr (val1);
   x2 = scm_to_mpfrptr (val2);
   mpfr_mul (x, x1, x2, MPFR_RNDN);
   return res;
}

SCM scm_mpfr_div (SCM val1, SCM val2)
{
   SCM res;
   mpfr_ptr x1, x2, x;
   res = mpfr_malloc ();
   x = scm_to_mpfrptr (res);
   x1 = scm_to_mpfrptr (val1);
   x2 = scm_to_mpfrptr (val2);
   mpfr_div (x, x1, x2, MPFR_RNDN);
   return res;
}

Automation

The excitement of writing wrapper functions by copy-pasting wears off very fast. Luckily we can use GCC macros to write the main content of the functions once and for all for a given prototype. For instance, the function scm_mpfr_div can be obtained from scm_mpfr_mul by replacing the two occurrences of mul by div. Since this happens in the middle of identifiers (scm_mpfr_mul itself and the function name mpfr_mul), we need a special construct, token concatenation through the ## preprocessing operator. The resulting code then becomes self-explaining; to keep the visual formatting, one just needs to also concatenate the different lines by ending them with a \.

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

GUILE_MPFR_F_FF(mul)
GUILE_MPFR_F_FF(div)

The idea behind the macro name GUILE_MPFR_F_FF is to encode the prototype of the MPFR function: It has one output and two inputs of MPFR type. Using analogous macros we complete the example by adding the sine function and the constant π as follows:

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

#define GUILE_MPFR_F(name) \
SCM scm_mpfr_ ## name (void) \
{ \
   SCM res; \
   mpfr_ptr x; \
   res = mpfr_malloc (); \
   x = scm_to_mpfrptr (res); \
   mpfr_const_ ## name (x, MPFR_RNDN); \
   return res; \
}

GUILE_MPFR_F_F(sin)
GUILE_MPFR_F(pi)

Registering C functions

Despite all the talk about bridging between C and Scheme and the use of the keyword SCM hinting at Scheme, so far we have remained firmly with C. It is time to register the C functions as Scheme variables. This is done by the scm_c_define_gsubr C function. For instance,

scm_c_define_gsubr ("mpfr-mul", 2, 0, 0, scm_mpfr_mul);

takes the scm_mpfr_mul C function and imports it into the Scheme REPL under the name mpfr-mul; it also specifies that the function takes two mandatory arguments (and no optional and remaining arguments).

We tie all the initialisation code (creation of the type and exporting of the C functions to Scheme) together into one function in guile-mpfr.c, which will serve as the entry point of the library:

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

Using the REPL

Running make we now obtain the compiled dynamic library guile-mpfr.so. It can be loaded using the foreign function interface as a foreign extension into a running Guile REPL by the Guile procedure load-extension. We go one step further and create a Guile module. For this, it is enough to place the following Guile code into the file mpfr.scm in the current directory:

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

(load-extension "guile_mpfr" "init_guile_mpfr")

The load-extension procedure opens the dynamic library and initialises it by calling the init_guile_mpfr function. As seen above, this function defines a number of Scheme procedures which, so far, in C parlance, are static to the module; but they are in fact #:exported through the module definition. To use the module, we run

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

in the terminal. The first line tells Guile to look for extension libraries in the current directory; the second line is there because it is better to call make too often than not often enough; the last line starts a Guile REPL and tells it to look for modules in the current directory.

In the REPL, we can now run the following code:

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

which should print the square root of 3 twice, once in quadruple precision, once (doubly) rounded to double precision.

Time to breathe before continuing with the next topic! The current state of the files is recorded as commit 486c622 of the guile-mpfr project in the multiprecision group on Codeberg.

Excursion: GOOPS!

Now that we have quadruple precision numbers, it is tempting to somehow insert them into the Guile numerical tower. Numbers are a special class of Guile objects in the GOOPS, the Guile Object Oriented Programming System, which works through polymorphic functions. So in addition to mpfr_type (which, somehow confusingly, is actually a C variable), we define an object type holding MPFR values in the file mpfr-goops.scm.

(define-module (mpfr-goops)
  #:use-module (mpfr)
  #:use-module (oop goops)
  #:export (<mpfr>
            mpfr-ref
            mul))

(define-class <mpfr> (<number>)
  (val #:init-value (integer->mpfr 0) #:init-keyword #:val #:accessor mpfr-ref))

Here we define the <mpfr> class as a subclass of the more general <number> class, itself a direct descendent of the most abstract class <top>. We provide it with one slot, val (again the name does not matter since we will access its content differently), which is supposed to hold one element of type mpfr_type. If an <mpfr> object is created without any further parameter, then our choice of #:init-value sets it to 0. But #:init-keyword adds the named parameter #:val to the creation procedure to choose a different initial value. Finally #:accessor adds a procedure with the name mpfr-ref to get and set the value. The following code snippet from the Guile REPL illustrates how to create, modify and print such objects.

(use-modules (mpfr) (oop goops) (mpfr-goops))
(define x (make <mpfr>))
(mpfr-print (mpfr-ref x))
(set! (mpfr-ref x) (integer->mpfr 1))
(mpfr-print (mpfr-ref x))
(define y (make <mpfr> #:val (integer->mpfr 2)))
(mpfr-print (mpfr-ref y))

First we create a variable x containing the common initial value 0, then replace this value by 1. The variable y is directly created holding the integral value 2.

We can now add a polymorphic multiplication function into the mpfr-goops module, which follows precisely the new-+ example from the Guile documentation:

(define-method (mul)
  (make <mpfr> #:val (integer->mpfr 1)))

(define-method (mul (x <mpfr>))
  (make <mpfr> #:val (mpfr-ref x)))

(define-method (mul (x <integer>))
  (make <mpfr> #:val (integer->mpfr x)))

(define-method (mul (x <mpfr>) . args)
  (make <mpfr>
    #:val (mpfr-mul (mpfr-ref x) (mpfr-ref (apply mul args)))))

(define-method (mul (x <integer>) . args)
  (mul (mul x) (apply mul args)))

It recursively defines a procedure mul that computes the product of its arguments, an arbitrary number of integers and <mpfr> objects in any order, by folding over the arguments. First, the product of no arguments is set to the neutral element 1. Then the product of one argument is set to its <mpfr> value; if it is an integer, one needs to add an explicit type cast. Finally the product of two or more arguments is obtained by multiplying the first argument with the product of all other ones.

We can now print the factorial of 4 by running

(mpfr-print (mpfr-ref (mul 3 y 4)))

in the above Guile REPL session.

We are not going to pursue this approach further; given that there are <integer>, <real> and <complex> numbers (and even though the last one is of no interest for real arithmetic, but would be relevant to Guile bindings for MPC), there is a combinatorial explosion for the different argument types. As is common when programming in C with the MPFR library, one can always add an explicit type conversion before calling a function, and then in fact the GOOPS-free approach of the previous section is sufficient. Nevertheless, the mpfr-goops module is available for illustrative purposes in commit 8121125 of the guile-mpfr project.

This post has covered the general steps of creating bindings from Guile to a C library, first by creating wrappers at the layer of C programming, which is also the language in which the foundations of Guile are written, and then on the level of calling the wrapped functions from a Guile REPL. The next post will cover facets that are particular to the MPFR library, namely the handling of floating point precisions and rounding modes. Nevertheless, it will also be an occasion to learn more Scheme in the process.