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/