February 18, 2008

Strange Phenomenon of Cosine function on Linux x86_64 + gcc env

The below article has been submitted to NA Digest (So, it would not be received ...). Japanese version of that can be read here.
----
Dear NA digest subscribers,

I am in trouble about the problem on cosine function in IEEE754 double precision.

On Fedora Core 4 or CentOS 4/5 x86_64 env, I ran the program which included the function below in 4 different rounding-modes (RN, RZ, RP, RM) in order to know amount of round-off error.

double cos_rmode(double x, int rmode)
{
int current_rmode;
double ret;

current_rmode = fegetround(); // get current rounding mode
fesetround(rmode); // change rounding mode
ret = cos(x); // cosine function
fesetround(current_rmode); // restore rounding mode

return ret;
}

As a result, very different values of cosine func at x >> 1 are obtained like:

cos( 5.04710873550435011e+01) =
RN, RZ: 9.78937668119415738e-01, 9.78937668119415627e-01
RP, RM: 5.55012195441045186e-01, 9.78937668119415627e-01

cos( 5.48598775598298900e+01) =
RN, RZ: -1.17720272477898139e-01, -3.04940663255730016e+00
RP, RM: -1.17720272477898139e-01, -3.04940663255730016e+00.

The underlined values above are incorrect.

In the RN mode, the cosine function always returns correct values, but 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 Mac OS X 10.4 + Xcode + Core2Duo.

I will appreciate your sending the information about this problem if you have anything.

Sincerely yours,

--
Tomonori KOUYA
http://na-inet.jp/

Posted by tkouya at February 18, 2008 5:31 PM