On 12/4/06, Howard Hinnant <[EMAIL PROTECTED]> wrote:
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)
1./3. is represented (round-to-nearest) as 0x1.5555555555555p-2 and pow (x, 1./3.) is of course not the same as the cube-root of x with exact arithmetic. cbrt as defined by C99 suggests that an approximation by pow (x, 1./3.) fullfils the requirements. The question is whether a correctly rounded "exact" cbrt differs from the pow replacement by more than 1ulp - it looks like this is not the case. For cbrt (100.) I get the same result as from pow (100., nextafter (1./3., 1)) for example. So, instead of only recognizing 1./3. rounded to nearest we should recognize both representable numbers that are nearest to 1./3.. Richard.
