From: Alberto Monteiro <albmont_at_centroin.com.br>

Date: Wed, 07 May 2008 09:49:56 -0200

}

M <- M %*% M

n <- n / 2

}

return(result)

}

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Wed 07 May 2008 - 12:27:41 GMT

Date: Wed, 07 May 2008 09:49:56 -0200

Jorge Ivan Velez wrote:

*>
*

> I think the function could be better but try this:

*>
**> # Function: M is your matrix and n MUST be an integer>0
**> mat.pow<-function(M,n) {
**> result<-M
**> if(n>1){
**> for ( iter in 2:n) result<-M%*%result
**> result
**> }
**> else {result}
**> result
**> }
**>
*

There are much more efficient ways to compute a power,
just check:

http://en.wikipedia.org/wiki/Exponentiation_by_squaring

Or, translating the pseudo-code to R:

mat.pow <- function(M, n) {

result <- diag(1, ncol(M))

while (n > 0) {

if (n %% 2 == 1) {

result <- result %*% M n <- n - 1

}

M <- M %*% M

n <- n / 2

}

return(result)

}

# The matrix

m <- rbind(c(1,1,0), c(0,1,1), c(0,0,1))

# Goal m^6 = m x m x m x m x m x m

goal= m %*% m %*% m %*% m %*% m %*% m

# matpow

res=mat.pow(m,6)

# Check point

all.equal(goal,res)

This algorithm would be fast, unless n is a _very_ big number.

Alberto Monteiro

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Wed 07 May 2008 - 12:27:41 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Wed 07 May 2008 - 12:30:35 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*