On Dec 4, 2006, at 4:57 PM, Richard Guenther wrote:
My inclination is that if pow(x, 1./3.) (computed
correctly to the last bit) ever differs from cbrt(x) (computed
correctly to the last bit) then this substitution should not be done.
C99 7.12.7.1 says "The cbrt functions compute the real cube root
of x." and "The cbrt functions return x**1/3". So it looks to me
that cbrt is required to return the same as pow(x, 1/3).
<shrug> For me, this:
#include <math.h>
#include <stdio.h>
int main()
{
printf("pow(1000000., 1./3.) = %a\ncbrt(1000000.) = %a\n",
pow(1000000., 1./3.), cbrt(1000000.));
}
prints out:
pow(1000000., 1./3.) = 0x1.8fffffffffffep+6
cbrt(1000000.) = 0x1.9p+6
I suspect that both are correct, rounded to the nearest least
significant bit. Admittedly I haven't checked the computation by
hand for pow(1000000., 1./3.). But I did the computation using gcc
4.0 on both PPC and Intel Mac hardware, and on PPC using CodeWarrior
math libs, and got the same results on all three platforms.
The pow function is not raising 1,000,000 to the power of 1/3. It is
raising 1,000,000 to some power which is very close to, but not equal
to 1/3.
Perhaps I've misunderstood and you have some way of exactly
representing the fraction 1/3 in Gfortran. In C and C++ we have no
way to exactly represent that fraction except implicitly using cbrt.
(or with a user-defined rational type)
-Howard
