Matrix-exponentiation operator: Difference between revisions
Content added Content deleted
ReeceGoding (talk | contribs) m (→Infix operator) |
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
||
Line 11: | Line 11: | ||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang="11l">F matrix_mul(m1, m2) |
||
assert(m1[0].len == m2.len) |
assert(m1[0].len == m2.len) |
||
V r = [[0] * m2[0].len] * m1.len |
V r = [[0] * m2[0].len] * m1.len |
||
Line 43: | Line 43: | ||
print("\n10:") |
print("\n10:") |
||
printtable(matrixExp(m, 10))</ |
printtable(matrixExp(m, 10))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 75: | Line 75: | ||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
This is a generic solution for any natural power exponent. It will work with any type that has +,*, additive and multiplicative 0s. The implementation factors out powers A<sup>2<sup>n</sup></sup>: |
This is a generic solution for any natural power exponent. It will work with any type that has +,*, additive and multiplicative 0s. The implementation factors out powers A<sup>2<sup>n</sup></sup>: |
||
< |
<syntaxhighlight lang="ada">with Ada.Text_IO; use Ada.Text_IO; |
||
procedure Test_Matrix is |
procedure Test_Matrix is |
||
Line 164: | Line 164: | ||
Put_Line ("M**10 ="); Put (M**10); |
Put_Line ("M**10 ="); Put (M**10); |
||
Put_Line ("M*M*M*M*M*M*M*M*M*M ="); Put (M*M*M*M*M*M*M*M*M*M); |
Put_Line ("M*M*M*M*M*M*M*M*M*M ="); Put (M*M*M*M*M*M*M*M*M*M); |
||
end Test_Matrix;</ |
end Test_Matrix;</syntaxhighlight> |
||
Sample output: |
Sample output: |
||
<pre> |
<pre> |
||
Line 202: | Line 202: | ||
</pre> |
</pre> |
||
The following program implements exponentiation of a square Hermitian complex matrix by any complex power. The limitation to be Hermitian is not essential and comes for the limitation of the standard Ada linear algebra library. |
The following program implements exponentiation of a square Hermitian complex matrix by any complex power. The limitation to be Hermitian is not essential and comes for the limitation of the standard Ada linear algebra library. |
||
< |
<syntaxhighlight lang="ada">with Ada.Text_IO; use Ada.Text_IO; |
||
with Ada.Complex_Text_IO; use Ada.Complex_Text_IO; |
with Ada.Complex_Text_IO; use Ada.Complex_Text_IO; |
||
with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types; |
with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types; |
||
Line 249: | Line 249: | ||
Put_Line ("M**1 ="); Put (M**(1.0,0.0)); |
Put_Line ("M**1 ="); Put (M**(1.0,0.0)); |
||
Put_Line ("M**0.5 ="); Put (M**(0.5,0.0)); |
Put_Line ("M**0.5 ="); Put (M**(0.5,0.0)); |
||
end Test_Matrix;</ |
end Test_Matrix;</syntaxhighlight> |
||
This solution is not tested, because the available version of GNAT GPL Ada compiler (20070405-41) does not provide an implementation of the standard library. |
This solution is not tested, because the available version of GNAT GPL Ada compiler (20070405-41) does not provide an implementation of the standard library. |
||
Line 257: | Line 257: | ||
{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of '''format'''[ted] ''transput''.}} |
{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of '''format'''[ted] ''transput''.}} |
||
'''File: Matrix_algebra.a68''' |
'''File: Matrix_algebra.a68''' |
||
< |
<syntaxhighlight lang="algol68">INT default upb=3; |
||
MODE VEC = [default upb]COSCAL; |
MODE VEC = [default upb]COSCAL; |
||
MODE MAT = [default upb,default upb]COSCAL; |
MODE MAT = [default upb,default upb]COSCAL; |
||
Line 287: | Line 287: | ||
OD; |
OD; |
||
out |
out |
||
);</ |
);</syntaxhighlight>'''File: Matrix-exponentiation_operator.a68''' |
||
< |
<syntaxhighlight lang="algol68">OP ** = (MAT base, INT exponent)MAT: ( |
||
BITS binary exponent:=BIN exponent ; |
BITS binary exponent:=BIN exponent ; |
||
MAT out := IF bits width ELEM binary exponent THEN base ELSE IDENTITY UPB base FI; |
MAT out := IF bits width ELEM binary exponent THEN base ELSE IDENTITY UPB base FI; |
||
Line 301: | Line 301: | ||
OD; |
OD; |
||
out |
out |
||
);</ |
);</syntaxhighlight>'''File: test_Matrix-exponentiation_operator.a68''' |
||
< |
<syntaxhighlight lang="algol68">#!/usr/local/bin/a68g --script # |
||
MODE COSCAL = COMPL; |
MODE COSCAL = COMPL; |
||
Line 325: | Line 325: | ||
printf(($" mat ** "g(0)":"l$,24)); |
printf(($" mat ** "g(0)":"l$,24)); |
||
compl mat printf(scal fmt, mat**24); |
compl mat printf(scal fmt, mat**24); |
||
print(newline)</ |
print(newline)</syntaxhighlight> |
||
Output: |
Output: |
||
<pre> |
<pre> |
||
Line 335: | Line 335: | ||
=={{header|BBC BASIC}}== |
=={{header|BBC BASIC}}== |
||
< |
<syntaxhighlight lang="bbcbasic"> DIM matrix(1,1), output(1,1) |
||
matrix() = 3, 2, 2, 1 |
matrix() = 3, 2, 2, 1 |
||
Line 359: | Line 359: | ||
NEXT |
NEXT |
||
ENDIF |
ENDIF |
||
ENDPROC</ |
ENDPROC</syntaxhighlight> |
||
Output: |
Output: |
||
<pre> |
<pre> |
||
Line 397: | Line 397: | ||
Matrix multiplication is a known idiom taken from BQN crate. Matrix exponentiation is simply doing Matrix multiplication n times. |
Matrix multiplication is a known idiom taken from BQN crate. Matrix exponentiation is simply doing Matrix multiplication n times. |
||
< |
<syntaxhighlight lang="bqn">MatMul ← +˝∘×⎉1‿∞ |
||
MatEx ← {𝕨 MatMul⍟(𝕩-1) 𝕨} |
MatEx ← {𝕨 MatMul⍟(𝕩-1) 𝕨} |
||
(>⟨3‿2 |
(>⟨3‿2 |
||
2‿1⟩) MatEx 1‿2‿3‿4‿10</ |
2‿1⟩) MatEx 1‿2‿3‿4‿10</syntaxhighlight><syntaxhighlight lang="bqn">┌─ |
||
· ┌─ ┌─ ┌─ ┌─ ┌─ |
· ┌─ ┌─ ┌─ ┌─ ┌─ |
||
╵ 3 2 ╵ 13 8 ╵ 55 34 ╵ 233 144 ╵ 1346269 832040 |
╵ 3 2 ╵ 13 8 ╵ 55 34 ╵ 233 144 ╵ 1346269 832040 |
||
2 1 8 5 34 21 144 89 832040 514229 |
2 1 8 5 34 21 144 89 832040 514229 |
||
┘ ┘ ┘ ┘ ┘ |
┘ ┘ ┘ ┘ ┘ |
||
┘</ |
┘</syntaxhighlight> |
||
For larger exponents it's more efficient to use a fast exponentiation pattern that builds large powers quickly with repeated squaring, then multiplies the appropriate power-of-two exponents together. |
For larger exponents it's more efficient to use a fast exponentiation pattern that builds large powers quickly with repeated squaring, then multiplies the appropriate power-of-two exponents together. |
||
< |
<syntaxhighlight lang="bqn">MatEx ← MatMul{𝔽´𝔽˜⍟(/2|⌊∘÷⟜2⍟(↕1+·⌊2⋆⁼⊢)𝕩)𝕨}</syntaxhighlight> |
||
=={{header|Burlesque}}== |
=={{header|Burlesque}}== |
||
< |
<syntaxhighlight lang="burlesque">blsq ) {{1 1} {1 0}} 10 .*{mm}r[ |
||
{{89 55} {55 34}}</ |
{{89 55} {55 34}}</syntaxhighlight> |
||
=={{header|C}}== |
=={{header|C}}== |
||
C doesn't support classes or allow operator overloading. The following is code that defines a function, <tt>SquareMtxPower</tt> that will raise a matrix to a positive integer power. |
C doesn't support classes or allow operator overloading. The following is code that defines a function, <tt>SquareMtxPower</tt> that will raise a matrix to a positive integer power. |
||
< |
<syntaxhighlight lang="c">#include <math.h> |
||
#include <stdio.h> |
#include <stdio.h> |
||
#include <stdlib.h> |
#include <stdlib.h> |
||
Line 602: | Line 602: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
Output: |
Output: |
||
<pre>m0 dim:3 = |
<pre>m0 dim:3 = |
||
Line 626: | Line 626: | ||
=={{header|C sharp}}== |
=={{header|C sharp}}== |
||
< |
<syntaxhighlight lang="csharp">using System; |
||
using System.Collections; |
using System.Collections; |
||
using System.Collections.Generic; |
using System.Collections.Generic; |
||
Line 684: | Line 684: | ||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre style="height:30ex;overflow:scroll"> |
<pre style="height:30ex;overflow:scroll"> |
||
Line 707: | Line 707: | ||
=={{header|C++}}== |
=={{header|C++}}== |
||
This is an implementation in C++. |
This is an implementation in C++. |
||
< |
<syntaxhighlight lang="cpp">#include <complex> |
||
#include <cmath> |
#include <cmath> |
||
#include <iostream> |
#include <iostream> |
||
Line 755: | Line 755: | ||
} |
} |
||
return d; |
return d; |
||
}</ |
}</syntaxhighlight> |
||
This is the task part. |
This is the task part. |
||
< |
<syntaxhighlight lang="cpp"> // C++ does not have a ** operator, instead, ^ (bitwise Xor) is used. |
||
Mx operator^(int n) { |
Mx operator^(int n) { |
||
if (n < 0) |
if (n < 0) |
||
Line 784: | Line 784: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 798: | Line 798: | ||
This uses the '*' operator for arrays as defined in [[Matrix_multiplication#Chapel]] |
This uses the '*' operator for arrays as defined in [[Matrix_multiplication#Chapel]] |
||
< |
<syntaxhighlight lang="chapel">proc **(a, e) { |
||
// create result matrix of same dimensions |
// create result matrix of same dimensions |
||
var r:[a.domain] a.eltType; |
var r:[a.domain] a.eltType; |
||
Line 809: | Line 809: | ||
return r; |
return r; |
||
}</ |
}</syntaxhighlight> |
||
Usage example (like Perl): |
Usage example (like Perl): |
||
< |
<syntaxhighlight lang="chapel">var m:[1..3, 1..3] int; |
||
m(1,1) = 1; m(1,2) = 2; m(1,3) = 0; |
m(1,1) = 1; m(1,2) = 2; m(1,3) = 0; |
||
m(2,1) = 0; m(2,2) = 3; m(2,3) = 1; |
m(2,1) = 0; m(2,2) = 3; m(2,3) = 1; |
||
Line 822: | Line 822: | ||
writeln("Order ", i); |
writeln("Order ", i); |
||
writeln(m ** i, "\n"); |
writeln(m ** i, "\n"); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 882: | Line 882: | ||
=={{header|Common Lisp}}== |
=={{header|Common Lisp}}== |
||
This Common Lisp implementation uses 2D Arrays to represent matrices, and checks to make sure that the arrays are the right dimensions for multiplication and square for exponentiation. |
This Common Lisp implementation uses 2D Arrays to represent matrices, and checks to make sure that the arrays are the right dimensions for multiplication and square for exponentiation. |
||
< |
<syntaxhighlight lang="lisp">(defun multiply-matrices (matrix-0 matrix-1) |
||
"Takes two 2D arrays and returns their product, or an error if they cannot be multiplied" |
"Takes two 2D arrays and returns their product, or an error if they cannot be multiplied" |
||
(let* ((m0-dims (array-dimensions matrix-0)) |
(let* ((m0-dims (array-dimensions matrix-0)) |
||
Line 940: | Line 940: | ||
(multiply-matrices me2 me2))) |
(multiply-matrices me2 me2))) |
||
(t (let ((me2 (matrix-expt matrix (/ (1- exp) 2)))) |
(t (let ((me2 (matrix-expt matrix (/ (1- exp) 2)))) |
||
(multiply-matrices matrix (multiply-matrices me2 me2)))))))</ |
(multiply-matrices matrix (multiply-matrices me2 me2)))))))</syntaxhighlight> |
||
Output (note that this lisp implementation uses single-precision floats for decimals by default). We can also use rationals: |
Output (note that this lisp implementation uses single-precision floats for decimals by default). We can also use rationals: |
||
CL-USER> (setf 5x5-matrix |
CL-USER> (setf 5x5-matrix |
||
Line 976: | Line 976: | ||
=={{header|D}}== |
=={{header|D}}== |
||
< |
<syntaxhighlight lang="d">import std.stdio, std.string, std.math, std.array, std.algorithm; |
||
struct SquareMat(T = creal) { |
struct SquareMat(T = creal) { |
||
Line 1,050: | Line 1,050: | ||
foreach (immutable p; [0, 1, 23, 24]) |
foreach (immutable p; [0, 1, 23, 24]) |
||
writefln("m ^^ %d =\n%s", p, m ^^ p); |
writefln("m ^^ %d =\n%s", p, m ^^ p); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>m ^^ 0 = |
<pre>m ^^ 0 = |
||
Line 1,070: | Line 1,070: | ||
=={{header|Delphi}}== |
=={{header|Delphi}}== |
||
<syntaxhighlight lang="delphi"> |
|||
<lang Delphi> |
|||
program Matrix_exponentiation_operator; |
program Matrix_exponentiation_operator; |
||
Line 1,212: | Line 1,212: | ||
Readln; |
Readln; |
||
end. |
end. |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,241: | Line 1,241: | ||
| 2 1 | |
| 2 1 | |
||
</pre> |
</pre> |
||
< |
<syntaxhighlight lang="erre">PROGRAM MAT_PROD |
||
!$MATRIX |
!$MATRIX |
||
Line 1,293: | Line 1,293: | ||
END FOR |
END FOR |
||
END PROGRAM</ |
END PROGRAM</syntaxhighlight> |
||
Sample output: |
Sample output: |
||
<pre> |
<pre> |
||
Line 1,303: | Line 1,303: | ||
There is already a built-in word (<code>m^n</code>) that implements exponentiation. Here is a simple and less efficient implementation. |
There is already a built-in word (<code>m^n</code>) that implements exponentiation. Here is a simple and less efficient implementation. |
||
< |
<syntaxhighlight lang="factor">USING: kernel math math.matrices sequences ; |
||
: my-m^n ( m n -- m' ) |
: my-m^n ( m n -- m' ) |
||
Line 1,309: | Line 1,309: | ||
[ drop length identity-matrix ] |
[ drop length identity-matrix ] |
||
[ swap '[ _ m. ] times ] 2bi |
[ swap '[ _ m. ] times ] 2bi |
||
] if ;</ |
] if ;</syntaxhighlight> |
||
( scratchpad ) { { 3 2 } { 2 1 } } 0 my-m^n . |
( scratchpad ) { { 3 2 } { 2 1 } } 0 my-m^n . |
||
Line 1,318: | Line 1,318: | ||
=={{header|Fermat}}== |
=={{header|Fermat}}== |
||
Matrix exponentiation for square matrices and integer powers is built in. |
Matrix exponentiation for square matrices and integer powers is built in. |
||
< |
<syntaxhighlight lang="fermat"> |
||
Array a[2,2]; {illustrate with a 2x2 matrix} |
Array a[2,2]; {illustrate with a 2x2 matrix} |
||
[a]:=[(2/3, 1/3, 4/5, 1/5)]; |
[a]:=[(2/3, 1/3, 4/5, 1/5)]; |
||
Line 1,326: | Line 1,326: | ||
[a]^3; |
[a]^3; |
||
[a]^10; |
[a]^10; |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,346: | Line 1,346: | ||
=={{header|Fortran}}== |
=={{header|Fortran}}== |
||
{{works with|Fortran|90 and later}} |
{{works with|Fortran|90 and later}} |
||
< |
<syntaxhighlight lang="fortran">module matmod |
||
implicit none |
implicit none |
||
Line 1,398: | Line 1,398: | ||
end do |
end do |
||
end program Matrix_exponentiation</ |
end program Matrix_exponentiation</syntaxhighlight> |
||
Output |
Output |
||
<pre> 1.00000 0.00000 0.00000 |
<pre> 1.00000 0.00000 0.00000 |
||
Line 1,425: | Line 1,425: | ||
This operator performs M^n for any square invertible matrix M and integer n, including negative powers. |
This operator performs M^n for any square invertible matrix M and integer n, including negative powers. |
||
< |
<syntaxhighlight lang="freebasic">#include once "matmult.bas" |
||
#include once "rowech.bas" |
#include once "rowech.bas" |
||
#include once "matinv.bas" |
#include once "matinv.bas" |
||
Line 1,450: | Line 1,450: | ||
next i |
next i |
||
print |
print |
||
next n</ |
next n</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> 308.9999999999998 -307.9999999999998 |
<pre> 308.9999999999998 -307.9999999999998 |
||
Line 1,474: | Line 1,474: | ||
=={{header|GAP}}== |
=={{header|GAP}}== |
||
< |
<syntaxhighlight lang="gap"># Matrix exponentiation is built-in |
||
A := [[0 , 1], [1, 1]]; |
A := [[0 , 1], [1, 1]]; |
||
PrintArray(A); |
PrintArray(A); |
||
Line 1,481: | Line 1,481: | ||
PrintArray(A^10); |
PrintArray(A^10); |
||
# [ [ 34, 55 ], |
# [ [ 34, 55 ], |
||
# [ 55, 89 ] ]</ |
# [ 55, 89 ] ]</syntaxhighlight> |
||
=={{header|Go}}== |
=={{header|Go}}== |
||
Line 1,487: | Line 1,487: | ||
<br> |
<br> |
||
Like some other languages here, Go doesn't have a symbolic operator for numeric exponentiation and even if it did doesn't support operator overloading. We therefore write the exponentiation operation for matrices as an equivalent 'pow' function. |
Like some other languages here, Go doesn't have a symbolic operator for numeric exponentiation and even if it did doesn't support operator overloading. We therefore write the exponentiation operation for matrices as an equivalent 'pow' function. |
||
< |
<syntaxhighlight lang="go">package main |
||
import "fmt" |
import "fmt" |
||
Line 1,557: | Line 1,557: | ||
fmt.Println() |
fmt.Println() |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,599: | Line 1,599: | ||
Instead of writing it directly, we can re-use the built-in [[exponentiation operator]] if we declare matrices as an instance of ''Num'', using [[matrix multiplication]] (and addition). For simplicity, we use the inefficient representation as list of lists. Note that we don't check the dimensions (there are several ways to do that on the type-level, for example with phantom types). |
Instead of writing it directly, we can re-use the built-in [[exponentiation operator]] if we declare matrices as an instance of ''Num'', using [[matrix multiplication]] (and addition). For simplicity, we use the inefficient representation as list of lists. Note that we don't check the dimensions (there are several ways to do that on the type-level, for example with phantom types). |
||
< |
<syntaxhighlight lang="haskell">import Data.List (transpose) |
||
(<+>) |
(<+>) |
||
Line 1,630: | Line 1,630: | ||
-- TEST ---------------------------------------------------------------------- |
-- TEST ---------------------------------------------------------------------- |
||
main :: IO () |
main :: IO () |
||
main = print $ Mat [[1, 2], [0, 1]] ^ 4</ |
main = print $ Mat [[1, 2], [0, 1]] ^ 4</syntaxhighlight> |
||
{{Out}} |
{{Out}} |
||
<pre>Mat [[1,8],[0,1]]</pre> |
<pre>Mat [[1,8],[0,1]]</pre> |
||
Line 1,640: | Line 1,640: | ||
===With Numeric.LinearAlgebra=== |
===With Numeric.LinearAlgebra=== |
||
< |
<syntaxhighlight lang="haskell">import Numeric.LinearAlgebra |
||
a :: Matrix I |
a :: Matrix I |
||
Line 1,650: | Line 1,650: | ||
print $ a^4 |
print $ a^4 |
||
putStrLn "power of zero: " |
putStrLn "power of zero: " |
||
print $ a^0</ |
print $ a^0</syntaxhighlight> |
||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
Line 1,663: | Line 1,663: | ||
< |
<syntaxhighlight lang="j">mp=: +/ .* NB. Matrix multiplication |
||
pow=: pow0=: 4 : 'mp&x^:y =i.#x'</ |
pow=: pow0=: 4 : 'mp&x^:y =i.#x'</syntaxhighlight> |
||
or, from [[j:Essays/Linear Recurrences|the J wiki]], and faster for large exponents: |
or, from [[j:Essays/Linear Recurrences|the J wiki]], and faster for large exponents: |
||
< |
<syntaxhighlight lang="j">pow=: pow1=: 4 : 'mp/ mp~^:(I.|.#:y) x'</syntaxhighlight> |
||
This implements an optimization where the exponent is represented in base 2, and repeated squaring is used to create a list of relevant powers of the base matrix, which are then combined using matrix multiplication. Note, however, that these two definitions treat a zero exponent differently (m pow0 0 gives an identity matrix whose shape matches m, while m pow1 0 gives a scalar 1). |
This implements an optimization where the exponent is represented in base 2, and repeated squaring is used to create a list of relevant powers of the base matrix, which are then combined using matrix multiplication. Note, however, that these two definitions treat a zero exponent differently (m pow0 0 gives an identity matrix whose shape matches m, while m pow1 0 gives a scalar 1). |
||
Line 1,682: | Line 1,682: | ||
Extends [[Matrix Transpose#JavaScript]] and [[Matrix multiplication#JavaScript]] |
Extends [[Matrix Transpose#JavaScript]] and [[Matrix multiplication#JavaScript]] |
||
< |
<syntaxhighlight lang="javascript">// IdentityMatrix is a "subclass" of Matrix |
||
function IdentityMatrix(n) { |
function IdentityMatrix(n) { |
||
this.height = n; |
this.height = n; |
||
Line 1,707: | Line 1,707: | ||
var m = new Matrix([[3, 2], [2, 1]]); |
var m = new Matrix([[3, 2], [2, 1]]); |
||
[0,1,2,3,4,10].forEach(function(e){print(m.exp(e)); print()})</ |
[0,1,2,3,4,10].forEach(function(e){print(m.exp(e)); print()})</syntaxhighlight> |
||
output |
output |
||
<pre>1,0 |
<pre>1,0 |
||
Line 1,735: | Line 1,735: | ||
matrix_exp(n) adopts a "divide-and-conquer" strategy to avoid unnecessarily many matrix multiplications. The implementation uses direct_matrix_exp(n) for small n; this function could be defined as an inner function, but is defined separately first for clarity, and second to simplify timing comparisons, as shown below. |
matrix_exp(n) adopts a "divide-and-conquer" strategy to avoid unnecessarily many matrix multiplications. The implementation uses direct_matrix_exp(n) for small n; this function could be defined as an inner function, but is defined separately first for clarity, and second to simplify timing comparisons, as shown below. |
||
< |
<syntaxhighlight lang="jq"># produce an array of length n that is 1 at i and 0 elsewhere |
||
def indicator(i;n): [range(0;n) | 0] | .[i] = 1; |
def indicator(i;n): [range(0;n) | 0] | .[i] = 1; |
||
Line 1,758: | Line 1,758: | ||
| multiply($ans; $residue ) |
| multiply($ans; $residue ) |
||
end |
end |
||
end;</ |
end;</syntaxhighlight> |
||
'''Examples''' |
'''Examples''' |
||
The execution speeds of matrix_exp and direct_matrix_exp are compared using a one-eighth-rotation matrix, which |
The execution speeds of matrix_exp and direct_matrix_exp are compared using a one-eighth-rotation matrix, which |
||
is raised to the 10,000th power. The direct method turns out to be almost as fast. |
is raised to the 10,000th power. The direct method turns out to be almost as fast. |
||
< |
<syntaxhighlight lang="jq">def pi: 4 * (1|atan); |
||
def rotation_matrix(theta): |
def rotation_matrix(theta): |
||
Line 1,771: | Line 1,771: | ||
def demo_direct_matrix_exp(n): |
def demo_direct_matrix_exp(n): |
||
rotation_matrix( pi / 4 ) | direct_matrix_exp(n) ;</ |
rotation_matrix( pi / 4 ) | direct_matrix_exp(n) ;</syntaxhighlight> |
||
'''Results''': |
'''Results''': |
||
< |
<syntaxhighlight lang="sh"># For demo_matrix_exp(10000) |
||
$ time jq -n -c -f Matrix-exponentiation_operator.rc |
$ time jq -n -c -f Matrix-exponentiation_operator.rc |
||
[[1,-1.1102230246251565e-12],[1.1102230246251565e-12,1]] |
[[1,-1.1102230246251565e-12],[1.1102230246251565e-12,1]] |
||
user 0m0.490s |
user 0m0.490s |
||
sys 0m0.008s</ |
sys 0m0.008s</syntaxhighlight> |
||
< |
<syntaxhighlight lang="sh"># For demo_direct_matrix_exp(10000) |
||
$ time jq -n -c -f Matrix-exponentiation_operator.rc |
$ time jq -n -c -f Matrix-exponentiation_operator.rc |
||
[[1,-7.849831895612169e-13],[7.849831895612169e-13,1]] |
[[1,-7.849831895612169e-13],[7.849831895612169e-13,1]] |
||
user 0m0.625s |
user 0m0.625s |
||
sys 0m0.006s</ |
sys 0m0.006s</syntaxhighlight> |
||
=={{header|Jsish}}== |
=={{header|Jsish}}== |
||
Line 1,789: | Line 1,789: | ||
Uses module listed in [[Matrix Transpose#Jsish]]. Fails the task spec actually, as Matrix.exp() is implemented as a method, not an operator. |
Uses module listed in [[Matrix Transpose#Jsish]]. Fails the task spec actually, as Matrix.exp() is implemented as a method, not an operator. |
||
< |
<syntaxhighlight lang="javascript">/* Matrix exponentiation, in Jsish */ |
||
require('Matrix'); |
require('Matrix'); |
||
Line 1,811: | Line 1,811: | ||
m.exp(10) ==> { height:2, mtx:[ [ 1346269, 832040 ], [ 832040, 514229 ] ], width:2 } |
m.exp(10) ==> { height:2, mtx:[ [ 1346269, 832040 ], [ 832040, 514229 ] ], width:2 } |
||
=!EXPECTEND!= |
=!EXPECTEND!= |
||
*/</ |
*/</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,819: | Line 1,819: | ||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
Matrix exponentiation is implemented by the built-in <code>^</code> operator. |
Matrix exponentiation is implemented by the built-in <code>^</code> operator. |
||
< |
<syntaxhighlight lang="julia">julia> [1 1 ; 1 0]^10 |
||
2x2 Array{Int64,2}: |
2x2 Array{Int64,2}: |
||
89 55 |
89 55 |
||
55 34</ |
55 34</syntaxhighlight> |
||
=={{header|K}}== |
=={{header|K}}== |
||
<syntaxhighlight lang="k"> |
|||
<lang K> |
|||
/Matrix Exponentiation |
/Matrix Exponentiation |
||
/mpow.k |
/mpow.k |
||
pow: {:[0=y; :({a=/:a:!x}(#x))];a: x; do[y-1; a: x _mul a]; :a} |
pow: {:[0=y; :({a=/:a:!x}(#x))];a: x; do[y-1; a: x _mul a]; :a} |
||
</syntaxhighlight> |
|||
</lang> |
|||
The output of a session is given below: |
The output of a session is given below: |
||
{{out}} |
{{out}} |
||
Line 1,863: | Line 1,863: | ||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
< |
<syntaxhighlight lang="scala">// version 1.1.3 |
||
typealias Vector = DoubleArray |
typealias Vector = DoubleArray |
||
Line 1,919: | Line 1,919: | ||
) |
) |
||
for (i in 0..10) printMatrix(m pow i, i) |
for (i in 0..10) printMatrix(m pow i, i) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,969: | Line 1,969: | ||
=={{header|Lambdatalk}}== |
=={{header|Lambdatalk}}== |
||
< |
<syntaxhighlight lang="scheme"> |
||
{require lib_matrix} |
{require lib_matrix} |
||
Line 1,993: | Line 1,993: | ||
M^4 = [[233,144],[144,89]] |
M^4 = [[233,144],[144,89]] |
||
M^10 = [[1346269,832040],[832040,514229]] |
M^10 = [[1346269,832040],[832040,514229]] |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Liberty BASIC}}== |
=={{header|Liberty BASIC}}== |
||
There is no native matrix capability. A set of functions is available at http://www.diga.me.uk/RCMatrixFuncs.bas implementing matrices of arbitrary dimension in a string format. |
There is no native matrix capability. A set of functions is available at http://www.diga.me.uk/RCMatrixFuncs.bas implementing matrices of arbitrary dimension in a string format. |
||
<syntaxhighlight lang="lb"> |
|||
<lang lb> |
|||
MatrixD$ ="3, 3, 0.86603, 0.50000, 0.00000, -0.50000, 0.86603, 0.00000, 0.00000, 0.00000, 1.00000" |
MatrixD$ ="3, 3, 0.86603, 0.50000, 0.00000, -0.50000, 0.86603, 0.00000, 0.00000, 0.00000, 1.00000" |
||
Line 2,009: | Line 2,009: | ||
MatrixE$ =MatrixToPower$( MatrixD$, 9) |
MatrixE$ =MatrixToPower$( MatrixD$, 9) |
||
call DisplayMatrix MatrixE$ |
call DisplayMatrix MatrixE$ |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
Line 2,028: | Line 2,028: | ||
=={{header|Lua}}== |
=={{header|Lua}}== |
||
< |
<syntaxhighlight lang="lua">Matrix = {} |
||
function Matrix.new( dim_y, dim_x ) |
function Matrix.new( dim_y, dim_x ) |
||
Line 2,123: | Line 2,123: | ||
n = m^4; |
n = m^4; |
||
Matrix.Show( n )</ |
Matrix.Show( n )</syntaxhighlight> |
||
=={{header|M2000 Interpreter}}== |
=={{header|M2000 Interpreter}}== |
||
<syntaxhighlight lang="m2000 interpreter"> |
|||
<lang M2000 Interpreter> |
|||
Module CheckIt { |
Module CheckIt { |
||
Class cArray { |
Class cArray { |
||
Line 2,188: | Line 2,188: | ||
} |
} |
||
Checkit |
Checkit |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
Line 2,229: | Line 2,229: | ||
=={{header|Maple}}== |
=={{header|Maple}}== |
||
Maple handles matrix powers implicitly with the built-in exponentiation operator: |
Maple handles matrix powers implicitly with the built-in exponentiation operator: |
||
< |
<syntaxhighlight lang="maple">> M := <<1,2>|<3,4>>; |
||
> M ^ 2;</ |
> M ^ 2;</syntaxhighlight> |
||
<math>\left[\begin{array}{cc} |
<math>\left[\begin{array}{cc} |
||
7 & 15 \\ |
7 & 15 \\ |
||
Line 2,237: | Line 2,237: | ||
If you want elementwise powers, you can use the elementwise <code>^~</code> operator: |
If you want elementwise powers, you can use the elementwise <code>^~</code> operator: |
||
< |
<syntaxhighlight lang="maple">> M := <<1,2>|<3,4>>; |
||
> M ^~ 2;</ |
> M ^~ 2;</syntaxhighlight> |
||
<math>\left[\begin{array}{cc} |
<math>\left[\begin{array}{cc} |
||
1 & 9 \\ |
1 & 9 \\ |
||
Line 2,246: | Line 2,246: | ||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
||
In Mathematica there is an distinction between powering elements wise and as a matrix. So m^2 will give m with each element squared. To do matrix exponentation we use the function MatrixPower. It can handle all types of numbers for the power (integers, floats, rationals, complex) but also symbols for the power, and all types for the matrix (numbers, symbols et cetera), and will always keep the result exact if the matrix and the exponent is exact. |
In Mathematica there is an distinction between powering elements wise and as a matrix. So m^2 will give m with each element squared. To do matrix exponentation we use the function MatrixPower. It can handle all types of numbers for the power (integers, floats, rationals, complex) but also symbols for the power, and all types for the matrix (numbers, symbols et cetera), and will always keep the result exact if the matrix and the exponent is exact. |
||
< |
<syntaxhighlight lang="mathematica">a = {{3, 2}, {4, 1}}; |
||
MatrixPower[a, 0] |
MatrixPower[a, 0] |
||
MatrixPower[a, 1] |
MatrixPower[a, 1] |
||
Line 2,252: | Line 2,252: | ||
MatrixPower[a, 4] |
MatrixPower[a, 4] |
||
MatrixPower[a, 1/2] |
MatrixPower[a, 1/2] |
||
MatrixPower[a, Pi]</ |
MatrixPower[a, Pi]</syntaxhighlight> |
||
gives back: |
gives back: |
||
Line 2,310: | Line 2,310: | ||
Symbolic matrices like {{i,j},{k,l}} to the power m give general solutions for all possible i,j,k,l, and m: |
Symbolic matrices like {{i,j},{k,l}} to the power m give general solutions for all possible i,j,k,l, and m: |
||
< |
<syntaxhighlight lang="mathematica">MatrixPower[{{i, j}, {k, l}}, m] // Simplify</syntaxhighlight> |
||
gives back (note that the simplification is not necessary for the evaluation, it just gives a shorter output): |
gives back (note that the simplification is not necessary for the evaluation, it just gives a shorter output): |
||
Line 2,332: | Line 2,332: | ||
=={{header|MATLAB}}== |
=={{header|MATLAB}}== |
||
For exponents in the form of A*A*A*A*...*A, A must be a square matrix: |
For exponents in the form of A*A*A*A*...*A, A must be a square matrix: |
||
< |
<syntaxhighlight lang="matlab">function [output] = matrixexponentiation(matrixA, exponent) |
||
output = matrixA^(exponent);</ |
output = matrixA^(exponent);</syntaxhighlight> |
||
Otherwise, to take the individual array elements to the power of an exponent (the matrix need not be square): |
Otherwise, to take the individual array elements to the power of an exponent (the matrix need not be square): |
||
< |
<syntaxhighlight lang="matlab">function [output] = matrixexponentiation(matrixA, exponent) |
||
output = matrixA.^(exponent);</ |
output = matrixA.^(exponent);</syntaxhighlight> |
||
=={{header|Maxima}}== |
=={{header|Maxima}}== |
||
< |
<syntaxhighlight lang="maxima">a: matrix([3, 2], |
||
[4, 1])$ |
[4, 1])$ |
||
Line 2,349: | Line 2,349: | ||
a ^^ -1; |
a ^^ -1; |
||
/* matrix([-1/5, 2/5], |
/* matrix([-1/5, 2/5], |
||
[4/5, -3/5]) */</ |
[4/5, -3/5]) */</syntaxhighlight> |
||
=={{header|Nim}}== |
=={{header|Nim}}== |
||
< |
<syntaxhighlight lang="nim">import sequtils, strutils |
||
type Matrix[N: static int; T] = array[1..N, array[1..N, T]] |
type Matrix[N: static int; T] = array[1..N, array[1..N, T]] |
||
Line 2,402: | Line 2,402: | ||
S30 = 1 / 2 |
S30 = 1 / 2 |
||
let m2: Matrix[2, float] = [[C30, -S30], [S30, C30]] # 30° rotation matrix. |
let m2: Matrix[2, float] = [[C30, -S30], [S30, C30]] # 30° rotation matrix. |
||
echo m2^12 # Nearly the identity matrix.</ |
echo m2^12 # Nearly the identity matrix.</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 2,416: | Line 2,416: | ||
We will use some auxiliary functions |
We will use some auxiliary functions |
||
< |
<syntaxhighlight lang="ocaml">(* identity matrix *) |
||
let eye n = |
let eye n = |
||
let a = Array.make_matrix n n 0.0 in |
let a = Array.make_matrix n n 0.0 in |
||
Line 2,473: | Line 2,473: | ||
(* example with integers *) |
(* example with integers *) |
||
pow 1 ( * ) 2 16;; |
pow 1 ( * ) 2 16;; |
||
(* - : int = 65536 *)</ |
(* - : int = 65536 *)</syntaxhighlight> |
||
Now matrix power is simply a special case of pow : |
Now matrix power is simply a special case of pow : |
||
< |
<syntaxhighlight lang="ocaml">let matpow a n = |
||
let p, q = dim a in |
let p, q = dim a in |
||
if p <> q then failwith "bad dimensions" else |
if p <> q then failwith "bad dimensions" else |
||
Line 2,489: | Line 2,489: | ||
[| [| 1.0; 1.0|]; [| 1.0; 0.0 |] |] ^^ 10;; |
[| [| 1.0; 1.0|]; [| 1.0; 0.0 |] |] ^^ 10;; |
||
(* - : float array array = [|[|89.; 55.|]; [|55.; 34.|]|] *)</ |
(* - : float array array = [|[|89.; 55.|]; [|55.; 34.|]|] *)</syntaxhighlight> |
||
=={{header|Octave}}== |
=={{header|Octave}}== |
||
Line 2,495: | Line 2,495: | ||
Of course GNU Octave handles matrix and operations on matrix "naturally". |
Of course GNU Octave handles matrix and operations on matrix "naturally". |
||
< |
<syntaxhighlight lang="octave">M = [ 3, 2; 2, 1 ]; |
||
M^0 |
M^0 |
||
M^1 |
M^1 |
||
M^2 |
M^2 |
||
M^(-1) |
M^(-1) |
||
M^0.5</ |
M^0.5</syntaxhighlight> |
||
Output: |
Output: |
||
Line 2,532: | Line 2,532: | ||
=={{header|PARI/GP}}== |
=={{header|PARI/GP}}== |
||
<lang |
<syntaxhighlight lang="parigp">M^n</syntaxhighlight> |
||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
< |
<syntaxhighlight lang="perl">use strict; |
||
package SquareMatrix; |
package SquareMatrix; |
||
use Carp; # standard, "it's not my fault" module |
use Carp; # standard, "it's not my fault" module |
||
Line 2,628: | Line 2,628: | ||
print "\n### WAY too big:\n", $m ** 1_000_000_000_000; |
print "\n### WAY too big:\n", $m ** 1_000_000_000_000; |
||
print "\n### But identity matrix can handle that\n", |
print "\n### But identity matrix can handle that\n", |
||
$m->identity ** 1_000_000_000_000;</ |
$m->identity ** 1_000_000_000_000;</syntaxhighlight> |
||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
Phix does not permit operator overloading, however here is a simple function to raise a square matrix to a non-negative integer power.<br> |
Phix does not permit operator overloading, however here is a simple function to raise a square matrix to a non-negative integer power.<br> |
||
First two routines copied straight from the [[Identity_matrix#Phix|Identity_matrix]] and [[Matrix_multiplication#Phix|Matrix_multiplication]] tasks. |
First two routines copied straight from the [[Identity_matrix#Phix|Identity_matrix]] and [[Matrix_multiplication#Phix|Matrix_multiplication]] tasks. |
||
<!--< |
<!--<syntaxhighlight lang="phix">(phixonline)--> |
||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
||
<span style="color: #008080;">function</span> <span style="color: #000000;">identity</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
<span style="color: #008080;">function</span> <span style="color: #000000;">identity</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
||
Line 2,687: | Line 2,687: | ||
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"==\n"</span><span style="color: #0000FF;">)</span> |
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"==\n"</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">matrix_exponent</span><span style="color: #0000FF;">(</span><span style="color: #000000;">identity</span><span style="color: #0000FF;">(</span><span style="color: #000000;">4</span><span style="color: #0000FF;">),</span><span style="color: #000000;">5</span><span style="color: #0000FF;">))</span> |
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">matrix_exponent</span><span style="color: #0000FF;">(</span><span style="color: #000000;">identity</span><span style="color: #0000FF;">(</span><span style="color: #000000;">4</span><span style="color: #0000FF;">),</span><span style="color: #000000;">5</span><span style="color: #0000FF;">))</span> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,715: | Line 2,715: | ||
=={{header|PicoLisp}}== |
=={{header|PicoLisp}}== |
||
Uses the 'matMul' function from [[Matrix multiplication#PicoLisp]] |
Uses the 'matMul' function from [[Matrix multiplication#PicoLisp]] |
||
< |
<syntaxhighlight lang="picolisp">(de matIdent (N) |
||
(let L (need N (1) 0) |
(let L (need N (1) 0) |
||
(mapcar '(() (copy (rot L))) L) ) ) |
(mapcar '(() (copy (rot L))) L) ) ) |
||
Line 2,725: | Line 2,725: | ||
M ) ) |
M ) ) |
||
(matExp '((3 2) (2 1)) 3)</ |
(matExp '((3 2) (2 1)) 3)</syntaxhighlight> |
||
Output: |
Output: |
||
<pre>-> ((55 34) (34 21))</pre> |
<pre>-> ((55 34) (34 21))</pre> |
||
Line 2,731: | Line 2,731: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
Using matrixMul from [[Matrix multiplication#Python]] |
Using matrixMul from [[Matrix multiplication#Python]] |
||
< |
<syntaxhighlight lang="python">>>> from operator import mul |
||
>>> def matrixMul(m1, m2): |
>>> def matrixMul(m1, m2): |
||
return map( |
return map( |
||
Line 2,786: | Line 2,786: | ||
1346269 832040 |
1346269 832040 |
||
832040 514229 |
832040 514229 |
||
>>></ |
>>></syntaxhighlight> |
||
Alternative Based Upon @ operator of Python 3.5 PEP 465 and using Matrix exponentation for faster computation of powers |
Alternative Based Upon @ operator of Python 3.5 PEP 465 and using Matrix exponentation for faster computation of powers |
||
<lang> |
<syntaxhighlight lang="text"> |
||
class Mat(list) : |
class Mat(list) : |
||
def __matmul__(self, B) : |
def __matmul__(self, B) : |
||
Line 2,822: | Line 2,822: | ||
print('\n%i:' % i) |
print('\n%i:' % i) |
||
printtable(power(m, i)) |
printtable(power(m, i)) |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{Output}} |
{{Output}} |
||
<pre> |
<pre> |
||
Line 2,844: | Line 2,844: | ||
===Library function call=== |
===Library function call=== |
||
{{libheader|Biodem}} |
{{libheader|Biodem}} |
||
< |
<syntaxhighlight lang="rsplus">library(Biodem) |
||
m <- matrix(c(3,2,2,1), nrow=2) |
m <- matrix(c(3,2,2,1), nrow=2) |
||
mtx.exp(m, 0) |
mtx.exp(m, 0) |
||
Line 2,865: | Line 2,865: | ||
# [,1] [,2] |
# [,1] [,2] |
||
# [1,] 1346269 832040 |
# [1,] 1346269 832040 |
||
# [2,] 832040 514229</ |
# [2,] 832040 514229</syntaxhighlight> |
||
Note that non-integer powers are not supported with this function. |
Note that non-integer powers are not supported with this function. |
||
===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: |
||
< |
<syntaxhighlight lang="rsplus">a <- matrix(c(1, 2, 3, 4), 2, 2) |
||
a^1 |
a^1 |
||
a^2</ |
a^2</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>> a^1 |
<pre>> a^1 |
||
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? |
||
< |
<syntaxhighlight 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 |
||
Line 2,902: | Line 2,902: | ||
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</ |
nonSquareMatrix %^% 2#R's %*% will throw the error for us</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>> a %^% 0 |
<pre>> a %^% 0 |
||
Line 2,939: | Line 2,939: | ||
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: |
||
< |
<syntaxhighlight lang="rsplus">library(Biodem) |
||
`%^%` <- function(mat, n) Biodem::mtx.exp(mat, n)</ |
`%^%` <- function(mat, n) Biodem::mtx.exp(mat, n)</syntaxhighlight> |
||
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}}== |
||
<syntaxhighlight lang="racket"> |
|||
<lang Racket> |
|||
#lang racket |
#lang racket |
||
(require math) |
(require math) |
||
Line 2,975: | Line 2,975: | ||
(for ([i (in-range 1 11)]) |
(for ([i (in-range 1 11)]) |
||
(printf "a^~a = ~s\n" i (matrix-expt a i))) |
(printf "a^~a = ~s\n" i (matrix-expt a i))) |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Raku}}== |
=={{header|Raku}}== |
||
(formerly Perl 6) |
(formerly Perl 6) |
||
<lang |
<syntaxhighlight lang="raku" line>subset SqMat of Array where { .elems == all(.[]».elems) } |
||
multi infix:<*>(SqMat $a, SqMat $b) {[ |
multi infix:<*>(SqMat $a, SqMat $b) {[ |
||
Line 3,013: | Line 3,013: | ||
say "### Order $order"; |
say "### Order $order"; |
||
show @m ** $order; |
show @m ** $order; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>### Order 0 |
<pre>### Order 0 |
||
Line 3,090: | Line 3,090: | ||
=={{header|Rust}}== |
=={{header|Rust}}== |
||
Rust (1.37.0) does not allow to overload the ** operator, instead ^ (bitwise xor) is used. |
Rust (1.37.0) does not allow to overload the ** operator, instead ^ (bitwise xor) is used. |
||
< |
<syntaxhighlight lang="rust">use std::fmt; |
||
use std::ops; |
use std::ops; |
||
const WIDTH: usize = 6; |
const WIDTH: usize = 6; |
||
Line 3,164: | Line 3,164: | ||
println!("Power of {}:\n{:?}", i, sm.clone() ^ i); |
println!("Power of {}:\n{:?}", i, sm.clone() ^ i); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 3,224: | Line 3,224: | ||
=={{header|Scala}}== |
=={{header|Scala}}== |
||
< |
<syntaxhighlight lang="scala">class Matrix[T](matrix:Array[Array[T]])(implicit n: Numeric[T], m: ClassManifest[T]) |
||
{ |
{ |
||
import n._ |
import n._ |
||
Line 3,261: | Line 3,261: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>-- m -- |
<pre>-- m -- |
||
Line 3,295: | Line 3,295: | ||
For simplicity, the matrix is represented as a list of lists, and no dimension checking occurs. This implementation does not work when the exponent is 0. |
For simplicity, the matrix is represented as a list of lists, and no dimension checking occurs. This implementation does not work when the exponent is 0. |
||
< |
<syntaxhighlight lang="scheme"> |
||
(define (dec x) |
(define (dec x) |
||
(- x 1)) |
(- x 1)) |
||
Line 3,319: | Line 3,319: | ||
(define (square-matrix mat) |
(define (square-matrix mat) |
||
(matrix-multiply mat mat)) |
(matrix-multiply mat mat)) |
||
</syntaxhighlight> |
|||
</lang> |
|||
Line 3,335: | Line 3,335: | ||
*A ''for'' loop which loops over values listed in an array literal |
*A ''for'' loop which loops over values listed in an array literal |
||
< |
<syntaxhighlight lang="seed7">$ include "seed7_05.s7i"; |
||
include "float.s7i"; |
include "float.s7i"; |
||
Line 3,426: | Line 3,426: | ||
writeln(m ** exponent); |
writeln(m ** exponent); |
||
end for; |
end for; |
||
end func;</ |
end func;</syntaxhighlight> |
||
Original source of matrix exponentiation: [http://seed7.sourceforge.net/algorith/math.htm#matrix_exponentiation] |
Original source of matrix exponentiation: [http://seed7.sourceforge.net/algorith/math.htm#matrix_exponentiation] |
||
Line 3,468: | Line 3,468: | ||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
< |
<syntaxhighlight lang="ruby">class Array { |
||
method ** (Number n { .>= 0 }) { |
method ** (Number n { .>= 0 }) { |
||
var tmp = self |
var tmp = self |
||
Line 3,489: | Line 3,489: | ||
var t = (m ** order) |
var t = (m ** order) |
||
say (' ', t.join("\n ")) |
say (' ', t.join("\n ")) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 3,522: | Line 3,522: | ||
{{works with|OpenAxiom}} |
{{works with|OpenAxiom}} |
||
{{works with|Axiom}} |
{{works with|Axiom}} |
||
< |
<syntaxhighlight lang="spad">(1) -> A:=matrix [[0,-%i],[%i,0]] |
||
+0 - %i+ |
+0 - %i+ |
||
Line 3,545: | Line 3,545: | ||
(4) | | |
(4) | | |
||
+%i 0 + |
+%i 0 + |
||
Type: Union(Matrix(Fraction(Complex(Integer))),...)</ |
Type: Union(Matrix(Fraction(Complex(Integer))),...)</syntaxhighlight> |
||
Domain:[http://fricas.github.io/api/Matrix.html?highlight=matrix Matrix(R)] |
Domain:[http://fricas.github.io/api/Matrix.html?highlight=matrix Matrix(R)] |
||
Line 3,553: | Line 3,553: | ||
This implementation uses [https://en.wikipedia.org/wiki/Exponentiation_by_squaring Exponentiation by squaring] to compute a^n for a matrix a and an integer n (which may be positive, negative or zero). |
This implementation uses [https://en.wikipedia.org/wiki/Exponentiation_by_squaring Exponentiation by squaring] to compute a^n for a matrix a and an integer n (which may be positive, negative or zero). |
||
< |
<syntaxhighlight lang="stata">real matrix matpow(real matrix a, real scalar n) { |
||
real matrix p, x |
real matrix p, x |
||
real scalar i, s |
real scalar i, s |
||
Line 3,565: | Line 3,565: | ||
} |
} |
||
return(s?luinv(p):p) |
return(s?luinv(p):p) |
||
}</ |
}</syntaxhighlight> |
||
Here is an example to compute Fibonacci numbers: |
Here is an example to compute Fibonacci numbers: |
||
< |
<syntaxhighlight lang="stata">: matpow((0,1\1,1),10) |
||
[symmetric] |
[symmetric] |
||
1 2 |
1 2 |
||
Line 3,575: | Line 3,575: | ||
1 | 34 | |
1 | 34 | |
||
2 | 55 89 | |
2 | 55 89 | |
||
+-----------+</ |
+-----------+</syntaxhighlight> |
||
=={{header|Tcl}}== |
=={{header|Tcl}}== |
||
Using code at [[Matrix multiplication#Tcl]] and [[Matrix Transpose#Tcl]] |
Using code at [[Matrix multiplication#Tcl]] and [[Matrix Transpose#Tcl]] |
||
< |
<syntaxhighlight lang="tcl">package require Tcl 8.5 |
||
namespace path {::tcl::mathop ::tcl::mathfunc} |
namespace path {::tcl::mathop ::tcl::mathfunc} |
||
Line 3,602: | Line 3,602: | ||
for {set n 0} {$n < $size} {incr n} {lset i $n $n 1} |
for {set n 0} {$n < $size} {incr n} {lset i $n $n 1} |
||
return $i |
return $i |
||
}</ |
}</syntaxhighlight> |
||
<pre>% print_matrix [matrix_exp {{3 2} {2 1}} 1] |
<pre>% print_matrix [matrix_exp {{3 2} {2 1}} 1] |
||
3 2 |
3 2 |
||
Line 3,625: | Line 3,625: | ||
=={{header|TI-89 BASIC}}== |
=={{header|TI-89 BASIC}}== |
||
Built-in exponentiation: |
Built-in exponentiation: |
||
< |
<syntaxhighlight lang="ti89b">[3,2;4,1]^4</syntaxhighlight> |
||
Output: <math>\begin{bmatrix}417 & 208 \\ 416 & 209\end{bmatrix}</math> |
Output: <math>\begin{bmatrix}417 & 208 \\ 416 & 209\end{bmatrix}</math> |
||
Line 3,631: | Line 3,631: | ||
For matrices of floating point numbers, the library function <code>mmult</code> can be used as shown. The user-defined <code>id</code> function takes a square matrix to the identity matrix of the same dimensions. The <code>mex</code> function takes a pair <math>(A,n)</math> |
For matrices of floating point numbers, the library function <code>mmult</code> can be used as shown. The user-defined <code>id</code> function takes a square matrix to the identity matrix of the same dimensions. The <code>mex</code> function takes a pair <math>(A,n)</math> |
||
representing a real matrix <math>A</math> and a natural exponent <math>n</math> to the exponentiation <math>A^n</math> using the naive algorithm. |
representing a real matrix <math>A</math> and a natural exponent <math>n</math> to the exponentiation <math>A^n</math> using the naive algorithm. |
||
< |
<syntaxhighlight lang="ursala">#import nat |
||
#import lin |
#import lin |
||
id = @h ^|CzyCK33/1.! 0.!* |
id = @h ^|CzyCK33/1.! 0.!* |
||
mex = ||id@l mmult:-0^|DlS/~& iota</ |
mex = ||id@l mmult:-0^|DlS/~& iota</syntaxhighlight> |
||
Alternatively, this version uses the fast binary algorithm. |
Alternatively, this version uses the fast binary algorithm. |
||
< |
<syntaxhighlight lang="ursala">mex = ~&ar^?\id@al (~&lr?/mmult@llPrX ~&r)^/~&alrhPX mmult@falrtPXPRiiX</syntaxhighlight> |
||
This test program raises a 2 by 2 matrix to a selection of powers. |
This test program raises a 2 by 2 matrix to a selection of powers. |
||
< |
<syntaxhighlight lang="ursala">#cast %eLLL |
||
test = mex/*<<3.,2.>,<2.,1.>> <0,1,2,3,4,10></ |
test = mex/*<<3.,2.>,<2.,1.>> <0,1,2,3,4,10></syntaxhighlight> |
||
output: |
output: |
||
<pre>< |
<pre>< |
||
Line 3,665: | Line 3,665: | ||
=={{header|VBA}}== |
=={{header|VBA}}== |
||
No operator overloading in VBA. Implemented as a function. Can not handle scalars. Requires matrix size greater than one. Does allow for negative exponents. |
No operator overloading in VBA. Implemented as a function. Can not handle scalars. Requires matrix size greater than one. Does allow for negative exponents. |
||
< |
<syntaxhighlight lang="vb">Option Base 1 |
||
Private Function Identity(n As Integer) As Variant |
Private Function Identity(n As Integer) As Variant |
||
Dim I() As Variant |
Dim I() As Variant |
||
Line 3,720: | Line 3,720: | ||
Debug.Print |
Debug.Print |
||
pp MatrixExponentiation(M3, 10) |
pp MatrixExponentiation(M3, 10) |
||
End Sub</ |
End Sub</syntaxhighlight>{{out}} |
||
<pre>-1 2 |
<pre>-1 2 |
||
2 -3 |
2 -3 |
||
Line 3,740: | Line 3,740: | ||
The Matrix class in the above module also has a 'pow' method but, as an alternative, overloads the otherwise unused '^' operator to provide the same functionality. |
The Matrix class in the above module also has a 'pow' method but, as an alternative, overloads the otherwise unused '^' operator to provide the same functionality. |
||
< |
<syntaxhighlight lang="ecmascript">import "/matrix" for Matrix |
||
import "/fmt" for Fmt |
import "/fmt" for Fmt |
||
Line 3,747: | Line 3,747: | ||
Fmt.mprint(m, 2, 0) |
Fmt.mprint(m, 2, 0) |
||
System.print("\nRaised to power of 10:\n") |
System.print("\nRaised to power of 10:\n") |
||
Fmt.mprint(m ^ 10, 3, 0)</ |
Fmt.mprint(m ^ 10, 3, 0)</syntaxhighlight> |
||
{{out}} |
{{out}} |