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

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Fabian Schriever
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


Reply | Threaded
Open this post in threaded view
|

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

Corinna Vinschen
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 (849 bytes) Download Attachment