forgot two files

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@9793 280ebfd0-de03-0410-8827-d642c229c3f4

2 files changed, 242 insertions, 0 deletions

diff --git a/src/fmma.c b/src/fmma.c
new file mode 100644
index 000000000..de156ed9b
--- /dev/null
+++ b/src/fmma.c
@@ -0,0 +1,135 @@
+/* mpfr_fmma, mpfr_fmms -- Compute a*b +/- c*d
+
+Copyright (C) 2002-2016 INRIA.
+Copied from the GNU MPC project.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library 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.
+
+The GNU MPFR Library 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 the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
+http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
+51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
+
+#include "mpfr-impl.h"
+
+static int
+mpfr_fmma_slow (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
+                   mpfr_srcptr d, mpfr_rnd_t rnd)
+{
+  mpfr_t ab, cd;
+  int inex;
+
+  mpfr_init2 (ab, MPFR_PREC(a) + MPFR_PREC(b));
+  mpfr_init2 (cd, MPFR_PREC(c) + MPFR_PREC(d));
+  mpfr_mul (ab, a, b, MPFR_RNDZ); /* exact */
+  mpfr_mul (cd, c, d, MPFR_RNDZ); /* exact */
+  inex = mpfr_add (z, ab, cd, rnd);
+  mpfr_clear (ab);
+  mpfr_clear (cd);
+  return inex;
+}
+
+/* z <- a*b + c*d */
+static int
+mpfr_fmma_fast (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
+                mpfr_srcptr d, mpfr_rnd_t rnd)
+{
+   /* Assumes that a, b, c, d are finite and non-zero; so any multiplication
+      of two of them yielding an infinity is an overflow, and a
+      multiplication yielding 0 is an underflow.
+      Assumes further that z is distinct from a, b, c, d. */
+
+   int inex;
+   mpfr_t u, v;
+   mp_size_t an, bn, cn, dn;
+   mp_ptr up, vp;
+   MPFR_TMP_DECL(marker);
+   MPFR_SAVE_EXPO_DECL (expo);
+
+   MPFR_TMP_MARK(marker);
+   MPFR_SAVE_EXPO_MARK (expo);
+
+   /* u=a*b, v=c*d exactly */
+   an = MPFR_LIMB_SIZE(a);
+   bn = MPFR_LIMB_SIZE(b);
+   cn = MPFR_LIMB_SIZE(c);
+   dn = MPFR_LIMB_SIZE(d);
+   MPFR_TMP_INIT (up, u, (an + bn) * GMP_NUMB_BITS, an + bn);
+   MPFR_TMP_INIT (vp, v, (cn + dn) * GMP_NUMB_BITS, cn + dn);
+
+   /* u <- a*b */
+   if (an >= bn)
+     mpn_mul (up, MPFR_MANT(a), an, MPFR_MANT(b), bn);
+   else
+     mpn_mul (up, MPFR_MANT(b), bn, MPFR_MANT(a), an);
+   if ((up[an + bn - 1] & MPFR_LIMB_HIGHBIT) == 0)
+     {
+       mpn_lshift (up, up, an + bn, 1);
+       MPFR_EXP(u) = MPFR_EXP(a) + MPFR_EXP(b) - 1;
+     }
+   else
+     MPFR_EXP(u) = MPFR_EXP(a) + MPFR_EXP(b);
+
+   /* v <- c*d */
+   if (cn >= dn)
+     mpn_mul (vp, MPFR_MANT(c), cn, MPFR_MANT(d), dn);
+   else
+     mpn_mul (vp, MPFR_MANT(d), dn, MPFR_MANT(c), cn);
+   if ((vp[cn + dn - 1] & MPFR_LIMB_HIGHBIT) == 0)
+     {
+       mpn_lshift (vp, vp, cn + dn, 1);
+       MPFR_EXP(v) = MPFR_EXP(c) + MPFR_EXP(d) - 1;
+     }
+   else
+     MPFR_EXP(v) = MPFR_EXP(c) + MPFR_EXP(d);
+
+   MPFR_PREC(u) = (an + bn) * GMP_NUMB_BITS;
+   MPFR_PREC(v) = (cn + dn) * GMP_NUMB_BITS;
+   MPFR_SIGN(u) = MPFR_MULT_SIGN(MPFR_SIGN(a), MPFR_SIGN(b));
+   MPFR_SIGN(v) = MPFR_MULT_SIGN(MPFR_SIGN(c), MPFR_SIGN(d));
+
+   /* tentatively compute z as u+v; here we need z to be distinct
+      from a, b, c, d to avoid losing the input values in case we
+      need to call mpfr_fmma_slow */
+   inex = mpfr_add (z, u, v, rnd);
+
+   MPFR_TMP_FREE(marker);
+   MPFR_SAVE_EXPO_FREE (expo);
+
+   return mpfr_check_range (z, inex, rnd);
+}
+
+/* z <- a*b + c*d */
+int
+mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
+           mpfr_srcptr d, mpfr_rnd_t rnd)
+{
+  mp_ptr zp = MPFR_MANT(z);
+
+  return (mpfr_regular_p (a) && mpfr_regular_p (b) && mpfr_regular_p (c) &&
+          mpfr_regular_p (d) && zp != MPFR_MANT(a) && zp != MPFR_MANT(b) &&
+          zp != MPFR_MANT(c) && zp != MPFR_MANT(d))
+    ? mpfr_fmma_fast (z, a, b, c, d, rnd)
+    : mpfr_fmma_slow (z, a, b, c, d, rnd);
+}
+
+/* z <- a*b - c*d */
+int
+mpfr_fmms (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
+           mpfr_srcptr d, mpfr_rnd_t rnd)
+{
+  mpfr_t minus_c;
+
+  MPFR_ALIAS (minus_c, c, -MPFR_SIGN(c), MPFR_EXP(c));
+  return mpfr_fmma (z, a, b, minus_c, d, rnd);
+}
diff --git a/tests/tfmma.c b/tests/tfmma.c
new file mode 100644
index 000000000..651d145ff
--- /dev/null
+++ b/tests/tfmma.c
@@ -0,0 +1,107 @@
+/* Test file for mpfr_fmma and mpfr_fmms.
+
+Copyright 2016 Free Software Foundation, Inc.
+Contributed by the AriC and Caramel projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library 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.
+
+The GNU MPFR Library 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 the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
+http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
+51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
+
+#include "mpfr-test.h"
+
+/* check both mpfr_fmma and mpfr_fmms */
+static void
+random_test (mpfr_t a, mpfr_t b, mpfr_t c, mpfr_t d, mpfr_rnd_t rnd)
+{
+  mpfr_t ref, res, ab, cd;
+  int inex_ref, inex_res;
+  mpfr_prec_t p = MPFR_PREC(a);
+
+  mpfr_init2 (res, p);
+  mpfr_init2 (ref, p);
+  mpfr_init2 (ab, mpfr_get_prec (a) +  mpfr_get_prec (b));
+  mpfr_init2 (cd, mpfr_get_prec (c) +  mpfr_get_prec (d));
+
+  /* first check fmma */
+  inex_res = mpfr_fmma (res, a, b, c, d, rnd);
+  mpfr_mul (ab, a, b, rnd);
+  mpfr_mul (cd, c, d, rnd);
+  inex_ref = mpfr_add (ref, ab, cd, rnd);
+  if (inex_res != inex_ref)
+    {
+      mpfr_printf ("mpfr_fmma failed for p=%u a=%Rf b=%Rf c=%Rf d=%Rf rnd=%d\n",
+                   p, a, b, c, d, rnd);
+      printf ("expected inex %d, got %d\n", inex_ref, inex_res);
+      exit (1);
+    }
+  if (mpfr_nan_p (res))
+    MPFR_ASSERTN (mpfr_nan_p (ref));
+  else
+    MPFR_ASSERTN (mpfr_cmp (res, ref) == 0);
+
+  /* then check fmms */
+  inex_res = mpfr_fmms (res, a, b, c, d, rnd);
+  mpfr_mul (ab, a, b, rnd);
+  mpfr_mul (cd, c, d, rnd);
+  inex_ref = mpfr_sub (ref, ab, cd, rnd);
+  if (inex_res != inex_ref)
+    {
+      mpfr_printf ("mpfr_fmms failed for p=%u a=%Rf b=%Rf c=%Rf d=%Rf rnd=%d\n",
+                   p, a, b, c, d, rnd);
+      printf ("expected inex %d, got %d\n", inex_ref, inex_res);
+      exit (1);
+    }
+  if (mpfr_nan_p (res))
+    MPFR_ASSERTN (mpfr_nan_p (ref));
+  else
+    MPFR_ASSERTN (mpfr_cmp (res, ref) == 0);
+
+  mpfr_clear (ab);
+  mpfr_clear (cd);
+  mpfr_clear (res);
+  mpfr_clear (ref);
+}
+
+static void
+random_tests (void)
+{
+  mpfr_prec_t p;
+  mpfr_rnd_t r;
+  mpfr_t a, b, c, d;
+
+  for (p = MPFR_PREC_MIN; p <= 4096; p++)
+    {
+      mpfr_inits2 (p, a, b, c, d, (mpfr_ptr) 0);
+      mpfr_urandomb (a, RANDS);
+      mpfr_urandomb (b, RANDS);
+      mpfr_urandomb (c, RANDS);
+      mpfr_urandomb (d, RANDS);
+      for (r = 0; r < 4; r++)
+        random_test (a, b, c, d, r);
+      mpfr_clears (a, b, c, d, (mpfr_ptr) 0);
+    }
+}
+
+int
+main (int argc, char *argv[])
+{
+  tests_start_mpfr ();
+
+  random_tests ();
+
+  tests_end_mpfr ();
+  return 0;
+}