From d4bcecb3e91f73091195adbfdcb0d4287cdf1753 Mon Sep 17 00:00:00 2001 From: Fabian Schriever Date: Tue, 3 Mar 2020 14:49:06 +0100 Subject: Fix error in float trig. function range reduction 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 1573ca9..3ca7cc2 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 */ -- cgit v1.1