Interfacing gmp and mpfr/mpc

[Andreas Enge <andreas.enge@inria.fr>]

With my PhD student, we wish to call library functions using floating point
numbers from the GNU MPFR and GNU MPC libraries inside Pari/GP.
Thanks to lots of help by Bill Allombert (who essentially dictated the
solution), we ended up with the attached code.

The functions mp{fr,c}_init_set_GEN take as argument an uninitialised
variable of type mp{fr,c}, allocate it on the pari stack and set its value
to that of the GEN given as second argument. The content of the variable
is copied, so they may be modified separately. Notice that unlike "normal"
mp{fr,c} variables, here they should not be initialised or cleared; they
may be modified, but their precision cannot be changed. If pari floating point
numbers are given as inputs, their pari precision (which is a multiple of
the bit length of long) is preserved; if integers are given, the current
pari default precision is used (thanks to additional parameters "long prec"
and the magic involved when installing C functions into the GP interpreter
with the type "p" of this variable).

Thanks to the GNU MPFR developers for providing the custom interface to make
this easily possible!

Since usually functions return mpfr and mpc results by side-effects, the
variables holding such results need to be allocated on the pari stack (instead
of by mp{fr,c}_init); this is done via the functions mp{fr,c}_pari_init.
As their argument prec is a pari precision (number of words of the signifi-
cand + 2), their mp{fr,c} precision is automatically a multiple of the word
length.

The functions mp{fr,c}_get_GEN return conversely a pari object of type GEN
from an mp{fr,c} number, the precision of which is assumed to be a multiple
of the word length (which is the case if they are obtained by the functions
explained above).

The function mul is given merely as an example of use. It takes two pari
floating point numbers and returns their product, computed by a call to
mpfr_mul or mpc_mul. To use it inside a GP session, execute the command
   install ("mul", "GGp", , "./mpfrc-pari.so");
The letters "GGp" stand for the types of the arguments (the return type GEN
being implicit, and the "p" being automatically replaced by the current
floating point precision, so that the function indeed takes only two
arguments). This assumes that the C file has been compiled into a shared
library by the command
   gcc mpfrc-pari.c -o mpfrc-pari.so -shared -fPIC -lmpc -lmpfr -lpari -I/usr/local/pari-dev/include -L/usr/local/pari-dev/lib
(with appropriately chosen include and library locations to the development
version of pari).

Comments and suggestions for improvement are welcome.

Andreas

#include <stdio.h>
#include <mpc.h>
#include <pari/pari.h>

/*
GP; install ("mul", "GGp", , "./mpfrc-pari.so");
*/

static void xmpn_mirrorcopy(mp_limb_t *z, mp_limb_t *x, long n)
{
  long i;
  for(i=0;i<n;i++)
    z[i]=x[n-1-i];
}


void mpfr_init_set_real (mpfr_ptr x, GEN r)
   /* creates space for x on the pari stack and sets it to r, which is
      supposed to be a t_REAL */
{
   mp_prec_t prec = realprec (r) - 2;
   GEN rc = new_chunk (prec);
   xmpn_mirrorcopy (rc, r+2, prec);
   mpfr_custom_init_set (x,
                         signe (r) ? signe (r) * MPFR_REGULAR_KIND : MPFR_ZERO_KIND,
                         expo (r) + 1,
                         BITS_IN_LONG * prec,
                         rc);
}


void mpfr_init_set_GEN (mpfr_ptr x, GEN r, long prec)
   /* creates space for x on the pari stack and sets it to r.
      If r is a t_REAL, its own precision is used;
      if r is a t_INT, prec is used as pari precision. */
{
   switch (typ (r)) {
      case t_REAL:
         mpfr_init_set_real (x, r);
         break;
      case t_INT:
         mpfr_init_set_real (x, itor (r, prec));
         break;
      default:
         pari_err_TYPE ("mpfr_init_set_GEN", r);
   }
}


void mpc_init_set_GEN (mpc_ptr z, GEN c, long prec)
   /* creates space for z on the pari stack and sets it to c.
      prec is passed through to mpfr_init_set_GEN. */
{
   if (typ (c) != t_COMPLEX) {
      mpfr_init_set_GEN (mpc_realref (z), c, prec);
      mpfr_init_set_GEN (mpc_imagref (z), gen_0, prec);
   }
   else {
      mpfr_init_set_GEN (mpc_realref (z), gel (c, 1), prec);
      mpfr_init_set_GEN (mpc_imagref (z), gel (c, 2), prec);
   }
}


/*
void
mpfr_print (GEN r, long prec)
{
   pari_sp av = avma;
   mpfr_t x;

   mpfr_init_set_GEN (x, r, prec);
   mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
   printf ("\n");

   avma = av;
}


void mpc_print (GEN c, long prec)
{
   pari_sp av = avma;
   mpc_t z;

   mpc_init_set_GEN (z, c, prec);
   mpc_out_str (stdout, 10, 0, z, MPC_RNDNN);
   printf ("\n");

   avma = av;
}
*/


void mpfr_pari_init (mpfr_ptr x, long prec)
   /* initialises x on the pari stack with the pari precision prec */
{
   mpfr_init_set_GEN (x, gen_0, prec);
}

void mpc_pari_init (mpc_ptr z, long prec)
   /* initialises z on the pari stack with the pari precision prec */
{
   mpc_init_set_GEN (z, gen_I (), prec);
}


GEN mpfr_get_GEN (mpfr_srcptr x)
   /* returns x as a t_REAL */
{
   long prec = mpfr_get_prec (x) / BITS_IN_LONG + 2;
   if (mpfr_zero_p (x))
      return real_0 (prec);
   else {
      GEN r = cgetr (prec);
      xmpn_mirrorcopy (r+2, mpfr_custom_get_significand (x), prec - 2);
      setsigne (r, (mpfr_signbit (x) ? -1 : 1));
      setexpo (r, mpfr_get_exp (x) - 1);
      return r;
   }
}


GEN mpc_get_GEN (mpc_srcptr z)
   /* returns z as a t_COMPLEX */
{
   GEN c = cgetg (3, t_COMPLEX);
   gel (c, 1) = mpfr_get_GEN (mpc_realref (z));
   gel (c, 2) = mpfr_get_GEN (mpc_imagref (z));
   return c;
}


GEN mul (GEN c1, GEN c2, long prec)
   /* test function */
{
   if (typ (c1) != t_COMPLEX && typ (c2) != t_COMPLEX) {
      mpfr_t z, z1, z2;
      mpfr_pari_init (z, prec);
      mpfr_init_set_GEN (z1, c1, prec);
      mpfr_init_set_GEN (z2, c2, prec);
      mpfr_mul (z, z1, z2, MPC_RNDNN);
      return mpfr_get_GEN (z);
   }
   else {
      mpc_t z, z1, z2;
      mpc_pari_init (z, prec);
      mpc_init_set_GEN (z1, c1, prec);
      mpc_init_set_GEN (z2, c2, prec);
      mpc_mul (z, z1, z2, MPC_RNDNN);
      return mpc_get_GEN (z);
   }
}