Re: [PATCH] Fix error in float trig. function range reduction

Corinna Vinschen <[email protected]> Tue, 3 Mar 2020 16:46:00 +0100

Newsgroups	gmane.comp.lib.newlib
Message-ID	<[email protected]>

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

signature.asc (application/pgp-signature, 833 B)