[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] [bug #45746] Incorrect results of trigonometric functions gsl_
From: |
Enyi Tang |
Subject: |
[Bug-gsl] [bug #45746] Incorrect results of trigonometric functions gsl_sf_sin and gsl_sf_cos |
Date: |
Thu, 13 Aug 2015 06:04:10 +0000 |
User-agent: |
Mozilla/5.0 (Macintosh; Intel Mac OS X 10_8_3) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/44.0.2403.89 Safari/537.36 |
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/
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Bug-gsl] [bug #45746] Incorrect results of trigonometric functions gsl_sf_sin and gsl_sf_cos,
Enyi Tang <=