1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
|
/* mpfr_fma -- Floating multiply-add
Copyright 2001, 2002, 2004, 2006, 2007 Free Software Foundation, Inc.
Contributed by the Arenaire and Cacao projects, INRIA.
This file is part of the MPFR Library.
The MPFR Library is free software; you can redistribute
it and/or modify it under the terms of the GNU Lesser
General Public License as published by the Free Software
Foundation; either version 2.1 of the License, or (at your
option) any later version.
The MPFR Library is distributed in the hope that it will
be useful, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser
General Public License along with the MPFR Library; see
the file COPYING.LIB. If not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
MA 02110-1301, USA. */
#include "mpfr-impl.h"
/* The fused-multiply-add (fma) of x, y and z is defined by:
fma(x,y,z)= x*y + z
*/
int
mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z,
mp_rnd_t rnd_mode)
{
int inexact;
mpfr_t u;
MPFR_SAVE_EXPO_DECL (expo);
/* particular cases */
if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ||
MPFR_IS_SINGULAR(y) ||
MPFR_IS_SINGULAR(z) ))
{
if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y) || MPFR_IS_NAN(z))
{
MPFR_SET_NAN(s);
MPFR_RET_NAN;
}
/* now neither x, y or z is NaN */
else if (MPFR_IS_INF(x) || MPFR_IS_INF(y))
{
/* cases Inf*0+z, 0*Inf+z, Inf-Inf */
if ((MPFR_IS_ZERO(y)) ||
(MPFR_IS_ZERO(x)) ||
(MPFR_IS_INF(z) &&
((MPFR_MULT_SIGN(MPFR_SIGN(x), MPFR_SIGN(y))) != MPFR_SIGN(z))))
{
MPFR_SET_NAN(s);
MPFR_RET_NAN;
}
else if (MPFR_IS_INF(z)) /* case Inf-Inf already checked above */
{
MPFR_SET_INF(s);
MPFR_SET_SAME_SIGN(s, z);
MPFR_RET(0);
}
else /* z is finite */
{
MPFR_SET_INF(s);
MPFR_SET_SIGN(s, MPFR_MULT_SIGN(MPFR_SIGN(x) , MPFR_SIGN(y)));
MPFR_RET(0);
}
}
/* now x and y are finite */
else if (MPFR_IS_INF(z))
{
MPFR_SET_INF(s);
MPFR_SET_SAME_SIGN(s, z);
MPFR_RET(0);
}
else if (MPFR_IS_ZERO(x) || MPFR_IS_ZERO(y))
{
if (MPFR_IS_ZERO(z))
{
int sign_p;
sign_p = MPFR_MULT_SIGN( MPFR_SIGN(x) , MPFR_SIGN(y) );
MPFR_SET_SIGN(s,(rnd_mode != GMP_RNDD ?
((MPFR_IS_NEG_SIGN(sign_p) && MPFR_IS_NEG(z))
? -1 : 1) :
((MPFR_IS_POS_SIGN(sign_p) && MPFR_IS_POS(z))
? 1 : -1)));
MPFR_SET_ZERO(s);
MPFR_RET(0);
}
else
return mpfr_set (s, z, rnd_mode);
}
else /* necessarily z is zero here */
{
MPFR_ASSERTD(MPFR_IS_ZERO(z));
return mpfr_mul (s, x, y, rnd_mode);
}
}
/* Useless since it is done by mpfr_add
* MPFR_CLEAR_FLAGS(s); */
/* If we take prec(u) >= prec(x) + prec(y), the product u <- x*y
is exact, except in case of overflow or underflow. */
MPFR_SAVE_EXPO_MARK (expo);
mpfr_init2 (u, MPFR_PREC(x) + MPFR_PREC(y));
if (MPFR_UNLIKELY (mpfr_mul (u, x, y, GMP_RNDN)))
{
/* overflow or underflow - this case is regarded as rare, thus
does not need to be very efficient (even if some tests below
could have been done earlier).
It is an overflow iff u is an infinity (since GMP_RNDN was used).
Alternatively, we could test the overflow flag, but in this case,
mpfr_clear_flags would have been necessary. */
if (MPFR_IS_INF (u)) /* overflow */
{
/* Let's eliminate the obvious case where x*y and z have the
same sign. No possible cancellation -> real overflow.
Also, we know that |z| < 2^emax. If E(x) + E(y) >= emax+3,
then |x*y| >= 2^(emax+1), and |x*y + z| >= 2^emax. This case
is also an overflow. */
if (MPFR_SIGN (u) == MPFR_SIGN (z) ||
MPFR_GET_EXP (x) + MPFR_GET_EXP (y) >= __gmpfr_emax + 3)
{
mpfr_clear (u);
MPFR_SAVE_EXPO_FREE (expo);
return mpfr_overflow (s, rnd_mode, MPFR_SIGN (z));
}
/* E(x) + E(y) <= emax+2, therefore |x*y| < 2^(emax+2), and
(x/4)*y does not overflow (let's recall that the result
is exact with an unbounded exponent range). It does not
underflow either because x*y overflows and the exponent
range is large enough. */
inexact = mpfr_div_2ui (u, x, 2, GMP_RNDN);
MPFR_ASSERTN (inexact == 0);
inexact = mpfr_mul (u, u, y, GMP_RNDN);
MPFR_ASSERTN (inexact == 0);
/* Now, we need to add z/4... But it may underflow! */
{
mpfr_t zo4;
mpfr_init2 (zo4, MPFR_PREC (z) + 2);
if (mpfr_div_2ui (zo4, z, 2, GMP_RNDZ))
{
/* The division by 4 underflowed! This probably means that
|z/4| < ulp(u), but this is not guaranteed by the current
MPFR_PREC_MAX definition (and internal computations can
significantly increase the precision).
Let z2 = sign(z) * 2^(E(z)-1), and z4 = z2 + z/4, which
is representable if one takes 2 more precision bits (see
the + 2 above). Then we compute u + z4 with the provided
rounding mode. */
MPFR_ASSERTN (0); /* TODO... */
mpfr_clears (zo4, u, (void *) 0);
}
else
{
/* The division by 4 didn't overflow (and was exact). */
mpfr_clear_flags ();
/* Let's recall that u = x*y/4 and zo4 = z/4 exactly. */
inexact = mpfr_add (s, u, zo4, rnd_mode);
/* u and zo4 have different signs, so that an overflow
is not possible. But an underflow is theoretically
possible! */
if (mpfr_underflow_p ())
{
MPFR_ASSERTN (0); /* TODO... */
mpfr_clears (zo4, u, (void *) 0);
}
else
{
int inex2;
mpfr_clears (zo4, u, (void *) 0);
inex2 = mpfr_mul_2ui (s, s, 2, rnd_mode);
if (inex2) /* overflow */
{
inexact = inex2;
MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
}
MPFR_SAVE_EXPO_FREE (expo);
return mpfr_check_range (s, inexact, rnd_mode);
}
}
}
}
else /* underflow */
{
MPFR_ASSERTN (0); /* TODO... */
mpfr_clear (u);
}
}
inexact = mpfr_add (s, z, u, rnd_mode);
mpfr_clear (u);
MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
MPFR_SAVE_EXPO_FREE (expo);
return mpfr_check_range (s, inexact, rnd_mode);
}
|
