Fortran - Strange Phenomenon of Cosine Function on Linux x86

Author

Strange Phenomenon of Cosine Function on Linux x86_64 env

Gordon Sande

2008-02-25, 8:18 am

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 <tkouya@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 <tkouya@na-net.ornl.gov>

Greg Lindahl

2008-02-25, 7:17 pm

>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.

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

2008-02-25, 7:17 pm

> 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?

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

2008-02-25, 7:17 pm

"Gordon Sande" <g.sande@worldnet.att.net> wrote in message
news:2008022509144616807-gsande@worldnetattnet...
>
> 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 <tkouya@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 <tkouya@na-net.ornl.gov>
>
>
>

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

2008-02-26, 7:23 pm

Gordon Sande <g.sande@worldnet.att.net> wrote:

> From: Tomonori Kouya <tkouya@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

2008-02-26, 7:23 pm

In article <fq18e4$oau$1@reader2.panix.com>,
pa@see.signature.invalid (Pierre Asselin) writes:
> Gordon Sande <g.sande@worldnet.att.net> wrote:
>
>
>
>
>
> Confirmed here.

Which compiler(s)?

--
Steve

Sponsored Links