bugGNU Scientific Library - Bugs: bug #45746, Incorrect results of trigonometric...


You are not allowed to post comments on this tracker with your current authentication level.

bug #45746: Incorrect results of trigonometric functions gsl_sf_sin and gsl_sf_cos

Submitter:  Enyi Tang <eytang>
Submitted:  Thu 13 Aug 2015 06:04:09 AM UTC
Category:  Accuracy problem Severity:  3 - Normal
Operating System:  Status:  None
Assigned to:  None Open/Closed:  Open
Release:  GSL-1.16

Tue 11 Oct 2016 04:35:38 PM UTC, comment #1: 

So gsl_sf_angle_restrict_symm_err_e is basically equivalent to e_rem_pio2.  As you can see in that file, this technique only works until about x ~ 2^18 * M_PI (~823549.6645, so smaller than you think). Past that, you need to implement Payne-Hanek reduction (which is the entirety of k_rem_pio2), which is somewhat of a royal payne.

http://stackoverflow.com/questions/30463616/payne-hanek-algorithm-implementation-in-c is a pretty comprehensive reference. I might throw something together myself later if I have time.

There's really no other alternative other than that or the current implementation, which is "return very invalid results for large values of fabs(x)", sadly.

Peter Barfuss <bofh453>
Thu 13 Aug 2015 06:04:09 AM UTC, original submission:  

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.

Enyi Tang <eytang>


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attached Files
file #34633:  gsltrigtest.c added by eytang (876B - application/octet-stream)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by bofh453 (Posted a comment)
  • -email is unavailable- added by eytang (Submitted the item)

    Follows 1 latest change.

    Date Changed by Updated Field Previous Value => Replaced by
    2015-08-13 eytang Attached File- Added gsltrigtest.c, #34633

    Back to the top

    Powered by Savane 3.13-3736.
    Corresponding source code