Add a toy implementation of complex ball arithmetic with relative error,wip-balls

using doubles for the ball radius. The test gives an example for which
small exponents are not enough.

6 files changed, 300 insertions, 3 deletions

diff --git a/configure.ac b/configure.ac
index a97485c..55b49d8 100644
--- a/configure.ac
+++ b/configure.ac
@@ -148,6 +148,7 @@ AC_CHECK_FUNCS([gettimeofday localeconv setlocale getrusage])
 AC_CHECK_FUNCS([dup dup2],,
         [AC_DEFINE([MPC_NO_STREAM_REDIRECTION],1,[Do not check mpc_out_str on stdout])])
 
+AC_CHECK_LIB([m], [fesetround])
 AC_CHECK_LIB([gmp], [__gmpz_init],
              [LIBS="-lgmp $LIBS"],
              [AC_MSG_ERROR([libgmp not found or uses a different ABI (including static vs shared).])])
diff --git a/src/Makefile.am b/src/Makefile.am
index f489219..4334ada 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -31,6 +31,7 @@ libmpc_la_SOURCES = mpc-impl.h abs.c acos.c acosh.c add.c add_fr.c         \
   urandom.c set.c 							   \
   set_prec.c set_str.c set_x.c set_x_x.c sin.c sin_cos.c sinh.c sqr.c      \
   sqrt.c strtoc.c sub.c sub_fr.c sub_ui.c sum.c swap.c tan.c tanh.c        \
-  uceil_log2.c ui_div.c ui_ui_sub.c
+  uceil_log2.c ui_div.c ui_ui_sub.c \
+  balls.c
 
 libmpc_la_LIBADD = @LTLIBOBJS@
diff --git a/src/balls.c b/src/balls.c
new file mode 100644
index 0000000..25da3a4
--- /dev/null
+++ b/src/balls.c
@@ -0,0 +1,157 @@
+/* balls -- Functions for complex ball arithmetic.
+
+Copyright (C) 2018 INRIA
+
+This file is part of GNU MPC.
+
+GNU MPC is free software; you can redistribute it and/or modify it under
+the terms of the GNU Lesser General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
+more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program. If not, see http://www.gnu.org/licenses/ .
+*/
+
+#include <stdio.h>
+#include <math.h>
+#include <fenv.h>
+#include "mpc-impl.h"
+
+void mpcb_print (mpcb_srcptr op)
+{
+   mpc_out_str (stdout, 10, 0, op->c, MPC_RNDNN);
+   printf (" %g\n", op->r);
+}
+
+
+void
+mpcb_init (mpcb_ptr rop)
+{
+   mpc_init2 (rop->c, 2);
+   rop->r = INFINITY;
+}
+
+
+void
+mpcb_clear (mpcb_ptr rop)
+{
+   mpc_clear (rop->c);
+}
+
+
+mpfr_prec_t
+mpcb_get_prec (mpcb_srcptr op)
+{
+   return mpc_get_prec (op->c);
+}
+
+
+void
+mpcb_set_prec (mpcb_ptr rop, mpfr_prec_t prec)
+{
+   mpc_set_prec (rop->c, prec);
+   rop->r = INFINITY;
+}
+
+
+void
+mpcb_set_c (mpcb_ptr rop, mpc_srcptr op)
+{
+   mpc_set_prec (rop->c, MPC_MAX_PREC (op));
+   mpc_set (rop->c, op, MPC_RNDNN);
+   rop->r = 0.0;
+}
+
+
+void
+mpcb_set (mpcb_ptr rop, mpcb_srcptr op)
+{
+   mpc_set_prec (rop->c, mpc_get_prec (op->c));
+   mpc_set (rop->c, op->c, MPC_RNDNN);
+   rop->r = op->r;
+}
+
+
+void
+mpcb_init_set_c (mpcb_ptr rop, mpc_srcptr op)
+{
+   mpc_init2 (rop->c, MPC_MAX_PREC (op));
+   mpc_set (rop->c, op, MPC_RNDNN);
+   rop->r = 0.0;
+}
+
+
+void
+mpcb_mul (mpcb_ptr z, mpcb_srcptr z1, mpcb_srcptr z2)
+/* FIXME: For the time being, we assume that z is different from z1 and from z2 */
+{
+   double r;
+   mpfr_prec_t p = mpc_get_prec (z1->c);
+
+   mpcb_set_prec (z, p);
+   mpc_mul (z->c, z1->c, z2->c, MPC_RNDNN);
+
+   fesetround (FE_UPWARD);
+   /* generic error of multiplication */
+   r = z1->r + z2->r + z1->r * z2->r;
+   /* error of rounding to nearest */
+   r += ldexp (1 + r, -p);
+   z->r = r;
+}
+
+
+void
+mpcb_add (mpcb_ptr z, mpcb_srcptr z1, mpcb_srcptr z2)
+/* FIXME: For the time being, we assume that z is different from z1 and from z2 */
+{
+   double r, r1, r2;
+   mpfr_prec_t p = mpc_get_prec (z1->c);
+
+   mpcb_set_prec (z, p);
+   mpc_add (z->c, z1->c, z2->c, MPC_RNDNN);
+
+   fesetround (FE_UPWARD);
+   /* absolute error of z1 and z2, divided by sqrt(2) */
+   r1 = ldexp (z1->r, MPC_MAX (mpfr_get_exp (mpc_realref (z1->c)),
+                               mpfr_get_exp (mpc_imagref (z1->c))));
+   r2 = ldexp (z2->r, MPC_MAX (mpfr_get_exp (mpc_realref (z2->c)),
+                               mpfr_get_exp (mpc_imagref (z2->c))));
+   /* absolute error of z, divided by sqrt(2) */
+   r = r1 + r2;
+   /* relative error of z; the sqrt(2), which come from the approximation
+      of complex norms by powers of 2, cancel out */
+   r = ldexp (r, - (MPC_MIN (mpfr_get_exp (mpc_realref (z->c)),
+                             mpfr_get_exp (mpc_imagref (z->c)))
+                    - 1));
+      /* The -1 comes from taking the lower bound of 1/2 for the mantissa
+         of c, whereas above we used the upper bound of 1 for r1 and r2.
+         So we may lose a factor of 2 here, for instance when z1==z2.
+         To prevent this, we would need to work with an approximation
+         of z, z1 and z2 in the type of r. */
+   /* error of rounding to nearest */
+   r += ldexp (1 + r, -p);
+   z->r = r;
+}
+
+void
+mpcb_sqrt (mpcb_ptr z, mpcb_srcptr z1)
+/* FIXME: Actually compute the error instead of assuming it is the same... */
+/* FIXME: For the time being, we assume that z is different from z1 */
+{
+   mpc_sqrt (z->c, z1->c, MPC_RNDNN);
+   z->r = z1->r;
+}
+
+void
+mpcb_div_2ui (mpcb_ptr z, mpcb_srcptr z1, unsigned long int e)
+{
+   mpc_div_2ui (z->c, z1->c, e, MPC_RNDNN);
+   z->r = z1->r;
+}
+
diff --git a/src/mpc.h b/src/mpc.h
index db09e8b..58e375d 100644
--- a/src/mpc.h
+++ b/src/mpc.h
@@ -108,6 +108,16 @@ typedef __mpc_struct mpc_t[1];
 typedef __mpc_struct *mpc_ptr;
 typedef const __mpc_struct *mpc_srcptr;
 
