RepoMirrors/musl
1
0
Fork 0
mirror of git://git.musl-libc.org/musl synced 2024-12-26 00:22:35 +00:00
Code Issues Packages Projects Releases Wiki Activity

math: add fma TODO comments about the underflow issue

Browse Source 
The underflow exception is not raised correctly in some
cornercases (see previous fma commit), added comments
with examples for fmaf, fmal and non-x86 fma.

In fmaf store the result before returning so it has the
correct precision when FLT_EVAL_METHOD!=0
This commit is contained in:
Szabolcs Nagy 2013-05-19 14:43:32 +00:00
parent ffd8ac2dd5
commit 1e5eb73545
3 changed files with 14 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/math/fma.c
View File

@ -441,6 +441,8 @@ double fma(double x, double y, double z)
/*
* There is no need to worry about double rounding in directed
* rounding modes.
* TODO: underflow is not raised properly, example in downward rounding:
* fma(0x1.000000001p-1000, 0x1.000000001p-30, -0x1p-1066)
*/
fesetround(oround);
adj = r.lo + xy.lo;

12
src/math/fmaf.c
View File

@ -49,7 +49,14 @@ float fmaf(float x, float y, float z)
(hr & 0x7ff00000) == 0x7ff00000 || /* NaN */
result - xy == z || /* exact */
fegetround() != FE_TONEAREST) /* not round-to-nearest */
return (result);
{
/*
TODO: underflow is not raised correctly, example in
downward rouding: fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
*/
z = result;
return z;
}
/*
* If result is inexact, and exactly halfway between two float values,
@ -63,5 +70,6 @@ float fmaf(float x, float y, float z)
fesetround(FE_TONEAREST);
if (result == adjusted_result)
SET_LOW_WORD(adjusted_result, lr + 1);
return (adjusted_result);
z = adjusted_result;
return z;
}

2
src/math/fmal.c
View File

@ -262,6 +262,8 @@ long double fmal(long double x, long double y, long double z)
/*
* There is no need to worry about double rounding in directed
* rounding modes.
* TODO: underflow is not raised correctly, example in downward rounding:
* fmal(0x1.0000000001p-16000L, 0x1.0000000001p-400L, -0x1p-16440L)
*/
fesetround(oround);
adj = r.lo + xy.lo;