From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>

Date: Wed 29 Jun 2005 - 23:39:08 EST

Date: Wed 29 Jun 2005 - 23:39:08 EST

On Wed, 29 Jun 2005, Duncan Murdoch wrote:

> On 6/29/2005 7:32 AM, Robin Hankin wrote:

>> I have been wondering if there one can speed up calculating small powers

*>> of numbers such as x^8 using multiplication.
**>>
**>> In addition, one can be a bit clever and calculate x^8 using only 3
**>> multiplies.
**>>
**>> look at this:
**>>
**>>
**>> > f1 <- function(x){x*x*x*x*x*x*x*x}
**>> > f2 <- function(x){x^8}
**>> > f3 <- function(x){x2 <- x*x;x4 <- x2*x2;return(x4*x4)}
**>>
**>> [so f1() and f2() and f3() are algebraically identical]
**>>
**>>
**>> > a <- rnorm(1000000)
**>> > system.time(ignore <- f1(a))
**>> [1] 0.50 0.17 2.88 0.00 0.00
**>>
**>> > system.time(ignore <- f2(a))
**>> [1] 0.31 0.03 1.40 0.00 0.00
**>>
**>> > system.time(ignore <- f3(a))
**>> [1] 0.10 0.07 0.18 0.00 0.00
**>>
**>>
**>> [these figures show little variance from trial to trial]
**>>
**>>
**>> I was expecting f2() and f3() to be about the same.
**>> I was not expecting a factor of 3 there!
**>>
**>> anyone got any comments?
**>
**> If you look in src/main/arithmetic.c, you'll see that R does a general
**> real-valued power (using C's pow() function) whenever either one of the
**> args is real (except for a few special cases, e.g. non-numbers, or
**> powers of 2 or 0.5). There is an internal R function equivalent to your
**> f3, but it is not used in the situation of real^integer (and in any
**> case, x^8 is real^real).
**>
**> I suppose if you wanted to submit a patch someone would take a look, but
**> the question is whether there is really any calculation whose execution
**> time would be materially affected by this. Most computations are not
**> dominated by integer power calculations, so is this really worth the
**> trouble?
*

As Luke Tierney frequently points out, selecting special cases can take more time than you save. The assembler code used by modern OSes will have worked out that compromise pretty effectively: even for real^real it spots simple cases of the exponent.

Also, it depends on the CPU: I get on Athlon 2600

*> system.time(ignore <- f1(a), gcFirst=T)
*

[1] 0.08 0.05 0.14 0.00 0.00

*> system.time(ignore <- f2(a), gcFirst=T)
*

[1] 0.20 0.01 0.20 0.00 0.00

*> system.time(ignore <- f3(a), gcFirst=T)
*

[1] 0.03 0.02 0.05 0.00 0.00

and Opteron 248

*> system.time(ignore <- f1(a), gcFirst=T)
*

[1] 0.08 0.06 0.14 0.00 0.00

*> system.time(ignore <- f2(a), gcFirst=T)
*

[1] 0.19 0.01 0.20 0.00 0.00

*> system.time(ignore <- f3(a), gcFirst=T)
*

[1] 0.04 0.01 0.05 0.00 0.00

Note

- the use of gcFirst=T
- these need to be run several times to tune the gc() behaviour. After which f1(a) caused a couple of gc()s and the others none.
- the Opteron is in general a much faster machine so these are all surprisingly slow.

Occasionally the C-level pow() is slow: that happened in one MinGW update and for R Windows users we replaced it. Even there I was unsure that it would make enough difference on a real problem, but eventually found one where it did (MDS with Minkowski distances). It was because of that real problem that I added the special case for 0.5, after careful timing (and hearing Luke's comment alluded to above).

-- Brian D. Ripley, ripley@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595 ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.htmlReceived on Wed Jun 29 23:44:54 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:33:06 EST
*