From: <davidr_at_rhotrading.com>

Date: Thu 30 Jun 2005 - 04:55:10 EST

> 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

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.html Received on Thu Jun 30 05:01:04 2005

Date: Thu 30 Jun 2005 - 04:55:10 EST

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@stat.math.ethz.ch [mailto:r-help-
**> bounces@stat.math.ethz.ch] On Behalf Of David Reiner
*

> <davidr@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@stat.math.ethz.ch [mailto:r-help-
**> > bounces@stat.math.ethz.ch] On Behalf Of David Reiner
**> > <davidr@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@stat.math.ethz.ch [mailto:r-help-
**> > > bounces@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@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@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@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@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@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 Received on Thu Jun 30 05:01:04 2005

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