[R] x*x*x*... vs x^n

[davidr@rhotrading.com]

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 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
> >
> > ______________________________________________
> > 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