Home | History | Annotate | Line # | Download | only in src 
      1 /* mpc_fma -- Fused multiply-add of three complex numbers
      2 
      3 Copyright (C) 2011, 2012, 2022 INRIA
      4 
      5 This file is part of GNU MPC.
      6 
      7 GNU MPC is free software; you can redistribute it and/or modify it under
      8 the terms of the GNU Lesser General Public License as published by the
      9 Free Software Foundation; either version 3 of the License, or (at your
     10 option) any later version.
     11 
     12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
     13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
     14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
     15 more details.
     16 
     17 You should have received a copy of the GNU Lesser General Public License
     18 along with this program. If not, see http://www.gnu.org/licenses/ .
     19 */
     20 
     21 #include "mpc-impl.h"
     22 
     23 int
     24 mpc_fma_naive (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
     25 {
     26   mpfr_t rea_reb, rea_imb, ima_reb, ima_imb;
     27   mpfr_ptr sum [3];
     28   int inex_re, inex_im;
     29 
     30   mpfr_init2 (rea_reb, mpfr_get_prec (mpc_realref(a)) + mpfr_get_prec (mpc_realref(b)));
     31   mpfr_init2 (rea_imb, mpfr_get_prec (mpc_realref(a)) + mpfr_get_prec (mpc_imagref(b)));
     32   mpfr_init2 (ima_reb, mpfr_get_prec (mpc_imagref(a)) + mpfr_get_prec (mpc_realref(b)));
     33   mpfr_init2 (ima_imb, mpfr_get_prec (mpc_imagref(a)) + mpfr_get_prec (mpc_imagref(b)));
     34 
     35   mpfr_mul (rea_reb, mpc_realref(a), mpc_realref(b), MPFR_RNDZ); /* exact */
     36   mpfr_mul (rea_imb, mpc_realref(a), mpc_imagref(b), MPFR_RNDZ); /* exact */
     37   mpfr_mul (ima_reb, mpc_imagref(a), mpc_realref(b), MPFR_RNDZ); /* exact */
     38   mpfr_mul (ima_imb, mpc_imagref(a), mpc_imagref(b), MPFR_RNDZ); /* exact */
     39 
     40   mpfr_neg (ima_imb, ima_imb, MPFR_RNDZ);
     41   sum [0] = rea_reb;
     42   sum [1] = ima_imb;
     43   sum [2] = (mpfr_ptr) mpc_realref (c);
     44   inex_re = mpfr_sum (mpc_realref (r), sum, 3, MPC_RND_RE (rnd));
     45   sum [0] = rea_imb;
     46   sum [1] = ima_reb;
     47   sum [2] = (mpfr_ptr) mpc_imagref (c);
     48   inex_im = mpfr_sum (mpc_imagref (r), sum, 3, MPC_RND_IM (rnd));
     49 
     50   mpfr_clear (rea_reb);
     51   mpfr_clear (rea_imb);
     52   mpfr_clear (ima_reb);
     53   mpfr_clear (ima_imb);
     54 
     55   return MPC_INEX(inex_re, inex_im);
     56 }
     57 
     58 /* The algorithm is as follows:
     59    - in a first pass, we use the target precision + some extra bits
     60    - if it fails, we add the number of cancelled bits when adding
     61      Re(a*b) and Re(c) [similarly for the imaginary part]
     62    - it is fails again, we call the mpc_fma_naive function, which also
     63      deals with the special cases */
     64 int
     65 mpc_fma (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
     66 {
     67   mpc_t ab;
     68   mpfr_prec_t pre, pim, wpre, wpim;
     69   mpfr_exp_t diffre, diffim;
     70   int i, inex = 0, okre = 0, okim = 0;
     71 
     72   if (mpc_fin_p (a) == 0 || mpc_fin_p (b) == 0 || mpc_fin_p (c) == 0)
     73     return mpc_fma_naive (r, a, b, c, rnd);
     74 
     75   pre = mpfr_get_prec (mpc_realref(r));
     76   pim = mpfr_get_prec (mpc_imagref(r));
     77   wpre = pre + mpc_ceil_log2 (pre) + 10;
     78   wpim = pim + mpc_ceil_log2 (pim) + 10;
     79   mpc_init3 (ab, wpre, wpim);
     80   for (i = 0; i < 2; ++i)
     81     {
     82       mpc_mul (ab, a, b, MPC_RNDZZ);
     83       if (mpfr_zero_p (mpc_realref(ab)) || mpfr_zero_p (mpc_imagref(ab)))
     84         break;
     85       diffre = mpfr_get_exp (mpc_realref(ab));
     86       diffim = mpfr_get_exp (mpc_imagref(ab));
     87       mpc_add (ab, ab, c, MPC_RNDZZ);
     88       if (mpfr_zero_p (mpc_realref(ab)) || mpfr_zero_p (mpc_imagref(ab)))
     89         break;
     90       diffre -= mpfr_get_exp (mpc_realref(ab));
     91       diffim -= mpfr_get_exp (mpc_imagref(ab));
     92       diffre = (diffre > 0 ? diffre + 1 : 1);
     93       diffim = (diffim > 0 ? diffim + 1 : 1);
     94       okre = diffre > (mpfr_exp_t) wpre ? 0 : mpfr_can_round (mpc_realref(ab),
     95                                  wpre - diffre, MPFR_RNDN, MPFR_RNDZ,
     96                                  pre + (MPC_RND_RE (rnd) == MPFR_RNDN));
     97       okim = diffim > (mpfr_exp_t) wpim ? 0 : mpfr_can_round (mpc_imagref(ab),
     98                                  wpim - diffim, MPFR_RNDN, MPFR_RNDZ,
     99                                  pim + (MPC_RND_IM (rnd) == MPFR_RNDN));
    100       if (okre && okim)
    101         {
    102           inex = mpc_set (r, ab, rnd);
    103           break;
    104         }
    105       if (i == 1)
    106         break;
    107       if (okre == 0 && diffre > 1)
    108         wpre += diffre;
    109       if (okim == 0 && diffim > 1)
    110         wpim += diffim;
    111       mpfr_set_prec (mpc_realref(ab), wpre);
    112       mpfr_set_prec (mpc_imagref(ab), wpim);
    113     }
    114   mpc_clear (ab);
    115 
    116   return (okre && okim) ? inex : mpc_fma_naive (r, a, b, c, rnd);
    117 }
    118