On Mon, 7 Apr 2008, Robert A LaBudde wrote:
> I'd say there is a problem with an iterative
> algorithm for a fractional power of a number near
> one. I ran your problem in R v.2.6.2 and got:
> ...
> > system.time(ttt.test.sqrt <- ttt.test^0.5)
> user system elapsed
> 0.06 0.01 0.07
>
> I suppose Insightful should look into their algorithm, if it takes 13 minutes.
>
> I would recommend you use "sqrt()" instead of
> "^0.5" as a workaround. It's usually better to do this anyway.
Splus uses that standard math library's pow() function. We do not
write our own code for things like that, although we do have to write
wrappers to work around known problems in the math library's functions
(e.g., they do know about our missing value flags and don't always do
the right things for Inf).
The root cause is that the 64-bit gcc 3.4.5 version of pow(x,y) is
extremely slow when x is 1 bit away from 1 in either direction (1+2^-52
or 1-2^-53) and y is quite close or equal to 0.5. Our tests did not
check those cases.
I just changed the C code in Splus to use sqrt(x) instead of
pow(x,0.5). Later I looked at R's C code (src/main/arithmetic.c) and
noticed that it does the same thing unless one is using gcc 4.3:
/* work around a bug in May 2007 snapshots of gcc pre-4.3.0 */
#if __GNUC__ == 4 && __GNUC_MINOR__ == 3
return (y == 2.0) ? x*x : pow(x, y);
#else
return (y == 2.0) ? x*x : ((y == 0.5) ? sqrt(x) : pow(x, y));
#endif
I don't know if it does this to avoid a speed problem or for some other
reason. The checkin comment says
"work around an optimization bug in currnt gcc pre-4.3.0"
You can still detect the problem in R (compiled with gcc 3.4.5
in 64-bit mode on Linux) if you use a power quite close to 0.5:
R64> unix.time(rep(1+2^-52,1000) ^ 0.5) # uses C sqrt()
user system elapsed
0 0 0
R64> unix.time(rep(1+2^-52,1000) ^ (0.5+2^-46)) # uses C pow()
user system elapsed
2.177 0.000 2.177
> >Allmost 13 minutes to take 826015 square roots of numbers near 1?
----------------------------------------------------------------------------
Bill Dunlap
Insightful Corporation
bill at insightful dot com
360-428-8146
"All statements in this message represent the opinions of the author and do
not necessarily reflect Insightful Corporation policy or position."
|