[PATCH] Fix error in exp in magnitude [2e-32,2e-28]

Fabian Schriever <[email protected]>
Newsgroups gmane.comp.lib.newlib
Message-ID <[email protected]>
While testing the exp function we noticed some errors at the specified
magnitude. Within this range the exp function returns the input value +1
as an output. We chose to run a test of 1m exponentially spaced values
in the ranges [-2^-27,-2^-32] and [2^-32,2^-27] which showed 7603 and
3912 results with an error of >=0.5 ULP (compared with MPFR in 128 bit)
with the highest being 0.56 ULP and 0.53 ULP.

It's easy to fix by changing the magnitude at which the input value +1
is returned from <2^-28 to <2^-32 and using the polynomial instead. This
reduces the number of results with an error of >=0.5 ULP to 485 and 479
in above tests, all of which are exactly 0.5 ULP.

As we were already checking on exp we also took a look at expf. For expf
the magnitude where the input value +1 is returned can be increased from
<2^-28 to <2^-23 without accuracy loss for a slight performance
improvement. To ensure this was the correct value we tested all values
in the ranges [-2^-17,-2^-28] and [2^-28,2^-17] (~92.3m values each).
---
 newlib/libm/math/e_exp.c  | 2 +-
 newlib/libm/math/ef_exp.c | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/newlib/libm/math/e_exp.c b/newlib/libm/math/e_exp.c
index 81ea64dfb..d23b1162b 100644
--- a/newlib/libm/math/e_exp.c
+++ b/newlib/libm/math/e_exp.c
@@ -142,7 +142,7 @@ P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
 	    }
 	    x  = hi - lo;
 	} 
-	else if(hx < 0x3e300000)  {	/* when |x|<2**-28 */
+	else if(hx < 0x3df00000)  {	/* when |x|<2**-32 */
 	    if(huge+x>one) return one+x;/* trigger inexact */
 	}
 
diff --git a/newlib/libm/math/ef_exp.c b/newlib/libm/math/ef_exp.c
index e817370ac..fb3e2ffe6 100644
--- a/newlib/libm/math/ef_exp.c
+++ b/newlib/libm/math/ef_exp.c
@@ -77,7 +77,7 @@ P5   =  4.1381369442e-08; /* 0x3331bb4c */
 	    }
 	    x  = hi - lo;
 	} 
-	else if(hx < 0x31800000)  {	/* when |x|<2**-28 */
+	else if(hx < 0x34000000)  {	/* when |x|<2**-23 */
 	    if(huge+x>one) return one+x;/* trigger inexact */
 	}
 
-- 
2.24.1.windows.2
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.