Fix error in float trig. function range reduction

Message ID 20200303134906.2011-1-fabian.schriever@gtd-gmbh.de
State Accepted
Commit d4bcecb3e91f73091195adbfdcb0d4287cdf1753
Headers show
Series
  • Fix error in float trig. function range reduction
Related show

Commit Message

Fabian Schriever March 3, 2020, 1:49 p.m.
The single-precision trigonometric functions show rather high errors in
specific ranges starting at about 30000 radians. For example the sinf
procedure produces an error of 7626.55 ULP with the input
5.195880078125e+04 (0x474AF6CD) (compared with MPFR in 128bit
precision). For the test we used 100k values evenly spaced in the range
of [30k, 70k]. The issues are periodic at higher ranges.

This error was introduced when the double precision range reduction was
first converted to float. The shift by 8 bits always returns 0 as iq is
never higher than 255.

The fix reduces the error of the example above to 0.45 ULP, highest
error within the test set fell to 1.31 ULP, which is not perfect, but
still a significant improvement. Testing other previously erroneous
ranges no longer show particularly large accuracy errors.
---
 newlib/libm/math/kf_rem_pio2.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

-- 
2.24.1.windows.2

Comments

Corinna Vinschen March 3, 2020, 3:46 p.m. | #1
On Mar  3 14:49, Fabian Schriever wrote:
> The single-precision trigonometric functions show rather high errors in

> specific ranges starting at about 30000 radians. For example the sinf

> procedure produces an error of 7626.55 ULP with the input

> 5.195880078125e+04 (0x474AF6CD) (compared with MPFR in 128bit

> precision). For the test we used 100k values evenly spaced in the range

> of [30k, 70k]. The issues are periodic at higher ranges.

> 

> This error was introduced when the double precision range reduction was

> first converted to float. The shift by 8 bits always returns 0 as iq is

> never higher than 255.

> 

> The fix reduces the error of the example above to 0.45 ULP, highest

> error within the test set fell to 1.31 ULP, which is not perfect, but

> still a significant improvement. Testing other previously erroneous

> ranges no longer show particularly large accuracy errors.

> ---

>  newlib/libm/math/kf_rem_pio2.c | 2 +-

>  1 file changed, 1 insertion(+), 1 deletion(-)

> 

> diff --git a/newlib/libm/math/kf_rem_pio2.c b/newlib/libm/math/kf_rem_pio2.c

> index 1573ca9f6..3ca7cc20d 100644

> --- a/newlib/libm/math/kf_rem_pio2.c

> +++ b/newlib/libm/math/kf_rem_pio2.c

> @@ -101,7 +101,7 @@ recompute:

>  	    iq[jz-1] -= i<<(8-q0);

>  	    ih = iq[jz-1]>>(7-q0);

>  	} 

> -	else if(q0==0) ih = iq[jz-1]>>8;

> +	else if(q0==0) ih = iq[jz-1]>>7;

>  	else if(z>=(float)0.5) ih=2;

>  

>  	if(ih>0) {	/* q > 0.5 */

> -- 

> 2.24.1.windows.2

> 


Pushed.


Thanks,
Corinna

-- 
Corinna Vinschen
Cygwin Maintainer
Red Hat

Patch

diff --git a/newlib/libm/math/kf_rem_pio2.c b/newlib/libm/math/kf_rem_pio2.c
index 1573ca9f6..3ca7cc20d 100644
--- a/newlib/libm/math/kf_rem_pio2.c
+++ b/newlib/libm/math/kf_rem_pio2.c
@@ -101,7 +101,7 @@  recompute:
 	    iq[jz-1] -= i<<(8-q0);
 	    ih = iq[jz-1]>>(7-q0);
 	} 
-	else if(q0==0) ih = iq[jz-1]>>8;
+	else if(q0==0) ih = iq[jz-1]>>7;
 	else if(z>=(float)0.5) ih=2;
 
 	if(ih>0) {	/* q > 0.5 */