[Bug-gsl] [bug #45746] Incorrect results of trigonometric functions gsl_sf_sin and gsl_sf_cos

[Enyi Tang]

URL:
  <http://savannah.gnu.org/bugs/?45746>

                 Summary: Incorrect results of trigonometric functions
gsl_sf_sin and gsl_sf_cos
                 Project: GNU Scientific Library
            Submitted by: eytang
            Submitted on: Thu 13 Aug 2015 06:04:09 AM GMT
                Category: Accuracy problem
                Severity: 3 - Normal
        Operating System: 
                  Status: None
             Assigned to: None
             Open/Closed: Open
                 Release: GSL-1.16
         Discussion Lock: Any

    _______________________________________________________

Details:

When I further investigate the reason of bug #45726, and the bug #45730. I
find the problem is from trigonometric function gsl_sf_sin and gsl_sf_cos. It
is cause by the cancellation in specfunc/trig.c, such as line 183 and line 199
in the following: 

183      double y = floor(abs_x/(0.25*M_PI));         // y=(int)|x|/(pi/4)
         …
199      z = ((abs_x - y * P1) - y * P2) - y * P3;    // z=|x|-y*(pi/4) In
theory 0 <= z <= (pi/4).

I comment the semantics of both the statements. Here when abs_x is a large
number such as 1.0e+22, line 183 cannot get the exact value of y. So y may be
a little bigger or a little smaller than its theoretical value. Then z will
not hold its theoretical properties. In my  testing, I found z can be
-6.53832e+03 in the execution.

Suggest fixing: I think the function gsl_sf_sin_e, gsl_sf_cos_e,
gsl_sf_angle_restrict_symm_err_e, gsl_sf_angle_restrict_pos_err_e should add a
branch as:
if(abs_x >= 0.00001/GSL_DBL_EPSILON) {
  // increase the precision of intermediate values, such as: y, pi, z etc.
  // use arbitrary precision arithmetic instead of double.
  // details can reference the source code of FdLibM, function
__ieee754_rem_pio2
  // http://www.netlib.org/fdlibm/e_rem_pio2.c
  // http://www.netlib.org/fdlibm/k_rem_pio2.c
  // http://www.netlib.org/fdlibm/s_sin.c
}

And I believe the problem also lies in the function
gsl_sf_angle_restrict_symm_err_e and gsl_sf_angle_restrict_pos_err_e, which
affects the precision of other functions, such as gsl_sf_clausen. A test case
is attached running on the latest version
867624b55b20de8da80d23d90549c74ec24cb3a6 on Aug. 7, 2015 in git. 


gsl_sf_sin(-1.000000e+22):
1835521230014198724225866739463153983718581136987376689610752.000000 
gsl_sf_cos(-1.000000e+22):
-5697492148800642151146773364290808164095751450412938747183104.000000 
gsl_sf_sin(-1.000000e+21):
-92157107960277830475099544379269515075741304051510280192.000000 
gsl_sf_cos(-1.000000e+21):
-286162998538776391616589535986753232120837467623860994048.000000 
gsl_sf_sin(-1.000000e+20): -915544002183042865166003045245607271727104.000000

gsl_sf_cos(-1.000000e+20): -2861662284539704961694391618394466579120128.000000

gsl_sf_sin(-1.000000e+19): 28883823730260844544.000000 
gsl_sf_cos(-1.000000e+19): -124644955674068746240.000000 
gsl_sf_sin(-1.000000e+18): 61.480345 
gsl_sf_cos(-1.000000e+18): -68.613184 
gsl_sf_sin(-1.000000e+17): 534.439930 
gsl_sf_cos(-1.000000e+17): -683.039086 
gsl_sf_sin(1.000000e+17): -534.439930 
gsl_sf_cos(1.000000e+17): -683.039086 
gsl_sf_sin(1.000000e+18): -61.480345 
gsl_sf_cos(1.000000e+18): -68.613184 
gsl_sf_sin(1.000000e+19): -28883823730260844544.000000 
gsl_sf_cos(1.000000e+19): -124644955674068746240.000000 
gsl_sf_sin(1.000000e+20): 915544002183042865166003045245607271727104.000000 
gsl_sf_cos(1.000000e+20): -2861662284539704961694391618394466579120128.000000

gsl_sf_sin(1.000000e+21):
92157107960277830475099544379269515075741304051510280192.000000 
gsl_sf_cos(1.000000e+21):
-286162998538776391616589535986753232120837467623860994048.000000 
gsl_sf_sin(1.000000e+22):
-1835521230014198724225866739463153983718581136987376689610752.000000 
gsl_sf_cos(1.000000e+22):
-5697492148800642151146773364290808164095751450412938747183104.000000 
gsl_sf_clausen(1.000000e+11): 1.006468684791068 exact_result: -1.0064688
gsl_sf_clausen(1.000000e+12): -0.937211344935619 exact_result: 0.9372076
gsl_sf_clausen(5.000000e+12): 0.101263304837350 exact_result: -0.1014365
gsl_sf_clausen(1.000000e+13): -0.652489993669150 exact_result: 0.6531079
gsl_sf_clausen(1.000000e+14): -0.150520089031612 exact_result: 0.1458414
gsl_sf_clausen(2.000000e+14): 0.800325512438204 exact_result: -0.7870726

The exact_result of the Clausen Function above is generated by arbitrary
precision arithmetic. 



    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Thu 13 Aug 2015 06:04:09 AM GMT  Name: gsltrigtest.c  Size: 876B   By:
eytang

<http://savannah.gnu.org/bugs/download.php?file_id=34633>

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?45746>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/