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: [EMAIL PROTECTED] [mailto:r-help- > [EMAIL PROTECTED] On Behalf Of David Reiner > <[EMAIL PROTECTED]> > 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: [EMAIL PROTECTED] [mailto:r-help- > > [EMAIL PROTECTED] On Behalf Of David Reiner > > <[EMAIL PROTECTED]> > > 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: [EMAIL PROTECTED] [mailto:r-help- > > > [EMAIL PROTECTED] 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