1 /* mpfr_fma -- Floating multiply-add 2 3 Copyright 2001-2002, 2004, 2006-2023 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramba projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 8 The GNU MPFR Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MPFR Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23 #define MPFR_NEED_LONGLONG_H 24 #include "mpfr-impl.h" 25 26 /* The fused-multiply-add (fma) of x, y and z is defined by: 27 fma(x,y,z)= x*y + z 28 */ 29 30 /* this function deals with all cases where inputs are singular, i.e., 31 either NaN, Inf or zero */ 32 static int 33 mpfr_fma_singular (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z, 34 mpfr_rnd_t rnd_mode) 35 { 36 if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y) || MPFR_IS_NAN(z)) 37 { 38 MPFR_SET_NAN(s); 39 MPFR_RET_NAN; 40 } 41 /* now neither x, y or z is NaN */ 42 else if (MPFR_IS_INF(x) || MPFR_IS_INF(y)) 43 { 44 /* cases Inf*0+z, 0*Inf+z, Inf-Inf */ 45 if ((MPFR_IS_ZERO(y)) || 46 (MPFR_IS_ZERO(x)) || 47 (MPFR_IS_INF(z) && 48 ((MPFR_MULT_SIGN(MPFR_SIGN(x), MPFR_SIGN(y))) != MPFR_SIGN(z)))) 49 { 50 MPFR_SET_NAN(s); 51 MPFR_RET_NAN; 52 } 53 else if (MPFR_IS_INF(z)) /* case Inf-Inf already checked above */ 54 { 55 MPFR_SET_INF(s); 56 MPFR_SET_SAME_SIGN(s, z); 57 MPFR_RET(0); 58 } 59 else /* z is finite */ 60 { 61 MPFR_SET_INF(s); 62 MPFR_SET_SIGN(s, MPFR_MULT_SIGN(MPFR_SIGN(x), MPFR_SIGN(y))); 63 MPFR_RET(0); 64 } 65 } 66 /* now x and y are finite */ 67 else if (MPFR_IS_INF(z)) 68 { 69 MPFR_SET_INF(s); 70 MPFR_SET_SAME_SIGN(s, z); 71 MPFR_RET(0); 72 } 73 else if (MPFR_IS_ZERO(x) || MPFR_IS_ZERO(y)) 74 { 75 if (MPFR_IS_ZERO(z)) 76 { 77 int sign_p; 78 sign_p = MPFR_MULT_SIGN(MPFR_SIGN(x), MPFR_SIGN(y)); 79 MPFR_SET_SIGN(s, (rnd_mode != MPFR_RNDD ? 80 (MPFR_IS_NEG_SIGN(sign_p) && MPFR_IS_NEG(z) ? 81 MPFR_SIGN_NEG : MPFR_SIGN_POS) : 82 (MPFR_IS_POS_SIGN(sign_p) && MPFR_IS_POS(z) ? 83 MPFR_SIGN_POS : MPFR_SIGN_NEG))); 84 MPFR_SET_ZERO(s); 85 MPFR_RET(0); 86 } 87 else 88 return mpfr_set (s, z, rnd_mode); 89 } 90 else /* necessarily z is zero here */ 91 { 92 MPFR_ASSERTD(MPFR_IS_ZERO(z)); 93 return (x == y) ? mpfr_sqr (s, x, rnd_mode) 94 : mpfr_mul (s, x, y, rnd_mode); 95 } 96 } 97 98 /* s <- x*y + z */ 99 int 100 mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z, 101 mpfr_rnd_t rnd_mode) 102 { 103 int inexact; 104 mpfr_t u; 105 mp_size_t n; 106 mpfr_exp_t e; 107 mpfr_prec_t precx, precy; 108 MPFR_SAVE_EXPO_DECL (expo); 109 MPFR_GROUP_DECL(group); 110 111 MPFR_LOG_FUNC 112 (("x[%Pd]=%.*Rg y[%Pd]=%.*Rg z[%Pd]=%.*Rg rnd=%d", 113 mpfr_get_prec (x), mpfr_log_prec, x, 114 mpfr_get_prec (y), mpfr_log_prec, y, 115 mpfr_get_prec (z), mpfr_log_prec, z, rnd_mode), 116 ("s[%Pd]=%.*Rg inexact=%d", 117 mpfr_get_prec (s), mpfr_log_prec, s, inexact)); 118 119 /* particular cases */ 120 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) || MPFR_IS_SINGULAR(y) || 121 MPFR_IS_SINGULAR(z) )) 122 return mpfr_fma_singular (s, x, y, z, rnd_mode); 123 124 e = MPFR_GET_EXP (x) + MPFR_GET_EXP (y); 125 126 precx = MPFR_PREC (x); 127 precy = MPFR_PREC (y); 128 129 /* First deal with special case where prec(x) = prec(y) and x*y does 130 not overflow nor underflow. Do it only for small sizes since for large 131 sizes x*y is faster using Mulders' algorithm (as a rule of thumb, 132 we assume mpn_mul_n is faster up to 4*MPFR_MUL_THRESHOLD). 133 Since |EXP(x)|, |EXP(y)| < 2^(k-2) on a k-bit computer, 134 |EXP(x)+EXP(y)| < 2^(k-1), thus cannot overflow nor underflow. */ 135 if (precx == precy && e <= __gmpfr_emax && e > __gmpfr_emin) 136 { 137 if (precx < GMP_NUMB_BITS && 138 MPFR_PREC(z) == precx && 139 MPFR_PREC(s) == precx) 140 { 141 mp_limb_t umant[2], zmant[2]; 142 mpfr_t zz; 143 int inex; 144 145 umul_ppmm (umant[1], umant[0], MPFR_MANT(x)[0], MPFR_MANT(y)[0]); 146 MPFR_PREC(u) = MPFR_PREC(zz) = 2 * precx; 147 MPFR_MANT(u) = umant; 148 MPFR_MANT(zz) = zmant; 149 MPFR_SIGN(u) = MPFR_MULT_SIGN(MPFR_SIGN(x), MPFR_SIGN(y)); 150 MPFR_SIGN(zz) = MPFR_SIGN(z); 151 MPFR_EXP(zz) = MPFR_EXP(z); 152 if (MPFR_PREC(zz) <= GMP_NUMB_BITS) /* zz fits in one limb */ 153 { 154 if ((umant[1] & MPFR_LIMB_HIGHBIT) == 0) 155 { 156 umant[0] = umant[1] << 1; 157 MPFR_EXP(u) = e - 1; 158 } 159 else 160 { 161 umant[0] = umant[1]; 162 MPFR_EXP(u) = e; 163 } 164 zmant[0] = MPFR_MANT(z)[0]; 165 } 166 else 167 { 168 zmant[1] = MPFR_MANT(z)[0]; 169 zmant[0] = MPFR_LIMB_ZERO; 170 if ((umant[1] & MPFR_LIMB_HIGHBIT) == 0) 171 { 172 umant[1] = (umant[1] << 1) | 173 (umant[0] >> (GMP_NUMB_BITS - 1)); 174 umant[0] = umant[0] << 1; 175 MPFR_EXP(u) = e - 1; 176 } 177 else 178 MPFR_EXP(u) = e; 179 } 180 inex = mpfr_add (u, u, zz, rnd_mode); 181 /* mpfr_set_1_2 requires PREC(u) = 2*PREC(s), 182 thus we need PREC(s) = PREC(x) = PREC(y) = PREC(z) */ 183 return mpfr_set_1_2 (s, u, rnd_mode, inex); 184 } 185 else if ((n = MPFR_LIMB_SIZE(x)) <= 4 * MPFR_MUL_THRESHOLD) 186 { 187 mpfr_limb_ptr up; 188 mp_size_t un = n + n; 189 MPFR_TMP_DECL(marker); 190 191 MPFR_TMP_MARK(marker); 192 MPFR_TMP_INIT (up, u, un * GMP_NUMB_BITS, un); 193 up = MPFR_MANT(u); 194 /* multiply x*y exactly into u */ 195 if (x == y) 196 mpn_sqr (up, MPFR_MANT(x), n); 197 else 198 mpn_mul_n (up, MPFR_MANT(x), MPFR_MANT(y), n); 199 if (MPFR_LIMB_MSB (up[un - 1]) == 0) 200 { 201 mpn_lshift (up, up, un, 1); 202 MPFR_EXP(u) = e - 1; 203 } 204 else 205 MPFR_EXP(u) = e; 206 MPFR_SIGN(u) = MPFR_MULT_SIGN(MPFR_SIGN(x), MPFR_SIGN(y)); 207 /* The above code does not generate any exception. 208 The exceptions will come only from mpfr_add. */ 209 inexact = mpfr_add (s, u, z, rnd_mode); 210 MPFR_TMP_FREE(marker); 211 return inexact; 212 } 213 } 214 215 /* If we take prec(u) >= prec(x) + prec(y), the product u <- x*y 216 is exact, except in case of overflow or underflow. */ 217 MPFR_ASSERTN (precx + precy <= MPFR_PREC_MAX); 218 MPFR_GROUP_INIT_1 (group, precx + precy, u); 219 MPFR_SAVE_EXPO_MARK (expo); 220 221 if (MPFR_UNLIKELY (mpfr_mul (u, x, y, MPFR_RNDN))) 222 { 223 /* overflow or underflow - this case is regarded as rare, thus 224 does not need to be very efficient (even if some tests below 225 could have been done earlier). 226 It is an overflow iff u is an infinity (since MPFR_RNDN was used). 227 Alternatively, we could test the overflow flag, but in this case, 228 mpfr_clear_flags would have been necessary. */ 229 230 if (MPFR_IS_INF (u)) /* overflow */ 231 { 232 int sign_u = MPFR_SIGN (u); 233 234 MPFR_LOG_MSG (("Overflow on x*y\n", 0)); 235 MPFR_GROUP_CLEAR (group); /* we no longer need u */ 236 237 /* Let's eliminate the obvious case where x*y and z have the 238 same sign. No possible cancellation -> real overflow. 239 Also, we know that |z| < 2^emax. If E(x) + E(y) >= emax+3, 240 then |x*y| >= 2^(emax+1), and |x*y + z| > 2^emax. This case 241 is also an overflow. */ 242 if (sign_u == MPFR_SIGN (z) || e >= __gmpfr_emax + 3) 243 { 244 MPFR_SAVE_EXPO_FREE (expo); 245 return mpfr_overflow (s, rnd_mode, sign_u); 246 } 247 } 248 else /* underflow: one has |x*y| < 2^(emin-1). */ 249 { 250 MPFR_LOG_MSG (("Underflow on x*y\n", 0)); 251 252 /* Easy cases: when 2^(emin-1) <= 1/2 * min(ulp(z),ulp(s)), 253 one can replace x*y by sign(x*y) * 2^(emin-1). Note that 254 this is even true in case of equality for MPFR_RNDN thanks 255 to the even-rounding rule. 256 The + 1 on MPFR_PREC (s) is necessary because the exponent 257 of the result can be EXP(z) - 1. */ 258 if (MPFR_GET_EXP (z) - __gmpfr_emin >= 259 MAX (MPFR_PREC (z), MPFR_PREC (s) + 1)) 260 { 261 MPFR_PREC (u) = MPFR_PREC_MIN; 262 mpfr_setmin (u, __gmpfr_emin); 263 MPFR_SET_SIGN (u, MPFR_MULT_SIGN (MPFR_SIGN (x), 264 MPFR_SIGN (y))); 265 mpfr_clear_flags (); 266 goto add; 267 } 268 269 MPFR_GROUP_CLEAR (group); /* we no longer need u */ 270 } 271 272 /* Let's use UBF to resolve the overflow/underflow issues. */ 273 { 274 mpfr_ubf_t uu; 275 mp_size_t un; 276 mpfr_limb_ptr up; 277 MPFR_TMP_DECL(marker); 278 279 MPFR_LOG_MSG (("Use UBF\n", 0)); 280 281 MPFR_TMP_MARK (marker); 282 un = MPFR_LIMB_SIZE (x) + MPFR_LIMB_SIZE (y); 283 MPFR_TMP_INIT (up, uu, (mpfr_prec_t) un * GMP_NUMB_BITS, un); 284 mpfr_ubf_mul_exact (uu, x, y); 285 mpfr_clear_flags (); 286 inexact = mpfr_add (s, (mpfr_srcptr) uu, z, rnd_mode); 287 MPFR_UBF_CLEAR_EXP (uu); 288 MPFR_TMP_FREE (marker); 289 } 290 } 291 else 292 { 293 add: 294 inexact = mpfr_add (s, u, z, rnd_mode); 295 MPFR_GROUP_CLEAR (group); 296 } 297 298 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); 299 MPFR_SAVE_EXPO_FREE (expo); 300 return mpfr_check_range (s, inexact, rnd_mode); 301 } 302
Indexes created Wed Mar 04 05:31:52 UTC 2026