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