diff options
author | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2016-01-12 12:29:04 +0000 |
---|---|---|
committer | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2016-01-12 12:29:04 +0000 |
commit | 737b3ca0582dab0db17534a558ae2b4f0eb66b6d (patch) | |
tree | 3ccff90e6d9406aeb22298e87a5021e6b27dff3e | |
parent | 30bec485e6fad5271d19c7999057457fe21edb02 (diff) | |
download | mpfr-737b3ca0582dab0db17534a558ae2b4f0eb66b6d.tar.gz |
forgot two files
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@9793 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r-- | src/fmma.c | 135 | ||||
-rw-r--r-- | tests/tfmma.c | 107 |
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; +} |