Hi
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?
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
Something like this is exploited very nicely in the mtx.exp
for matrix powers in the Malmig package, actually.
****************************
Hi
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?
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
____________________
Ken Knoblauch
Inserm U371, Cerveau et Vision
Department of Cognitive Neurosciences
18 avenue du Doyen Lepine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: 06 84 10 64 10
http://www.lyon.inserm.fr/371/
I tried your code and got different results:
system.time(ignore <- f1(a))
[1] 0.83 0.09 1.08 NA NA
> system.time(ignore <- f2(a))
[1] 0.38 0.01 0.41 NA NA
> system.time(ignore <- f3(a))
[1] 0.32 0.04 0.43 NA NA
So I tried it again but with a loop and got:
> for(i in 1:10) cat(system.time(ignore <- f2(a)), "\n")
0.36 0.04 0.44 NA NA
0.32 0.01 0.34 NA NA
0.28 0.03 0.32 NA NA
0.29 0.03 0.35 NA NA
0.3 0.02 0.32 NA NA
0.28 0.03 0.32 NA NA
0.3 0.02 0.32 NA NA
0.29 0.02 0.34 NA NA
0.23 0.03 0.32 NA NA
0.42 0 0.45 NA NA
> for(i in 1:10) cat(system.time(ignore <- f3(a)), "\n")
0.19 0.04 0.25 NA NA
0.17 0.04 0.25 NA NA
0.21 0.02 0.25 NA NA
0.21 0.02 0.23 NA NA
0.18 0.04 0.23 NA NA
0.18 0.05 0.23 NA NA
0.18 0.04 0.25 NA NA
0.17 0.06 0.23 NA NA
0.2 0.02 0.23 NA NA
0.14 0.06 0.25 NA NA
It seems to me that f3 is 50% slower than f2 not 300%.
My System is:
- R version: R 2.1.1
- Operating System: Win XP
- Compiler: mingw32-gcc-3.4.2
Jarek
====================================================\======
Jarek Tuszynski, PhD. o / \
Science Applications International Corporation <\__,|
(703) 676-4192 "> \
Jaroslaw.W.Tuszynski at saic.com ` \
-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Robin Hankin
Sent: Wednesday, June 29, 2005 7:32 AM
To: r-help
Subject: [R] x*x*x*... vs x^n
Hi
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?
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton European Way, Southampton SO14
3ZH, UK
tel 023-8059-7743
______________________________________________
R-help at 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.html
On 6/29/2005 7:32 AM, Robin Hankin wrote:> Hi > > 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? Duncan Murdoch
In general, the "Russian peasant algorithm", which requires only O(log
n) multiplications, is very good. Section 4.6.3 of Knuth's The Art of
Computer Programming. Volume 2: Seminumerical Algorithms has an in depth
discussion.
I have had to use this in the past, when computers were slower and
compilers were not so clever. It is also better when x is not just a
real number, say complex or matrix (as has been mentioned.)
In many cases though, one needs many powers sequentially, and then it
may be more efficient to just multiply the previous power by x and use
the power, etc. (unless you have a parallel computer.)
So
pows <- x^(1:1000)
# use pows in calculations
could be sped up by employing a faster algorithm, but probably a loop
will be faster:
pows <- 1
for(i in 1:1000) {
pows <- pows * x
# use this power
}
David L. Reiner, Ph.D.
Rho Trading
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-
> bounces at stat.math.ethz.ch] On Behalf Of Robin Hankin
> Sent: Wednesday, June 29, 2005 9:13 AM
> To: Duncan Murdoch
> Cc: r-help; Robin Hankin
> Subject: Re: [R] x*x*x*... vs x^n
>
>
> On Jun 29, 2005, at 02:47 pm, Duncan Murdoch wrote:
>
> > On 6/29/2005 9:31 AM, Robin Hankin wrote:
> >
> >> Hi Duncan
> >>
> >>
> >> library(gsl)
> >> system.time(ignore <- pow_int(a,8))
> >> [1] 1.07 1.11 3.08 0.00 0.00
> >>
> >> <why the slow execution time?>
> >>
> >
> > Shouldn't you ask the gsl maintainer that? :-)
> >
>
> well I did ask myself, but in this case the gsl maintainer
> told me to ask the gsl co-author, who
> is no doubt much better informed in these matters ;-)
>
> >>
> >> (Of course, I'm not suggesting that other programming tasks be
> >> suspended! All I'm pointing
> >> out is that there may exist a user to whom fast integer powers are
> >> very very important)
> >>
> >
> > Then that user should submit the patch, not you. But whoever does
it> > should include an argument to convince an R core member that the
> > change
> > is worth looking at, and do it well enough that the patch is
accepted.> > Changes to the way R does arithmetic affect everyone, so they had
> > better
> > be done right, and checking them takes time.
> >
>
> yes, that's a fair point.
> But including a native R command pow.int(), say, wouldn't affect
> anyone, would it?
> One could even use the (tested) GSL code, as it is GPL'ed.
>
> This would just be a new function that users could use at their
> discretion, and no
> existing code would break.
>
> I assume that such a function would not suffer whatever performance
> disadvantage
> that the GSL package approach had, so it may well be quite a
> significant improvement
> over the method used by R_pow_di() in arithmetic.c especially for
> large n.
>
>
> best wishes
>
> rksh
>
> > Duncan Murdoch
> >
> > ______________________________________________
> > R-help at 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.html
> >
>
> ______________________________________________
> R-help at 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.html
Looking at the code for gsl_pow_int, I see they do use that method. David L. Reiner> -----Original Message----- > From: r-help-bounces at stat.math.ethz.ch [mailto:r-help- > bounces at stat.math.ethz.ch] On Behalf Of David Reiner > <davidr at rhotrading.com> > Sent: Wednesday, June 29, 2005 9:50 AM > To: r-help > Subject: Re: [R] x*x*x*... vs x^n > > In general, the "Russian peasant algorithm", which requires only O(log > n) multiplications, is very good. Section 4.6.3 of Knuth's The Art of > Computer Programming. Volume 2: Seminumerical Algorithms has an indepth> discussion. > I have had to use this in the past, when computers were slower and > compilers were not so clever. It is also better when x is not just a > real number, say complex or matrix (as has been mentioned.) > In many cases though, one needs many powers sequentially, and then it > may be more efficient to just multiply the previous power by x and use > the power, etc. (unless you have a parallel computer.) > So > > pows <- x^(1:1000) > # use pows in calculations > > could be sped up by employing a faster algorithm, but probably a loop > will be faster: > > pows <- 1 > for(i in 1:1000) { > pows <- pows * x > # use this power > } > > David L. Reiner, Ph.D. > Rho Trading > > > > -----Original Message----- > > From: r-help-bounces at stat.math.ethz.ch [mailto:r-help- > > bounces at stat.math.ethz.ch] On Behalf Of Robin Hankin > > Sent: Wednesday, June 29, 2005 9:13 AM > > To: Duncan Murdoch > > Cc: r-help; Robin Hankin > > Subject: Re: [R] x*x*x*... vs x^n > > > > > > On Jun 29, 2005, at 02:47 pm, Duncan Murdoch wrote: > > > > > On 6/29/2005 9:31 AM, Robin Hankin wrote: > > > > > >> Hi Duncan > > >> > > >> > > >> library(gsl) > > >> system.time(ignore <- pow_int(a,8)) > > >> [1] 1.07 1.11 3.08 0.00 0.00 > > >> > > >> <why the slow execution time?> > > >> > > > > > > Shouldn't you ask the gsl maintainer that? :-) > > > > > > > well I did ask myself, but in this case the gsl maintainer > > told me to ask the gsl co-author, who > > is no doubt much better informed in these matters ;-) > > > > >> > > >> (Of course, I'm not suggesting that other programming tasks be > > >> suspended! All I'm pointing > > >> out is that there may exist a user to whom fast integer powersare> > >> very very important) > > >> > > > > > > Then that user should submit the patch, not you. But whoever does > it > > > should include an argument to convince an R core member that the > > > change > > > is worth looking at, and do it well enough that the patch is > accepted. > > > Changes to the way R does arithmetic affect everyone, so they had > > > better > > > be done right, and checking them takes time. > > > > > > > yes, that's a fair point. > > But including a native R command pow.int(), say, wouldn't affect > > anyone, would it? > > One could even use the (tested) GSL code, as it is GPL'ed. > > > > This would just be a new function that users could use at their > > discretion, and no > > existing code would break. > > > > I assume that such a function would not suffer whatever performance > > disadvantage > > that the GSL package approach had, so it may well be quite a > > significant improvement > > over the method used by R_pow_di() in arithmetic.c especially for > > large n. > > > > > > best wishes > > > > rksh > > > > > Duncan Murdoch > > > > > > ______________________________________________ > > > R-help at 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.html > > > > > > > ______________________________________________ > > R-help at 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.html > > ______________________________________________ > R-help at 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.html
I was surprised to find that I was wrong about powers of complexes:> seq.pow1 <- function(x,n) {+ y <- rep(x,n) + for(i in 2:n) y[i] <- y[i-1] * x + y + }> seq.pow2 <- function(x,n) x^(1:n) > x <- 1.001 + 1i * 0.999# several reps of the following> system.time(ignore <- seq.pow1(x,100000),gcFirst=TRUE)[1] 0.73 0.00 0.74 NA NA # several reps of the following> system.time(ignore <- seq.pow2(x,100000),gcFirst=TRUE)[1] 0.35 0.00 0.35 NA NA I apologize for using "probably" below when "not" was correct (modulo grammar). David L. Reiner> -----Original Message----- > From: r-help-bounces at stat.math.ethz.ch [mailto:r-help- > bounces at stat.math.ethz.ch] On Behalf Of David Reiner > <davidr at rhotrading.com> > Sent: Wednesday, June 29, 2005 10:24 AM > To: r-help > Subject: Re: [R] x*x*x*... vs x^n > > Looking at the code for gsl_pow_int, I see they do use that method. > > David L. Reiner > > > > -----Original Message----- > > From: r-help-bounces at stat.math.ethz.ch [mailto:r-help- > > bounces at stat.math.ethz.ch] On Behalf Of David Reiner > > <davidr at rhotrading.com> > > Sent: Wednesday, June 29, 2005 9:50 AM > > To: r-help > > Subject: Re: [R] x*x*x*... vs x^n > > > > In general, the "Russian peasant algorithm", which requires onlyO(log> > n) multiplications, is very good. Section 4.6.3 of Knuth's The Artof> > Computer Programming. Volume 2: Seminumerical Algorithms has an in > depth > > discussion. > > I have had to use this in the past, when computers were slower and > > compilers were not so clever. It is also better when x is not just a > > real number, say complex or matrix (as has been mentioned.) > > In many cases though, one needs many powers sequentially, and thenit> > may be more efficient to just multiply the previous power by x anduse> > the power, etc. (unless you have a parallel computer.) > > So > > > > pows <- x^(1:1000) > > # use pows in calculations > > > > could be sped up by employing a faster algorithm, but probably aloop> > will be faster: > > > > pows <- 1 > > for(i in 1:1000) { > > pows <- pows * x > > # use this power > > } > > > > David L. Reiner, Ph.D. > > Rho Trading > > > > > > > -----Original Message----- > > > From: r-help-bounces at stat.math.ethz.ch [mailto:r-help- > > > bounces at stat.math.ethz.ch] On Behalf Of Robin Hankin > > > Sent: Wednesday, June 29, 2005 9:13 AM > > > To: Duncan Murdoch > > > Cc: r-help; Robin Hankin > > > Subject: Re: [R] x*x*x*... vs x^n > > > > > > > > > On Jun 29, 2005, at 02:47 pm, Duncan Murdoch wrote: > > > > > > > On 6/29/2005 9:31 AM, Robin Hankin wrote: > > > > > > > >> Hi Duncan > > > >> > > > >> > > > >> library(gsl) > > > >> system.time(ignore <- pow_int(a,8)) > > > >> [1] 1.07 1.11 3.08 0.00 0.00 > > > >> > > > >> <why the slow execution time?> > > > >> > > > > > > > > Shouldn't you ask the gsl maintainer that? :-) > > > > > > > > > > well I did ask myself, but in this case the gsl maintainer > > > told me to ask the gsl co-author, who > > > is no doubt much better informed in these matters ;-) > > > > > > >> > > > >> (Of course, I'm not suggesting that other programming tasks be > > > >> suspended! All I'm pointing > > > >> out is that there may exist a user to whom fast integer powers > are > > > >> very very important) > > > >> > > > > > > > > Then that user should submit the patch, not you. But whoeverdoes> > it > > > > should include an argument to convince an R core member that the > > > > change > > > > is worth looking at, and do it well enough that the patch is > > accepted. > > > > Changes to the way R does arithmetic affect everyone, so theyhad> > > > better > > > > be done right, and checking them takes time. > > > > > > > > > > yes, that's a fair point. > > > But including a native R command pow.int(), say, wouldn't affect > > > anyone, would it? > > > One could even use the (tested) GSL code, as it is GPL'ed. > > > > > > This would just be a new function that users could use at their > > > discretion, and no > > > existing code would break. > > > > > > I assume that such a function would not suffer whateverperformance> > > disadvantage > > > that the GSL package approach had, so it may well be quite a > > > significant improvement > > > over the method used by R_pow_di() in arithmetic.c especially for > > > large n. > > > > > > > > > best wishes > > > > > > rksh > > > > > > > Duncan Murdoch > > > > > > > > ______________________________________________ > > > > R-help at 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.html > > > > > > > > > > ______________________________________________ > > > R-help at 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.html > > > > ______________________________________________ > > R-help at 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.html > > ______________________________________________ > R-help at 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.html