Matrix-exponentiation operator: Difference between revisions

Content added Content deleted
m (→‎Infix operator: Improved syntax.)
Line 2,869: Line 2,869:
===Infix operator===
===Infix operator===
The task wants the implementation to be "as an operator". Given that R lets us define new infix operators, it seems fitting to show how to do this. Ideally, for a matrix a and int n, we'd want to be able to use a^n. R actually has this already, but it's not what the task wants:
The task wants the implementation to be "as an operator". Given that R lets us define new infix operators, it seems fitting to show how to do this. Ideally, for a matrix a and int n, we'd want to be able to use a^n. R actually has this already, but it's not what the task wants:
<lang rsplus>a<-matrix(c(1,2,3,4),2,2)
<lang rsplus>a <- matrix(c(1, 2, 3, 4), 2, 2)
a^1
a^1
a^2</lang>
a^2</lang>
Line 2,882: Line 2,882:
[2,] 4 16</pre>
[2,] 4 16</pre>
As we can see, it instead returns the given matrix with its elements raised to the nth power. Overwriting the ^ operator would be dangerous and rude. However, R's base library suggests an alternative. %*% is already defined as matrix multiplication, so why not use %^% for exponentiation?
As we can see, it instead returns the given matrix with its elements raised to the nth power. Overwriting the ^ operator would be dangerous and rude. However, R's base library suggests an alternative. %*% is already defined as matrix multiplication, so why not use %^% for exponentiation?
<lang rsplus>`%^%`<-function(mat,n)
<lang rsplus>`%^%` <- function(mat, n)
{
{
is.wholenumber<-function(x,tol=.Machine$double.eps^0.5){abs(x - round(x))<tol}#See the docs for is.integer
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol#See the docs for is.integer
if(is.matrix(mat) && is.numeric(n) && is.wholenumber(n))
if(is.matrix(mat) && is.numeric(n) && is.wholenumber(n))
{
{
if(n==0){diag(nrow=nrow(mat))}#Identity matrix of mat's dimensions
if(n==0) diag(nrow = nrow(mat))#Identity matrix of mat's dimensions
else if(n==1){mat}
else if(n == 1) mat
else if(n>1){mat%*%(mat%^%(n-1))}
else if(n > 1) mat %*% (mat%^%(n - 1))
else stop("Invalid n.")
else stop("Invalid n.")
}
}
Line 2,895: Line 2,895:
}
}
#For output:
#For output:
a%^%0
a %^% 0
a%^%1
a %^% 1
a%^%2
a %^% 2
a%*%a%*%a#Base R's equivalent of a%^%3
a %*% a %*% a#Base R's equivalent of a %^% 3
a%^%3
a %^% 3
nonSquareMatrix<-matrix(c(1,2,3,4,5,6),nrow=2,ncol=3)
nonSquareMatrix <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3)
nonSquareMatrix%^%1
nonSquareMatrix %^% 1
nonSquareMatrix%^%2#R's %*% will throw the error for us</lang>
nonSquareMatrix %^% 2#R's %*% will throw the error for us</lang>
{{out}}
{{out}}
<pre>> a%^%0
<pre>> a %^% 0
[,1] [,2]
[,1] [,2]
[1,] 1 0
[1,] 1 0
[2,] 0 1
[2,] 0 1


> a%^%1
> a %^% 1
[,1] [,2]
[,1] [,2]
[1,] 1 3
[1,] 1 3
[2,] 2 4
[2,] 2 4


> a%^%2
> a %^% 2
[,1] [,2]
[,1] [,2]
[1,] 7 15
[1,] 7 15
[2,] 10 22
[2,] 10 22


> a%*%a%*%a#Base R's equivalent of a%^%3
> a %*% a %*% a#Base R's equivalent of a %^% 3
[,1] [,2]
[,1] [,2]
[1,] 37 81
[1,] 37 81
[2,] 54 118
[2,] 54 118


> a%^%3
> a %^% 3
[,1] [,2]
[,1] [,2]
[1,] 37 81
[1,] 37 81
[2,] 54 118
[2,] 54 118


> nonSquareMatrix<-matrix(c(1,2,3,4,5,6),nrow=2,ncol=3)
> nonSquareMatrix <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3)


> nonSquareMatrix%^%1
> nonSquareMatrix %^% 1
[,1] [,2] [,3]
[,1] [,2] [,3]
[1,] 1 3 5
[1,] 1 3 5
[2,] 2 4 6
[2,] 2 4 6


> nonSquareMatrix%^%2#R's %*% will throw the error for us
> nonSquareMatrix %^% 2#R's %*% will throw the error for us
Error in mat %*% (mat %^% (n - 1)) : non-conformable arguments</pre>
Error in mat %*% (mat %^% (n - 1)) : non-conformable arguments</pre>
Our code is far from efficient and could do with more error-checking, but it demonstrates the principle. If we wanted to do this properly, we'd use a library - ideally one that calls C code. Following the previous submission's example, we can just do this:
Our code is far from efficient and could do with more error-checking, but it demonstrates the principle. If we wanted to do this properly, we'd use a library - ideally one that calls C code. Following the previous submission's example, we can just do this:
<lang rsplus>library(Biodem)
<lang rsplus>library(Biodem)
`%^%`<-function(mat,n) Biodem::mtx.exp(mat,n)</lang>
`%^%` <- function(mat, n) Biodem::mtx.exp(mat, n)</lang>
And it will work just the same, except for being much faster and throwing an error on nonSquareMatrix%^%1.
And it will work just the same, except for being much faster and throwing an error on nonSquareMatrix %^% 1.


=={{header|Racket}}==
=={{header|Racket}}==
Retrieved from "https://rosettacode.org/wiki/Matrix-exponentiation_operator"