1 /* mpfr_fms -- Floating multiply-subtract
3 Copyright 2001, 2002, 2004, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao projects, INRIA.
6 This file is part of the GNU MPFR Library.
8 The GNU MPFR Library is free software; you can redistribute
9 it and/or modify it under the terms of the GNU Lesser
10 General Public License as published by the Free Software
11 Foundation; either version 2.1 of the License, or (at your
12 option) any later version.
14 The GNU MPFR Library is distributed in the hope that it will
15 be useful, but WITHOUT ANY WARRANTY; without even the
16 implied warranty of MERCHANTABILITY or FITNESS FOR A
17 PARTICULAR PURPOSE. See the GNU Lesser General Public
18 License for more details.
20 You should have received a copy of the GNU Lesser
21 General Public License along with the GNU MPFR Library; see
22 the file COPYING.LIB. If not, write to the Free Software
23 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
24 MA 02110-1301, USA. */
26 #include "mpfr-impl.h"
28 /* The fused-multiply-subtract (fms) of x, y and z is defined by:
30 Note: this is neither in IEEE754R, nor in LIA-2, but both the
31 PowerPC and the Itanium define fms as x*y - z.
35 mpfr_fms (mpfr_ptr s
, mpfr_srcptr x
, mpfr_srcptr y
, mpfr_srcptr z
,
40 MPFR_SAVE_EXPO_DECL (expo
);
42 /* particular cases */
43 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x
) ||
44 MPFR_IS_SINGULAR(y
) ||
45 MPFR_IS_SINGULAR(z
) ))
47 if (MPFR_IS_NAN(x
) || MPFR_IS_NAN(y
) || MPFR_IS_NAN(z
))
52 /* now neither x, y or z is NaN */
53 else if (MPFR_IS_INF(x
) || MPFR_IS_INF(y
))
55 /* cases Inf*0-z, 0*Inf-z, Inf-Inf */
56 if ((MPFR_IS_ZERO(y
)) ||
59 ((MPFR_MULT_SIGN(MPFR_SIGN(x
), MPFR_SIGN(y
))) == MPFR_SIGN(z
))))
64 else if (MPFR_IS_INF(z
)) /* case Inf-Inf already checked above */
67 MPFR_SET_OPPOSITE_SIGN(s
, z
);
70 else /* z is finite */
73 MPFR_SET_SIGN(s
, MPFR_MULT_SIGN(MPFR_SIGN(x
) , MPFR_SIGN(y
)));
77 /* now x and y are finite */
78 else if (MPFR_IS_INF(z
))
81 MPFR_SET_OPPOSITE_SIGN(s
, z
);
84 else if (MPFR_IS_ZERO(x
) || MPFR_IS_ZERO(y
))
89 sign_p
= MPFR_MULT_SIGN( MPFR_SIGN(x
) , MPFR_SIGN(y
) );
90 MPFR_SET_SIGN(s
,(rnd_mode
!= GMP_RNDD
?
91 ((MPFR_IS_NEG_SIGN(sign_p
) && MPFR_IS_POS(z
))
93 ((MPFR_IS_POS_SIGN(sign_p
) && MPFR_IS_NEG(z
))
99 return mpfr_neg (s
, z
, rnd_mode
);
101 else /* necessarily z is zero here */
103 MPFR_ASSERTD(MPFR_IS_ZERO(z
));
104 return mpfr_mul (s
, x
, y
, rnd_mode
);
108 /* If we take prec(u) >= prec(x) + prec(y), the product u <- x*y
109 is exact, except in case of overflow or underflow. */
110 MPFR_SAVE_EXPO_MARK (expo
);
111 mpfr_init2 (u
, MPFR_PREC(x
) + MPFR_PREC(y
));
113 if (MPFR_UNLIKELY (mpfr_mul (u
, x
, y
, GMP_RNDN
)))
115 /* overflow or underflow - this case is regarded as rare, thus
116 does not need to be very efficient (even if some tests below
117 could have been done earlier).
118 It is an overflow iff u is an infinity (since GMP_RNDN was used).
119 Alternatively, we could test the overflow flag, but in this case,
120 mpfr_clear_flags would have been necessary. */
121 if (MPFR_IS_INF (u
)) /* overflow */
123 /* Let's eliminate the obvious case where x*y and z have the
124 same sign. No possible cancellation -> real overflow.
125 Also, we know that |z| < 2^emax. If E(x) + E(y) >= emax+3,
126 then |x*y| >= 2^(emax+1), and |x*y - z| >= 2^emax. This case
127 is also an overflow. */
128 if (MPFR_SIGN (u
) != MPFR_SIGN (z
) ||
129 MPFR_GET_EXP (x
) + MPFR_GET_EXP (y
) >= __gmpfr_emax
+ 3)
132 MPFR_SAVE_EXPO_FREE (expo
);
133 return mpfr_overflow (s
, rnd_mode
, - MPFR_SIGN (z
));
136 /* E(x) + E(y) <= emax+2, therefore |x*y| < 2^(emax+2), and
137 (x/4)*y does not overflow (let's recall that the result
138 is exact with an unbounded exponent range). It does not
139 underflow either, because x*y overflows and the exponent
140 range is large enough. */
141 inexact
= mpfr_div_2ui (u
, x
, 2, GMP_RNDN
);
142 MPFR_ASSERTN (inexact
== 0);
143 inexact
= mpfr_mul (u
, u
, y
, GMP_RNDN
);
144 MPFR_ASSERTN (inexact
== 0);
146 /* Now, we need to subtract z/4... But it may underflow! */
150 MPFR_BLOCK_DECL (flags
);
152 if (MPFR_GET_EXP (u
) > MPFR_GET_EXP (z
) &&
153 MPFR_GET_EXP (u
) - MPFR_GET_EXP (z
) > MPFR_PREC (u
))
155 /* |z| < ulp(u)/2, therefore one can use z instead of z/4. */
160 mpfr_init2 (zo4
, MPFR_PREC (z
));
161 if (mpfr_div_2ui (zo4
, z
, 2, GMP_RNDZ
))
163 /* The division by 4 underflowed! */
164 MPFR_ASSERTN (0); /* TODO... */
169 /* Let's recall that u = x*y/4 and zz = z/4 (or z if the
170 following subtraction would give the same result). */
171 MPFR_BLOCK (flags
, inexact
= mpfr_sub (s
, u
, zz
, rnd_mode
));
172 /* u and zz have the same sign, so that an overflow
173 is not possible. But an underflow is theoretically
175 if (MPFR_UNDERFLOW (flags
))
177 MPFR_ASSERTN (zz
!= z
);
178 MPFR_ASSERTN (0); /* TODO... */
179 mpfr_clears (zo4
, u
, (mpfr_ptr
) 0);
188 MPFR_ASSERTN (! MPFR_OVERFLOW (flags
));
189 inex2
= mpfr_mul_2ui (s
, s
, 2, rnd_mode
);
190 if (inex2
) /* overflow */
193 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, __gmpfr_flags
);
199 else /* underflow: one has |xy| < 2^(emin-1). */
201 unsigned long scale
= 0;
208 /* Let's scale z so that ulp(z) > 2^emin and ulp(s) > 2^emin
209 (the + 1 on MPFR_PREC (s) is necessary because the exponent
210 of the result can be EXP(z) - 1). */
211 diffexp
= MPFR_GET_EXP (z
) - __gmpfr_emin
;
212 pzs
= MAX (MPFR_PREC (z
), MPFR_PREC (s
) + 1);
215 mp_exp_unsigned_t uscale
;
217 MPFR_BLOCK_DECL (flags
);
219 uscale
= (mp_exp_unsigned_t
) pzs
- diffexp
+ 1;
220 MPFR_ASSERTN (uscale
> 0);
221 MPFR_ASSERTN (uscale
<= ULONG_MAX
);
223 mpfr_init2 (scaled_z
, MPFR_PREC (z
));
224 inexact
= mpfr_mul_2ui (scaled_z
, z
, scale
, GMP_RNDN
);
225 MPFR_ASSERTN (inexact
== 0); /* TODO: overflow case */
227 /* Now we need to recompute u = xy * 2^scale. */
229 if (MPFR_GET_EXP (x
) < MPFR_GET_EXP (y
))
231 mpfr_init2 (scaled_v
, MPFR_PREC (x
));
232 mpfr_mul_2ui (scaled_v
, x
, scale
, GMP_RNDN
);
233 mpfr_mul (u
, scaled_v
, y
, GMP_RNDN
);
237 mpfr_init2 (scaled_v
, MPFR_PREC (y
));
238 mpfr_mul_2ui (scaled_v
, y
, scale
, GMP_RNDN
);
239 mpfr_mul (u
, x
, scaled_v
, GMP_RNDN
);
241 mpfr_clear (scaled_v
);
242 MPFR_ASSERTN (! MPFR_OVERFLOW (flags
));
243 xy_underflows
= MPFR_UNDERFLOW (flags
);
253 /* Let's replace xy by sign(xy) * 2^(emin-1). */
254 mpfr_set_prec (u
, MPFR_PREC_MIN
);
255 mpfr_setmin (u
, __gmpfr_emin
);
256 MPFR_SET_SIGN (u
, MPFR_MULT_SIGN (MPFR_SIGN (x
),
261 MPFR_BLOCK_DECL (flags
);
263 MPFR_BLOCK (flags
, inexact
= mpfr_sub (s
, u
, new_z
, rnd_mode
));
270 mpfr_clear (scaled_z
);
271 /* Here an overflow is theoretically possible, in which case
272 the result may be wrong, hence the assert. An underflow
273 is not possible, but let's check that anyway. */
274 MPFR_ASSERTN (! MPFR_OVERFLOW (flags
)); /* TODO... */
275 MPFR_ASSERTN (! MPFR_UNDERFLOW (flags
)); /* not possible */
276 inex2
= mpfr_div_2ui (s
, s
, scale
, GMP_RNDN
);
277 /* FIXME: this seems incorrect. GMP_RNDN -> rnd_mode?
278 Also, handle the double rounding case:
279 s / 2^scale = 2^(emin - 2) in GMP_RNDN. */
280 if (inex2
) /* underflow */
285 /* FIXME/TODO: I'm not sure that the following is correct.
286 Check for possible spurious exceptions due to intermediate
288 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, __gmpfr_flags
);
293 inexact
= mpfr_sub (s
, u
, z
, rnd_mode
);
295 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, __gmpfr_flags
);
297 MPFR_SAVE_EXPO_FREE (expo
);
298 return mpfr_check_range (s
, inexact
, rnd_mode
);