mpc-impl.h, fma.c: export function fma_naive for tests

tfma.c: add random tests comparing fma to fma_naive
fma.dat: test with exact 0 as result


git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@1168 211d60ee-9f03-0410-a15a-8952a2c7a4e4

5 files changed, 77 insertions, 4 deletions

diff --git a/configure.ac b/configure.ac
index a6b0d6a..907a168 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1,4 +1,4 @@
-# Copyright (C) 2008, 2009, 2010, 2011 INRIA
+# Copyright (C) 2008, 2009, 2010, 2011, 2012 INRIA
 #
 # This file is part of GNU MPC.
 #
diff --git a/src/fma.c b/src/fma.c
index 197cd65..b38f2be 100644
--- a/src/fma.c
+++ b/src/fma.c
@@ -1,6 +1,6 @@
 /* mpc_fma -- Fused multiply-add of three complex numbers
 
-Copyright (C) 2011 INRIA
+Copyright (C) 2011, 2012 INRIA
 
 This file is part of GNU MPC.
 
@@ -39,7 +39,7 @@ bound_prec_addsub (mpfr_srcptr x, mpfr_srcptr y)
 }
 
 /* r <- a*b+c */
-static int
+int
 mpc_fma_naive (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
 {
   mpfr_t rea_reb, rea_imb, ima_reb, ima_imb, tmp;
diff --git a/src/mpc-impl.h b/src/mpc-impl.h
index a9f9906..c3ca9af 100644
--- a/src/mpc-impl.h
+++ b/src/mpc-impl.h
@@ -157,6 +157,7 @@ extern "C" {
 
 __MPC_DECLSPEC int  mpc_mul_naive (mpc_ptr, mpc_srcptr, mpc_srcptr, mpc_rnd_t);
 __MPC_DECLSPEC int  mpc_mul_karatsuba (mpc_ptr, mpc_srcptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_fma_naive (mpc_ptr, mpc_srcptr, mpc_srcptr, mpc_srcptr, mpc_rnd_t);
 __MPC_DECLSPEC int  mpc_pow_usi (mpc_ptr, mpc_srcptr, unsigned long, int, mpc_rnd_t);
 __MPC_DECLSPEC char* mpc_alloc_str (size_t);
 __MPC_DECLSPEC char* mpc_realloc_str (char*, size_t, size_t);
diff --git a/tests/fma.dat b/tests/fma.dat
index 756cd2f..49e54c6 100644
--- a/tests/fma.dat
+++ b/tests/fma.dat
@@ -25,3 +25,4 @@
 # see sin.dat for precisions
 
 + - 53 -0x189281b52abc03p-55 53 0x178a1d3cd134e5p-49 53 0x16A09E667F3BCDp-52 53 0x9CC470A049097p-50 53 0x23C6EF372FE95p-48 53 0x9CC470A049097p-50 53 0xA953FD4E97C75p-50 53 0x16A09E667F3BCDp-51 N N
+0 0 10 0 10 0  10 7 10 5 10 3 10 -12 10 -81 10 69 N N
diff --git a/tests/tfma.c b/tests/tfma.c
index 79430d7..02eed66 100644
--- a/tests/tfma.c
+++ b/tests/tfma.c
@@ -1,6 +1,6 @@
 /* tfma -- test file for mpc_fma.
 
-Copyright (C) 2011 INRIA
+Copyright (C) 2011, 2012 INRIA
 
 This file is part of GNU MPC.
 
@@ -20,6 +20,75 @@ along with this program. If not, see http://www.gnu.org/licenses/ .
 
 #include "mpc-tests.h"
 
+static void
+cmpfma (mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
+   /* computes a*b+c with the naive and fast functions using the rounding
+      mode rnd and compares the results and return values.                                          
+      In our current test suite, all input precisions are the same, and we
+      use this precision also for the result.
+   */
+{
+   mpc_t z, t;
+   int   inex_z, inex_t;
+
+   mpc_init2 (z, MPC_MAX_PREC (a));
+   mpc_init2 (t, MPC_MAX_PREC (a));
+
+   inex_z = mpc_fma_naive (z, a, b, c, rnd);
+   inex_t = mpc_fma (t, a, b, c, rnd);
+
+   if (mpc_cmp (z, t) != 0 || inex_z != inex_t) {
+      fprintf (stderr, "fma_naive and fma differ for rnd=(%s,%s)\n",
+               mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+               mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+      MPC_OUT (a);
+      MPC_OUT (b);
+      MPC_OUT (c);
+      MPC_OUT (z);
+      MPC_OUT (t);
+      if (inex_z != inex_t) {
+         fprintf (stderr, "inex_re (z): %s\n", MPC_INEX_STR (inex_z));
+         fprintf (stderr, "inex_re (t): %s\n", MPC_INEX_STR (inex_t));
+      }
+      exit (1);
+   }
+
+  mpc_clear (z);
+  mpc_clear (t);
+}
+
+
+static void
+check_random (void)
+{
+   mpfr_prec_t prec;
+   mpfr_rnd_t rnd_re, rnd_im;
+   mpc_t a, b, c;
+
+   mpc_init2 (a, 1000);
+   mpc_init2 (b, 1000);
+   mpc_init2 (c, 1000);
+
+   for (prec = 2; prec < 1000; prec = (mpfr_prec_t) (prec * 1.1 + 1)) {
+      mpc_set_prec (a, prec);
+      mpc_set_prec (b, prec);
+      mpc_set_prec (c, prec);
+
+      test_default_random (a, -1024, 1024, 128, 0);
+      test_default_random (b, -1024, 1024, 128, 0);
+      test_default_random (c, -1024, 1024, 128, 0);
+
+      for (rnd_re = 0; rnd_re < 4; rnd_re ++)
+         for (rnd_im = 0; rnd_im < 4; rnd_im ++)
+            cmpfma (a, b, c, RNDC(rnd_re, rnd_im));
+   }
+
+   mpc_clear (a);
+   mpc_clear (b);
+   mpc_clear (c);
+}
+
+
 int
 main (void)
 {
@@ -27,6 +96,8 @@ main (void)
 
   test_start ();
 
+  check_random ();
+
   data_check (f, "fma.dat");
   tgeneric (f, 2, 1024, 1, 256);