+typedef struct {
+  mpc_t  c;
+  double r;
+}
+__mpcb_struct;
+
+typedef __mpcb_struct mpcb_t [1];
+typedef __mpcb_struct *mpcb_ptr;
+typedef const __mpcb_struct *mpcb_srcptr;
+
 /* Support for WINDOWS DLL, see
    http://lists.gforge.inria.fr/pipermail/mpc-discuss/2011-November/000990.html;
    when building the DLL, export symbols, otherwise behave as GMP           */
@@ -239,6 +249,19 @@ __MPC_DECLSPEC int mpc_inp_str    (mpc_ptr, FILE *, size_t *, int, mpc_rnd_t);
 __MPC_DECLSPEC size_t mpc_out_str (FILE *, int, size_t, mpc_srcptr, mpc_rnd_t);
 #endif
 
+__MPC_DECLSPEC void mpcb_print (mpcb_srcptr);
+__MPC_DECLSPEC void mpcb_init (mpcb_ptr);
+__MPC_DECLSPEC void mpcb_clear (mpcb_ptr);
+__MPC_DECLSPEC mpfr_prec_t mpcb_get_prec (mpcb_srcptr);
+__MPC_DECLSPEC void mpcb_set_prec (mpcb_ptr, mpfr_prec_t);
+__MPC_DECLSPEC void mpcb_set (mpcb_ptr, mpcb_srcptr);
+__MPC_DECLSPEC void mpcb_set_c (mpcb_ptr, mpc_srcptr);
+__MPC_DECLSPEC void mpcb_init_set_c (mpcb_ptr, mpc_srcptr);
+__MPC_DECLSPEC void mpcb_mul (mpcb_ptr, mpcb_srcptr, mpcb_srcptr);
+__MPC_DECLSPEC void mpcb_add (mpcb_ptr, mpcb_srcptr, mpcb_srcptr);
+__MPC_DECLSPEC void mpcb_sqrt (mpcb_ptr, mpcb_srcptr);
+__MPC_DECLSPEC void mpcb_div_2ui (mpcb_ptr, mpcb_srcptr, unsigned long int);
+
 #if defined (__cplusplus)
 }
 #endif
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 2b7d447..49f133f 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -1,6 +1,6 @@
 ## tests/Makefile.am -- Process this file with automake to produce Makefile.in
 ##
-## Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013, 2016 INRIA
+## Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013, 2016, 2018 INRIA
 ##
 ## This file is part of GNU MPC.
 ##
@@ -37,7 +37,8 @@ check_PROGRAMS = tabs tacos tacosh tadd tadd_fr tadd_si tadd_ui targ	\
   tpow_d tpow_fr tpow_ld tpow_si tpow_ui tpow_z tprec tproj treal	\
   treimref trootofunity                                                 \
   tset tsin tsin_cos tsinh tsqr tsqrt tstrtoc tsub tsub_fr	        \
-  tsub_ui tsum tswap ttan ttanh tui_div tui_ui_sub tget_version exceptions
+  tsub_ui tsum tswap ttan ttanh tui_div tui_ui_sub tget_version exceptions \
+  tballs
 
 check_LTLIBRARIES=libmpc-tests.la
 libmpc_tests_la_SOURCES = mpc-tests.h check_data.c clear_parameters.c	\
diff --git a/tests/tballs.c b/tests/tballs.c
new file mode 100644
index 0000000..110c45a
--- /dev/null
+++ b/tests/tballs.c
@@ -0,0 +1,114 @@
+/* tballs -- test file for complex ball arithmetic.
+
+Copyright (C) 2018 INRIA
+
+This file is part of GNU MPC.
+
+GNU MPC is free software; you can redistribute it and/or modify it under
+the terms of the GNU Lesser General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
+more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program. If not, see http://www.gnu.org/licenses/ .
+*/
+
+#include "mpc-tests.h"
+
+int
+test_exp (void)
+{
+   mpfr_prec_t p;
+   mpc_t c;
+   mpfr_t one;
+   mpcb_t z, minusone, factor, tmp;
+   int i, n;
+
+   p = 20;
+   n = 4;
+   mpfr_init2 (one, p);
+   mpc_init2 (c, p);
+   mpc_set_si_si (c, -1, -1, MPC_RNDNN);
+   mpcb_init_set_c (minusone, c);
+   mpfr_set_ui (one, 1, MPFR_RNDN);
+   mpfr_exp (mpc_realref (c), one, MPFR_RNDU);
+   mpfr_exp (mpc_imagref (c), one, MPFR_RNDD);
+   mpcb_init_set_c (z, c);
+   mpcb_init (tmp);
+   mpcb_init (factor);
+
+   mpcb_print (z);
+   for (i = 1; i <= n; i++) {
+      mpcb_add (tmp, z, minusone);
+      printf (" +");
+      mpcb_print (tmp);
+      mpc_set_ui_ui (c, i, 0, MPC_RNDNN);
+      mpcb_set_c (factor, c);
+      mpcb_mul (z, tmp, factor);
+      printf ("%3i ", i);
+      mpcb_print (z);
+   }
+
+   mpfr_clear (one);
+   mpc_clear (c);
+   mpcb_clear (z);
+   mpcb_clear (minusone);
+   mpcb_clear (factor);
+   mpcb_clear (tmp);
+
+   return -1;
+}
+
+int
+test_agm (void)
+{
+   mpfr_prec_t p;
+   mpc_t c;
+   mpcb_t a, b, a1, b1;
+   int i, n;
+
+   p = 20;
+   n = 10;
+   mpc_init2 (c, p);
+   mpc_set_si_si (c, 1, 0, MPC_RNDNN);
+   mpcb_init_set_c (a, c);
+   mpc_set_si_si (c, 0, 1, MPC_RNDNN);
+   mpcb_init_set_c (b, c);
+   mpcb_init (a1);
+   mpcb_init (b1);
+
+   printf ("0 ");
+   mpcb_print (a);
+   printf ("  ");
+   mpcb_print (b);
+   for (i = 1; i <= n; i++) {
+      mpcb_add (a1, a, b);
+      mpcb_mul (b1, a, b);
+      mpcb_div_2ui (a, a1, 1);
+      mpcb_sqrt (b, b1);
+      printf ("%i ", i);
+      mpcb_print (a);
+      printf ("  ");
+      mpcb_print (b);
+   }
+
+   mpc_clear (c);
+   mpcb_clear (a);
+   mpcb_clear (b);
+   mpcb_clear (a1);
+   mpcb_clear (b1);
+
+   return -1;
+}
+
+int
+main (void)
+{
+   return test_agm ();
+}
+