1 /* mpfr_fms -- Floating multiply-subtract
3 Copyright 2001, 2002, 2004, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 #include "mpfr-impl.h"
25 /* The fused-multiply-subtract (fms) of x, y and z is defined by:
27 Note: this is neither in IEEE754R, nor in LIA-2, but both the
28 PowerPC and the Itanium define fms as x*y - z.
32 mpfr_fms (mpfr_ptr s
, mpfr_srcptr x
, mpfr_srcptr y
, mpfr_srcptr z
,
37 MPFR_SAVE_EXPO_DECL (expo
);
38 MPFR_GROUP_DECL(group
);
41 (("x[%Pu]=%.*Rg y[%Pu]=%.*Rg z[%Pu]=%.*Rg rnd=%d",
42 mpfr_get_prec (x
), mpfr_log_prec
, x
,
43 mpfr_get_prec (y
), mpfr_log_prec
, y
,
44 mpfr_get_prec (z
), mpfr_log_prec
, z
, rnd_mode
),
45 ("s[%Pu]=%.*Rg inexact=%d",
46 mpfr_get_prec (s
), mpfr_log_prec
, s
, inexact
));
48 /* particular cases */
49 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x
) ||
50 MPFR_IS_SINGULAR(y
) ||
51 MPFR_IS_SINGULAR(z
) ))
53 if (MPFR_IS_NAN(x
) || MPFR_IS_NAN(y
) || MPFR_IS_NAN(z
))
58 /* now neither x, y or z is NaN */
59 else if (MPFR_IS_INF(x
) || MPFR_IS_INF(y
))
61 /* cases Inf*0-z, 0*Inf-z, Inf-Inf */
62 if ((MPFR_IS_ZERO(y
)) ||
65 ((MPFR_MULT_SIGN(MPFR_SIGN(x
), MPFR_SIGN(y
))) == MPFR_SIGN(z
))))
70 else if (MPFR_IS_INF(z
)) /* case Inf-Inf already checked above */
73 MPFR_SET_OPPOSITE_SIGN(s
, z
);
76 else /* z is finite */
79 MPFR_SET_SIGN(s
, MPFR_MULT_SIGN(MPFR_SIGN(x
) , MPFR_SIGN(y
)));
83 /* now x and y are finite */
84 else if (MPFR_IS_INF(z
))
87 MPFR_SET_OPPOSITE_SIGN(s
, z
);
90 else if (MPFR_IS_ZERO(x
) || MPFR_IS_ZERO(y
))
95 sign_p
= MPFR_MULT_SIGN( MPFR_SIGN(x
) , MPFR_SIGN(y
) );
96 MPFR_SET_SIGN(s
,(rnd_mode
!= MPFR_RNDD
?
97 ((MPFR_IS_NEG_SIGN(sign_p
) && MPFR_IS_POS(z
))
99 ((MPFR_IS_POS_SIGN(sign_p
) && MPFR_IS_NEG(z
))
105 return mpfr_neg (s
, z
, rnd_mode
);
107 else /* necessarily z is zero here */
109 MPFR_ASSERTD(MPFR_IS_ZERO(z
));
110 return mpfr_mul (s
, x
, y
, rnd_mode
);
114 /* If we take prec(u) >= prec(x) + prec(y), the product u <- x*y
115 is exact, except in case of overflow or underflow. */
116 MPFR_SAVE_EXPO_MARK (expo
);
117 MPFR_GROUP_INIT_1 (group
, MPFR_PREC(x
) + MPFR_PREC(y
), u
);
119 if (MPFR_UNLIKELY (mpfr_mul (u
, x
, y
, MPFR_RNDN
)))
121 /* overflow or underflow - this case is regarded as rare, thus
122 does not need to be very efficient (even if some tests below
123 could have been done earlier).
124 It is an overflow iff u is an infinity (since MPFR_RNDN was used).
125 Alternatively, we could test the overflow flag, but in this case,
126 mpfr_clear_flags would have been necessary. */
127 if (MPFR_IS_INF (u
)) /* overflow */
129 /* Let's eliminate the obvious case where x*y and z have the
130 same sign. No possible cancellation -> real overflow.
131 Also, we know that |z| < 2^emax. If E(x) + E(y) >= emax+3,
132 then |x*y| >= 2^(emax+1), and |x*y - z| >= 2^emax. This case
133 is also an overflow. */
134 if (MPFR_SIGN (u
) != MPFR_SIGN (z
) ||
135 MPFR_GET_EXP (x
) + MPFR_GET_EXP (y
) >= __gmpfr_emax
+ 3)
137 MPFR_GROUP_CLEAR (group
);
138 MPFR_SAVE_EXPO_FREE (expo
);
139 return mpfr_overflow (s
, rnd_mode
, - MPFR_SIGN (z
));
142 /* E(x) + E(y) <= emax+2, therefore |x*y| < 2^(emax+2), and
143 (x/4)*y does not overflow (let's recall that the result
144 is exact with an unbounded exponent range). It does not
145 underflow either, because x*y overflows and the exponent
146 range is large enough. */
147 inexact
= mpfr_div_2ui (u
, x
, 2, MPFR_RNDN
);
148 MPFR_ASSERTN (inexact
== 0);
149 inexact
= mpfr_mul (u
, u
, y
, MPFR_RNDN
);
150 MPFR_ASSERTN (inexact
== 0);
152 /* Now, we need to subtract z/4... But it may underflow! */
156 MPFR_BLOCK_DECL (flags
);
158 if (MPFR_GET_EXP (u
) > MPFR_GET_EXP (z
) &&
159 MPFR_GET_EXP (u
) - MPFR_GET_EXP (z
) > MPFR_PREC (u
))
161 /* |z| < ulp(u)/2, therefore one can use z instead of z/4. */
166 mpfr_init2 (zo4
, MPFR_PREC (z
));
167 if (mpfr_div_2ui (zo4
, z
, 2, MPFR_RNDZ
))
169 /* The division by 4 underflowed! */
170 MPFR_ASSERTN (0); /* TODO... */
175 /* Let's recall that u = x*y/4 and zz = z/4 (or z if the
176 following subtraction would give the same result). */
177 MPFR_BLOCK (flags
, inexact
= mpfr_sub (s
, u
, zz
, rnd_mode
));
178 /* u and zz have the same sign, so that an overflow
179 is not possible. But an underflow is theoretically
181 if (MPFR_UNDERFLOW (flags
))
183 MPFR_ASSERTN (zz
!= z
);
184 MPFR_ASSERTN (0); /* TODO... */
185 mpfr_clears (zo4
, u
, (mpfr_ptr
) 0);
193 MPFR_GROUP_CLEAR (group
);
194 MPFR_ASSERTN (! MPFR_OVERFLOW (flags
));
195 inex2
= mpfr_mul_2ui (s
, s
, 2, rnd_mode
);
196 if (inex2
) /* overflow */
199 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, __gmpfr_flags
);
205 else /* underflow: one has |xy| < 2^(emin-1). */
207 unsigned long scale
= 0;
214 /* Let's scale z so that ulp(z) > 2^emin and ulp(s) > 2^emin
215 (the + 1 on MPFR_PREC (s) is necessary because the exponent
216 of the result can be EXP(z) - 1). */
217 diffexp
= MPFR_GET_EXP (z
) - __gmpfr_emin
;
218 pzs
= MAX (MPFR_PREC (z
), MPFR_PREC (s
) + 1);
223 MPFR_BLOCK_DECL (flags
);
225 uscale
= (mpfr_uexp_t
) pzs
- diffexp
+ 1;
226 MPFR_ASSERTN (uscale
> 0);
227 MPFR_ASSERTN (uscale
<= ULONG_MAX
);
229 mpfr_init2 (scaled_z
, MPFR_PREC (z
));
230 inexact
= mpfr_mul_2ui (scaled_z
, z
, scale
, MPFR_RNDN
);
231 MPFR_ASSERTN (inexact
== 0); /* TODO: overflow case */
233 /* Now we need to recompute u = xy * 2^scale. */
235 if (MPFR_GET_EXP (x
) < MPFR_GET_EXP (y
))
237 mpfr_init2 (scaled_v
, MPFR_PREC (x
));
238 mpfr_mul_2ui (scaled_v
, x
, scale
, MPFR_RNDN
);
239 mpfr_mul (u
, scaled_v
, y
, MPFR_RNDN
);
243 mpfr_init2 (scaled_v
, MPFR_PREC (y
));
244 mpfr_mul_2ui (scaled_v
, y
, scale
, MPFR_RNDN
);
245 mpfr_mul (u
, x
, scaled_v
, MPFR_RNDN
);
247 mpfr_clear (scaled_v
);
248 MPFR_ASSERTN (! MPFR_OVERFLOW (flags
));
249 xy_underflows
= MPFR_UNDERFLOW (flags
);
259 /* Let's replace xy by sign(xy) * 2^(emin-1). */
260 MPFR_PREC (u
) = MPFR_PREC_MIN
;
261 mpfr_setmin (u
, __gmpfr_emin
);
262 MPFR_SET_SIGN (u
, MPFR_MULT_SIGN (MPFR_SIGN (x
),
267 MPFR_BLOCK_DECL (flags
);
269 MPFR_BLOCK (flags
, inexact
= mpfr_sub (s
, u
, new_z
, rnd_mode
));
270 MPFR_GROUP_CLEAR (group
);
275 mpfr_clear (scaled_z
);
276 /* Here an overflow is theoretically possible, in which case
277 the result may be wrong, hence the assert. An underflow
278 is not possible, but let's check that anyway. */
279 MPFR_ASSERTN (! MPFR_OVERFLOW (flags
)); /* TODO... */
280 MPFR_ASSERTN (! MPFR_UNDERFLOW (flags
)); /* not possible */
281 inex2
= mpfr_div_2ui (s
, s
, scale
, MPFR_RNDN
);
282 /* FIXME: this seems incorrect. MPFR_RNDN -> rnd_mode?
283 Also, handle the double rounding case:
284 s / 2^scale = 2^(emin - 2) in MPFR_RNDN. */
285 if (inex2
) /* underflow */
290 /* FIXME/TODO: I'm not sure that the following is correct.
291 Check for possible spurious exceptions due to intermediate
293 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, __gmpfr_flags
);
298 inexact
= mpfr_sub (s
, u
, z
, rnd_mode
);
299 MPFR_GROUP_CLEAR (group
);
300 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, __gmpfr_flags
);
302 MPFR_SAVE_EXPO_FREE (expo
);
303 return mpfr_check_range (s
, inexact
, rnd_mode
);