git
/
gcc
mirror of git://gcc.gnu.org/git/gcc.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

fmaq.c (fmaq): Merge from GLIBC. 

2012-11-15  Tobias Burnus  <burnus@net-b.de>
            Joseph Myers  <joseph@codesourcery.com>

        * math/fmaq.c (fmaq): Merge from GLIBC. Fix fma
        underflows with small x * y; Fix overflow results
        outside round-to-nearest mode; make use of Dekker
        and Knuth algorithms use round-to-nearest.


Co-Authored-By: Joseph Myers <joseph@codesourcery.com>

From-SVN: r193538
This commit is contained in:
Tobias Burnus 2012-11-15 18:22:21 +01:00 committed by Tobias Burnus
parent 0c604a61a3
commit 7ee2eb8277
2 changed files with 52 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
libquadmath/ChangeLog
View File

@ -1,3 +1,11 @@
2012-11-15 Tobias Burnus <burnus@net-b.de>
Joseph Myers <joseph@codesourcery.com>
* math/fmaq.c (fmaq): Merge from GLIBC. Fix fma
underflows with small x * y; Fix overflow results
outside round-to-nearest mode; make use of Dekker
and Knuth algorithms use round-to-nearest.
2012-11-01 Tobias Burnus <burnus@net-b.de> 2012-11-01 Tobias Burnus <burnus@net-b.de>
* math/fmaq.c (fmaq): Fix build. * math/fmaq.c (fmaq): Fix build.

58
libquadmath/math/fmaq.c
View File

@ -14,9 +14,8 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
#include <math.h> #include <math.h>
@ -62,17 +61,18 @@ fmaq (__float128 x, __float128 y, __float128 z)
underflows to 0. */ underflows to 0. */
if (z == 0 && x != 0 && y != 0) if (z == 0 && x != 0 && y != 0)
return x * y; return x * y;
/* If x or y or z is Inf/NaN, or if fma will certainly overflow, /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
or if x * y is less than half of FLT128_DENORM_MIN, x * y + z. */
compute as x * y + z. */
if (u.ieee.exponent == 0x7fff if (u.ieee.exponent == 0x7fff
|| v.ieee.exponent == 0x7fff || v.ieee.exponent == 0x7fff
|| w.ieee.exponent == 0x7fff || w.ieee.exponent == 0x7fff
|| u.ieee.exponent + v.ieee.exponent || x == 0
> 0x7fff + IEEE854_FLOAT128_BIAS || y == 0)
|| u.ieee.exponent + v.ieee.exponent
< IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG - 2)
return x * y + z; return x * y + z;
/* If fma will certainly overflow, compute as x * y. */
if (u.ieee.exponent + v.ieee.exponent
> 0x7fff + IEEE854_FLOAT128_BIAS)
return x * y;
/* If x * y is less than 1/4 of FLT128_DENORM_MIN, neither the /* If x * y is less than 1/4 of FLT128_DENORM_MIN, neither the
result nor whether there is underflow depends on its exact result nor whether there is underflow depends on its exact
value, only on its sign. */ value, only on its sign. */
@ -121,8 +121,17 @@ fmaq (__float128 x, __float128 y, __float128 z)
{ {
/* Similarly. /* Similarly.
If z exponent is very large and x and y exponents are If z exponent is very large and x and y exponents are
very small, it doesn't matter if we don't adjust it. */ very small, adjust them up to avoid spurious underflows,
if (u.ieee.exponent > v.ieee.exponent) rather than down. */
if (u.ieee.exponent + v.ieee.exponent
<= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG)
{
if (u.ieee.exponent > v.ieee.exponent)
u.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
else
v.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
}
else if (u.ieee.exponent > v.ieee.exponent)
{ {
if (u.ieee.exponent > FLT128_MANT_DIG) if (u.ieee.exponent > FLT128_MANT_DIG)
u.ieee.exponent -= FLT128_MANT_DIG; u.ieee.exponent -= FLT128_MANT_DIG;
@ -175,6 +184,12 @@ fmaq (__float128 x, __float128 y, __float128 z)
if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0)) if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0))
return x * y + z; return x * y + z;
#ifdef USE_FENV_H
fenv_t env;
feholdexcept (&env);
fesetround (FE_TONEAREST);
#endif
/* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */
#define C ((1LL << (FLT128_MANT_DIG + 1) / 2) + 1) #define C ((1LL << (FLT128_MANT_DIG + 1) / 2) + 1)
__float128 x1 = x * C; __float128 x1 = x * C;
@ -193,10 +208,25 @@ fmaq (__float128 x, __float128 y, __float128 z)
t1 = m1 - t1; t1 = m1 - t1;
t2 = z - t2; t2 = z - t2;
__float128 a2 = t1 + t2; __float128 a2 = t1 + t2;
#ifdef USE_FENV_H
feclearexcept (FE_INEXACT);
#endif
/* If the result is an exact zero, ensure it has the correct
sign. */
if (a1 == 0 && m2 == 0)
{
#ifdef USE_FENV_H
feupdateenv (&env);
#endif
/* Ensure that round-to-nearest value of z + m1 is not
reused. */
asm volatile ("" : "=m" (z) : "m" (z));
return z + m1;
}
#ifdef USE_FENV_H #ifdef USE_FENV_H
fenv_t env;
feholdexcept (&env);
fesetround (FE_TOWARDZERO); fesetround (FE_TOWARDZERO);
#endif #endif
/* Perform m2 + a2 addition with round to odd. */ /* Perform m2 + a2 addition with round to odd. */