GNU 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: |
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)
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:
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
The exact_result of the Clausen Function above is generated by arbitrary precision arithmetic. |
Enyi Tang <eytang> |
Depends on the following items: None found
Items that depend on this one: None found
Follows 1 latest change.
Date | Changed by | Updated Field | Previous Value | => | Replaced by |
---|---|---|---|---|---|
2015-08-13 | eytang | Attached File | - | ![]() |
Added gsltrigtest.c, #34633 |
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.