Strange Phenomenon of Cosine Function on Linux x86_64 env
Gordon Sande
The following was in Numerical Analysis Digest. The test program given
is in C but the problem may also apply to Fortran as some Fortrans
use the C runtime library extensively.
=======================================================================
From: Tomonori Kouya <tko...@na-net.ornl.gov>
Date: Sun, 17 Feb 2008 20:45:30 -0500
Subject: Strange Phenomenon of Cosine Function on Linux x86_64 env
Dear NA digest subscribers,
I am in trouble about the problem on cosine function in IEEE754 double
precision. On Fedora Core 4/8 or CentOS 4/5 x86_64 + gcc 4 env, I ran the
program which included the functions below in 4 different rounding-modes (RN,
RZ, RP, RM) in order to know amount of round-off error.
current_rmode = fegetround(); // get current rounding mode
fesetround(rmode); // change rounding mode
ret = cos(x); // cosine function
fesetround(current_rmode); // restore rounding mode
As a result, very different values of cosine func are obtained like:
cos( 5.04710873550435011e+01) =
RN, RZ: 9.78937668119415738e-01, 9.78937668119415627e-01
RP, RM: 5.55012195441045186e-01, 9.78937668119415627e-01
^^^^^^^^^^^^^^^^^^^^^^^ incorrect!
In the RN mode, the cosine function always returns correct values, but often
incorrect and very different values in other 3 modes. If you want to check
this phenomenon on your env, you can download my sample program from
http://na-inet.jp/weblog/archives/test_cos.c. By trying to run the above
"test_cos.c" on 32bit and 64bit Linux environments around me, this phenomenon
could not be found in 32bit env such as Pentium 3/4 and in Mac OS X 10.4 +
Xcode + Core2Duo.
I will appreciate your sending the information about this problem.
Sincerely yours,
Tomonori KOUYA <tko...@na-net.ornl.gov>
Greg Lindahl
>precision. On Fedora Core 4/8 or CentOS 4/5 x86_64 + gcc 4 env, I ran the
>program which included the functions below in 4 different rounding-modes (RN,
>RZ, RP, RM) in order to know amount of round-off error.
I wonder why he assumed that the math library worked correctly in all
rounding modes? Many of them don't. Has anyone seen one that does?
Without penalizing the typical case?
-- greg
FX
> rounding modes? Many of them don't. Has anyone seen one that does?
> Without penalizing the typical case?
CRlibm (Correctly Rounded mathematical library):
http://lipforge.ens-lyon.fr/www/crlibm/
They have benchmarks available, and it is not too bad.
--
FX
Gerry Ford
"Gordon Sande" <g.s...@worldnet.att.net> wrote in message
news:2008022509144616807-gsande@worldnetattnet...
If cos is in trouble, is sin?
--
Gerry Ford
"Er hat sich georgiert." Der Spiegel, 2008, sich auf Chimpy Eins komma null
beziehend.
Pierre Asselin
> From: Tomonori Kouya <tko...@na-net.ornl.gov>
> Date: Sun, 17 Feb 2008 20:45:30 -0500
> Subject: Strange Phenomenon of Cosine Function on Linux x86_64 env
> [ ... ]
> current_rmode = fegetround(); // get current rounding mode
> fesetround(rmode); // change rounding mode
> ret = cos(x); // cosine function
> fesetround(current_rmode); // restore rounding mode
> As a result, very different values of cosine func are obtained like:
> cos( 5.04710873550435011e+01) =
> RN, RZ: 9.78937668119415738e-01, 9.78937668119415627e-01
> RP, RM: 5.55012195441045186e-01, 9.78937668119415627e-01
> ^^^^^^^^^^^^^^^^^^^^^^^ incorrect!
Confirmed here. The problem goes away with any optimization, even
-Os. If you isolate his cos_rmode() function, the generated assembly
looks very different with and without optimization, but they all
call fesetround() around the cosine. I am not fluent enough
in x86_64 assembly to spot the error.
Still, this looks like a compiler bug rather than a libm bug.
--
pa at panix dot com
Steven G. Kargl
p...@see.signature.invalid (Pierre Asselin) writes:
> Gordon Sande <g.s...@worldnet.att.net> wrote:
>
>> From: Tomonori Kouya <tko...@na-net.ornl.gov>
>> Date: Sun, 17 Feb 2008 20:45:30 -0500
>> Subject: Strange Phenomenon of Cosine Function on Linux x86_64 env
>
>> [ ... ]
>> current_rmode = fegetround(); // get current rounding mode
>> fesetround(rmode); // change rounding mode
>> ret = cos(x); // cosine function
>> fesetround(current_rmode); // restore rounding mode
>
>> As a result, very different values of cosine func are obtained like:
>
>> cos( 5.04710873550435011e+01) =
>> RN, RZ: 9.78937668119415738e-01, 9.78937668119415627e-01
>> RP, RM: 5.55012195441045186e-01, 9.78937668119415627e-01
>> ^^^^^^^^^^^^^^^^^^^^^^^ incorrect!
>
> Confirmed here.
Which compiler(s)?
--
Steve