|
Packit Service |
2e9770 |
/* mpc_fma -- Fused multiply-add of three complex numbers
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Copyright (C) 2011, 2012 INRIA
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
This file is part of GNU MPC.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
GNU MPC is free software; you can redistribute it and/or modify it under
|
|
Packit Service |
2e9770 |
the terms of the GNU Lesser General Public License as published by the
|
|
Packit Service |
2e9770 |
Free Software Foundation; either version 3 of the License, or (at your
|
|
Packit Service |
2e9770 |
option) any later version.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
Packit Service |
2e9770 |
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
Packit Service |
2e9770 |
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
|
|
Packit Service |
2e9770 |
more details.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
You should have received a copy of the GNU Lesser General Public License
|
|
Packit Service |
2e9770 |
along with this program. If not, see http://www.gnu.org/licenses/ .
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
#include "mpc-impl.h"
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* return a bound on the precision needed to add or subtract x and y exactly */
|
|
Packit Service |
2e9770 |
static mpfr_prec_t
|
|
Packit Service |
2e9770 |
bound_prec_addsub (mpfr_srcptr x, mpfr_srcptr y)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (!mpfr_regular_p (x))
|
|
Packit Service |
2e9770 |
return mpfr_get_prec (y);
|
|
Packit Service |
2e9770 |
else if (!mpfr_regular_p (y))
|
|
Packit Service |
2e9770 |
return mpfr_get_prec (x);
|
|
Packit Service |
2e9770 |
else /* neither x nor y are NaN, Inf or zero */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_exp_t ex = mpfr_get_exp (x);
|
|
Packit Service |
2e9770 |
mpfr_exp_t ey = mpfr_get_exp (y);
|
|
Packit Service |
2e9770 |
mpfr_exp_t ulpx = ex - mpfr_get_prec (x);
|
|
Packit Service |
2e9770 |
mpfr_exp_t ulpy = ey - mpfr_get_prec (y);
|
|
Packit Service |
2e9770 |
return ((ex >= ey) ? ex : ey) + 1 - ((ulpx <= ulpy) ? ulpx : ulpy);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* r <- a*b+c */
|
|
Packit Service |
2e9770 |
int
|
|
Packit Service |
2e9770 |
mpc_fma_naive (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_t rea_reb, rea_imb, ima_reb, ima_imb, tmp;
|
|
Packit Service |
2e9770 |
mpfr_prec_t pre12, pre13, pre23, pim12, pim13, pim23;
|
|
Packit Service |
2e9770 |
int inex_re, inex_im;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_init2 (rea_reb, mpfr_get_prec (mpc_realref(a)) + mpfr_get_prec (mpc_realref(b)));
|
|
Packit Service |
2e9770 |
mpfr_init2 (rea_imb, mpfr_get_prec (mpc_realref(a)) + mpfr_get_prec (mpc_imagref(b)));
|
|
Packit Service |
2e9770 |
mpfr_init2 (ima_reb, mpfr_get_prec (mpc_imagref(a)) + mpfr_get_prec (mpc_realref(b)));
|
|
Packit Service |
2e9770 |
mpfr_init2 (ima_imb, mpfr_get_prec (mpc_imagref(a)) + mpfr_get_prec (mpc_imagref(b)));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_mul (rea_reb, mpc_realref(a), mpc_realref(b), MPFR_RNDZ); /* exact */
|
|
Packit Service |
2e9770 |
mpfr_mul (rea_imb, mpc_realref(a), mpc_imagref(b), MPFR_RNDZ); /* exact */
|
|
Packit Service |
2e9770 |
mpfr_mul (ima_reb, mpc_imagref(a), mpc_realref(b), MPFR_RNDZ); /* exact */
|
|
Packit Service |
2e9770 |
mpfr_mul (ima_imb, mpc_imagref(a), mpc_imagref(b), MPFR_RNDZ); /* exact */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* Re(r) <- rea_reb - ima_imb + Re(c) */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
pre12 = bound_prec_addsub (rea_reb, ima_imb); /* bound on exact precision for
|
|
Packit Service |
2e9770 |
rea_reb - ima_imb */
|
|
Packit Service |
2e9770 |
pre13 = bound_prec_addsub (rea_reb, mpc_realref(c));
|
|
Packit Service |
2e9770 |
/* bound for rea_reb + Re(c) */
|
|
Packit Service |
2e9770 |
pre23 = bound_prec_addsub (ima_imb, mpc_realref(c));
|
|
Packit Service |
2e9770 |
/* bound for ima_imb - Re(c) */
|
|
Packit Service |
2e9770 |
if (pre12 <= pre13 && pre12 <= pre23) /* (rea_reb - ima_imb) + Re(c) */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_init2 (tmp, pre12);
|
|
Packit Service |
2e9770 |
mpfr_sub (tmp, rea_reb, ima_imb, MPFR_RNDZ); /* exact */
|
|
Packit Service |
2e9770 |
inex_re = mpfr_add (mpc_realref(r), tmp, mpc_realref(c), MPC_RND_RE(rnd));
|
|
Packit Service |
2e9770 |
/* the only possible bad overlap is between r and c, but since we are
|
|
Packit Service |
2e9770 |
only touching the real part of both, it is ok */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (pre13 <= pre23) /* (rea_reb + Re(c)) - ima_imb */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_init2 (tmp, pre13);
|
|
Packit Service |
2e9770 |
mpfr_add (tmp, rea_reb, mpc_realref(c), MPFR_RNDZ); /* exact */
|
|
Packit Service |
2e9770 |
inex_re = mpfr_sub (mpc_realref(r), tmp, ima_imb, MPC_RND_RE(rnd));
|
|
Packit Service |
2e9770 |
/* the only possible bad overlap is between r and c, but since we are
|
|
Packit Service |
2e9770 |
only touching the real part of both, it is ok */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else /* rea_reb + (Re(c) - ima_imb) */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_init2 (tmp, pre23);
|
|
Packit Service |
2e9770 |
mpfr_sub (tmp, mpc_realref(c), ima_imb, MPFR_RNDZ); /* exact */
|
|
Packit Service |
2e9770 |
inex_re = mpfr_add (mpc_realref(r), tmp, rea_reb, MPC_RND_RE(rnd));
|
|
Packit Service |
2e9770 |
/* the only possible bad overlap is between r and c, but since we are
|
|
Packit Service |
2e9770 |
only touching the real part of both, it is ok */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* Im(r) <- rea_imb + ima_reb + Im(c) */
|
|
Packit Service |
2e9770 |
pim12 = bound_prec_addsub (rea_imb, ima_reb); /* bound on exact precision for
|
|
Packit Service |
2e9770 |
rea_imb + ima_reb */
|
|
Packit Service |
2e9770 |
pim13 = bound_prec_addsub (rea_imb, mpc_imagref(c));
|
|
Packit Service |
2e9770 |
/* bound for rea_imb + Im(c) */
|
|
Packit Service |
2e9770 |
pim23 = bound_prec_addsub (ima_reb, mpc_imagref(c));
|
|
Packit Service |
2e9770 |
/* bound for ima_reb + Im(c) */
|
|
Packit Service |
2e9770 |
if (pim12 <= pim13 && pim12 <= pim23) /* (rea_imb + ima_reb) + Im(c) */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_prec (tmp, pim12);
|
|
Packit Service |
2e9770 |
mpfr_add (tmp, rea_imb, ima_reb, MPFR_RNDZ); /* exact */
|
|
Packit Service |
2e9770 |
inex_im = mpfr_add (mpc_imagref(r), tmp, mpc_imagref(c), MPC_RND_IM(rnd));
|
|
Packit Service |
2e9770 |
/* the only possible bad overlap is between r and c, but since we are
|
|
Packit Service |
2e9770 |
only touching the imaginary part of both, it is ok */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (pim13 <= pim23) /* (rea_imb + Im(c)) + ima_reb */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_prec (tmp, pim13);
|
|
Packit Service |
2e9770 |
mpfr_add (tmp, rea_imb, mpc_imagref(c), MPFR_RNDZ); /* exact */
|
|
Packit Service |
2e9770 |
inex_im = mpfr_add (mpc_imagref(r), tmp, ima_reb, MPC_RND_IM(rnd));
|
|
Packit Service |
2e9770 |
/* the only possible bad overlap is between r and c, but since we are
|
|
Packit Service |
2e9770 |
only touching the imaginary part of both, it is ok */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else /* rea_imb + (Im(c) + ima_reb) */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_prec (tmp, pre23);
|
|
Packit Service |
2e9770 |
mpfr_add (tmp, mpc_imagref(c), ima_reb, MPFR_RNDZ); /* exact */
|
|
Packit Service |
2e9770 |
inex_im = mpfr_add (mpc_imagref(r), tmp, rea_imb, MPC_RND_IM(rnd));
|
|
Packit Service |
2e9770 |
/* the only possible bad overlap is between r and c, but since we are
|
|
Packit Service |
2e9770 |
only touching the imaginary part of both, it is ok */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_clear (rea_reb);
|
|
Packit Service |
2e9770 |
mpfr_clear (rea_imb);
|
|
Packit Service |
2e9770 |
mpfr_clear (ima_reb);
|
|
Packit Service |
2e9770 |
mpfr_clear (ima_imb);
|
|
Packit Service |
2e9770 |
mpfr_clear (tmp);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return MPC_INEX(inex_re, inex_im);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* The algorithm is as follows:
|
|
Packit Service |
2e9770 |
- in a first pass, we use the target precision + some extra bits
|
|
Packit Service |
2e9770 |
- if it fails, we add the number of cancelled bits when adding
|
|
Packit Service |
2e9770 |
Re(a*b) and Re(c) [similarly for the imaginary part]
|
|
Packit Service |
2e9770 |
- it is fails again, we call the mpc_fma_naive function, which also
|
|
Packit Service |
2e9770 |
deals with the special cases */
|
|
Packit Service |
2e9770 |
int
|
|
Packit Service |
2e9770 |
mpc_fma (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpc_t ab;
|
|
Packit Service |
2e9770 |
mpfr_prec_t pre, pim, wpre, wpim;
|
|
Packit Service |
2e9770 |
mpfr_exp_t diffre, diffim;
|
|
Packit Service |
2e9770 |
int i, inex = 0, okre = 0, okim = 0;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (mpc_fin_p (a) == 0 || mpc_fin_p (b) == 0 || mpc_fin_p (c) == 0)
|
|
Packit Service |
2e9770 |
return mpc_fma_naive (r, a, b, c, rnd);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
pre = mpfr_get_prec (mpc_realref(r));
|
|
Packit Service |
2e9770 |
pim = mpfr_get_prec (mpc_imagref(r));
|
|
Packit Service |
2e9770 |
wpre = pre + mpc_ceil_log2 (pre) + 10;
|
|
Packit Service |
2e9770 |
wpim = pim + mpc_ceil_log2 (pim) + 10;
|
|
Packit Service |
2e9770 |
mpc_init3 (ab, wpre, wpim);
|
|
Packit Service |
2e9770 |
for (i = 0; i < 2; ++i)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpc_mul (ab, a, b, MPC_RNDZZ);
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_realref(ab)) || mpfr_zero_p (mpc_imagref(ab)))
|
|
Packit Service |
2e9770 |
break;
|
|
Packit Service |
2e9770 |
diffre = mpfr_get_exp (mpc_realref(ab));
|
|
Packit Service |
2e9770 |
diffim = mpfr_get_exp (mpc_imagref(ab));
|
|
Packit Service |
2e9770 |
mpc_add (ab, ab, c, MPC_RNDZZ);
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_realref(ab)) || mpfr_zero_p (mpc_imagref(ab)))
|
|
Packit Service |
2e9770 |
break;
|
|
Packit Service |
2e9770 |
diffre -= mpfr_get_exp (mpc_realref(ab));
|
|
Packit Service |
2e9770 |
diffim -= mpfr_get_exp (mpc_imagref(ab));
|
|
Packit Service |
2e9770 |
diffre = (diffre > 0 ? diffre + 1 : 1);
|
|
Packit Service |
2e9770 |
diffim = (diffim > 0 ? diffim + 1 : 1);
|
|
Packit Service |
2e9770 |
okre = diffre > (mpfr_exp_t) wpre ? 0 : mpfr_can_round (mpc_realref(ab),
|
|
Packit Service |
2e9770 |
wpre - diffre, MPFR_RNDN, MPFR_RNDZ,
|
|
Packit Service |
2e9770 |
pre + (MPC_RND_RE (rnd) == MPFR_RNDN));
|
|
Packit Service |
2e9770 |
okim = diffim > (mpfr_exp_t) wpim ? 0 : mpfr_can_round (mpc_imagref(ab),
|
|
Packit Service |
2e9770 |
wpim - diffim, MPFR_RNDN, MPFR_RNDZ,
|
|
Packit Service |
2e9770 |
pim + (MPC_RND_IM (rnd) == MPFR_RNDN));
|
|
Packit Service |
2e9770 |
if (okre && okim)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
inex = mpc_set (r, ab, rnd);
|
|
Packit Service |
2e9770 |
break;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
if (i == 1)
|
|
Packit Service |
2e9770 |
break;
|
|
Packit Service |
2e9770 |
if (okre == 0 && diffre > 1)
|
|
Packit Service |
2e9770 |
wpre += diffre;
|
|
Packit Service |
2e9770 |
if (okim == 0 && diffim > 1)
|
|
Packit Service |
2e9770 |
wpim += diffim;
|
|
Packit Service |
2e9770 |
mpfr_set_prec (mpc_realref(ab), wpre);
|
|
Packit Service |
2e9770 |
mpfr_set_prec (mpc_imagref(ab), wpim);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
mpc_clear (ab);
|
|
Packit Service |
2e9770 |
return okre && okim ? inex : mpc_fma_naive (r, a, b, c, rnd);
|
|
Packit Service |
2e9770 |
}
|