Matrix-exponentiation operator: Difference between revisions

{{trans|Python}}

 

assert(m1[0].len == m2.len)

V r = [[0] * m2[0].len] * m1.len

 

print("\n10:")

 

{{out}}

=={{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^2ⁿ:

procedure Test_Matrix is

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);

Sample output:

<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.

with Ada.Complex_Text_IO; use Ada.Complex_Text_IO;

with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types;

Put_Line ("M**1 ="); Put (M**(1.0,0.0));

Put_Line ("M**0.5 ="); Put (M**(0.5,0.0));

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.

 

{{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'''

MODE VEC = [default upb]COSCAL;

MODE MAT = [default upb,default upb]COSCAL;

OD;

out

BITS binary exponent:=BIN exponent ;

MAT out := IF bits width ELEM binary exponent THEN base ELSE IDENTITY UPB base FI;

OD;

out

 

MODE COSCAL = COMPL;

printf(($" mat ** "g(0)":"l$,24));

compl mat printf(scal fmt, mat**24);

Output:

<pre>

 

=={{header|BBC BASIC}}==

matrix() = 3, 2, 2, 1

NEXT

ENDIF

Output:

<pre>

 

Matrix multiplication is a known idiom taken from BQN crate. Matrix exponentiation is simply doing Matrix multiplication n times.

 

MatEx ← {𝕨 MatMul⍟(𝕩-1) 𝕨}

 

(>⟨3‿2

· ┌─ ┌─ ┌─ ┌─ ┌─

╵ 3 2 ╵ 13 8 ╵ 55 34 ╵ 233 144 ╵ 1346269 832040

2 1 8 5 34 21 144 89 832040 514229

┘ ┘ ┘ ┘ ┘

 

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.

 

 

=={{header|Burlesque}}==

 

=={{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.

#include <stdio.h>

#include <stdlib.h>

 

return 0;

Output:

<pre>m0 dim:3 =

 

=={{header|C sharp}}==

using System.Collections;

using System.Collections.Generic;

}

 

{{out}}

<pre style="height:30ex;overflow:scroll">

=={{header|C++}}==

This is an implementation in C++.

#include <cmath>

#include <iostream>

}

return d;

This is the task part.

Mx operator^(int n) {

if (n < 0)

 

return 0;

{{out}}

<pre>

 

This uses the '*' operator for arrays as defined in [[Matrix_multiplication#Chapel]]

// create result matrix of same dimensions

var r:[a.domain] a.eltType;

 

return r;

 

Usage example (like Perl):

m(1,1) = 1; m(1,2) = 2; m(1,3) = 0;

m(2,1) = 0; m(2,2) = 3; m(2,3) = 1;

writeln("Order ", i);

writeln(m ** i, "\n");

 

{{out}}

=={{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.

"Takes two 2D arrays and returns their product, or an error if they cannot be multiplied"

(let* ((m0-dims (array-dimensions matrix-0))

(multiply-matrices me2 me2)))

(t (let ((me2 (matrix-expt matrix (/ (1- exp) 2))))

Output (note that this lisp implementation uses single-precision floats for decimals by default). We can also use rationals:

CL-USER> (setf 5x5-matrix

 

=={{header|D}}==

 

struct SquareMat(T = creal) {

foreach (immutable p; [0, 1, 23, 24])

writefln("m ^^ %d =\n%s", p, m ^^ p);

{{out}}

<pre>m ^^ 0 =

 

=={{header|Delphi}}==

program Matrix_exponentiation_operator;

 

Readln;

end.

{{out}}

<pre>

| 2 1 |

</pre>

 

!$MATRIX

END FOR

 

Sample output:

<pre>

There is already a built-in word (<code>m^n</code>) that implements exponentiation. Here is a simple and less efficient implementation.

 

 

: my-m^n ( m n -- m' )

[ drop length identity-matrix ]

[ swap '[ _ m. ] times ] 2bi

 

( scratchpad ) { { 3 2 } { 2 1 } } 0 my-m^n .

=={{header|Fermat}}==

Matrix exponentiation for square matrices and integer powers is built in.

Array a[2,2]; {illustrate with a 2x2 matrix}

[a]:=[(2/3, 1/3, 4/5, 1/5)];

[a]^3;

[a]^10;

{{out}}

<pre>

=={{header|Fortran}}==

{{works with|Fortran|90 and later}}

implicit none

end do

 

Output

<pre> 1.00000 0.00000 0.00000

This operator performs M^n for any square invertible matrix M and integer n, including negative powers.

 

#include once "rowech.bas"

#include once "matinv.bas"

next i

print

{{out}}

<pre> 308.9999999999998 -307.9999999999998

 

=={{header|GAP}}==

A := [[0 , 1], [1, 1]];

PrintArray(A);

PrintArray(A^10);

# [ [ 34, 55 ],

 

=={{header|Go}}==

<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.

 

import "fmt"

fmt.Println()

}

 

{{out}}

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).

 

 

(<+>)

-- TEST ----------------------------------------------------------------------

main :: IO ()

{{Out}}

<pre>Mat [[1,8],[0,1]]</pre>

===With Numeric.LinearAlgebra===

 

 

a :: Matrix I

print $ a^4

putStrLn "power of zero: "

{{Out}}

<pre>

 

 

 

or, from [[j:Essays/Linear Recurrences|the J wiki]], and faster for large exponents:

 

 

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).

 

Extends [[Matrix Transpose#JavaScript]] and [[Matrix multiplication#JavaScript]]

function IdentityMatrix(n) {

this.height = n;

 

var m = new Matrix([[3, 2], [2, 1]]);

output

<pre>1,0

 

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.

def indicator(i;n): [range(0;n) | 0] | .[i] = 1;

 

| multiply($ans; $residue )

end

'''Examples'''

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.

 

def rotation_matrix(theta):

 

def demo_direct_matrix_exp(n):

'''Results''':

$ time jq -n -c -f Matrix-exponentiation_operator.rc

[[1,-1.1102230246251565e-12],[1.1102230246251565e-12,1]]

user 0m0.490s

$ time jq -n -c -f Matrix-exponentiation_operator.rc

[[1,-7.849831895612169e-13],[7.849831895612169e-13,1]]

user 0m0.625s

 

=={{header|Jsish}}==

Uses module listed in [[Matrix Transpose#Jsish]]. Fails the task spec actually, as Matrix.exp() is implemented as a method, not an operator.

 

require('Matrix');

 

m.exp(10) ==> { height:2, mtx:[ [ 1346269, 832040 ], [ 832040, 514229 ] ], width:2 }

=!EXPECTEND!=

 

{{out}}

=={{header|Julia}}==

Matrix exponentiation is implemented by the built-in <code>^</code> operator.

2x2 Array{Int64,2}:

89 55

 

=={{header|K}}==

/Matrix Exponentiation

/mpow.k

pow: {:[0=y; :({a=/:a:!x}(#x))];a: x; do[y-1; a: x _mul a]; :a}

 

The output of a session is given below:

{{out}}

 

=={{header|Kotlin}}==

 

typealias Vector = DoubleArray

)

for (i in 0..10) printMatrix(m pow i, i)

 

{{out}}

 

=={{header|Lambdatalk}}==

{require lib_matrix}

 

M^4 = [[233,144],[144,89]]

M^10 = [[1346269,832040],[832040,514229]]

 

=={{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.

MatrixD$ ="3, 3, 0.86603, 0.50000, 0.00000, -0.50000, 0.86603, 0.00000, 0.00000, 0.00000, 1.00000"

 

MatrixE$ =MatrixToPower$( MatrixD$, 9)

call DisplayMatrix MatrixE$

 

{{out}}

 

=={{header|Lua}}==

 

function Matrix.new( dim_y, dim_x )

n = m^4;

 

 

=={{header|M2000 Interpreter}}==

Module CheckIt {

Class cArray {

}

Checkit

 

{{out}}

=={{header|Maple}}==

Maple handles matrix powers implicitly with the built-in exponentiation operator:

<math>\left[\begin{array}{cc}

7 & 15 \\

 

If you want elementwise powers, you can use the elementwise <code>^~</code> operator:

<math>\left[\begin{array}{cc}

1 & 9 \\

=={{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.

MatrixPower[a, 0]

MatrixPower[a, 1]

MatrixPower[a, 4]

MatrixPower[a, 1/2]

gives back:

 

 

Symbolic matrices like {{i,j},{k,l}} to the power m give general solutions for all possible i,j,k,l, and m:

gives back (note that the simplification is not necessary for the evaluation, it just gives a shorter output):

 

=={{header|MATLAB}}==

For exponents in the form of A*A*A*A*...*A, A must be a square matrix:

 

Otherwise, to take the individual array elements to the power of an exponent (the matrix need not be square):

 

=={{header|Maxima}}==

[4, 1])$

 

a ^^ -1;

/* matrix([-1/5, 2/5],

 

=={{header|Nim}}==

 

type Matrix[N: static int; T] = array[1..N, array[1..N, T]]

S30 = 1 / 2

let m2: Matrix[2, float] = [[C30, -S30], [S30, C30]] # 30° rotation matrix.

 

{{out}}

We will use some auxiliary functions

 

let eye n =

let a = Array.make_matrix n n 0.0 in

(* example with integers *)

pow 1 ( * ) 2 16;;

 

Now matrix power is simply a special case of pow :

 

let p, q = dim a in

if p <> q then failwith "bad dimensions" else

 

[| [| 1.0; 1.0|]; [| 1.0; 0.0 |] |] ^^ 10;;

 

=={{header|Octave}}==

Of course GNU Octave handles matrix and operations on matrix "naturally".

 

M^0

M^1

M^2

M^(-1)

 

Output:

 

=={{header|PARI/GP}}==

 

=={{header|Perl}}==

package SquareMatrix;

use Carp; # standard, "it's not my fault" module

print "\n### WAY too big:\n", $m ** 1_000_000_000_000;

print "\n### But identity matrix can handle that\n",

 

=={{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>

First two routines copied straight from the [[Identity_matrix#Phix|Identity_matrix]] and [[Matrix_multiplication#Phix|Matrix_multiplication]] tasks.

<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: #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>

{{out}}

<pre>

=={{header|PicoLisp}}==

Uses the 'matMul' function from [[Matrix multiplication#PicoLisp]]

(let L (need N (1) 0)

(mapcar '(() (copy (rot L))) L) ) )

M ) )

 

Output:

<pre>-> ((55 34) (34 21))</pre>

=={{header|Python}}==

Using matrixMul from [[Matrix multiplication#Python]]

>>> def matrixMul(m1, m2):

return map(

1346269 832040

832040 514229

 

Alternative Based Upon @ operator of Python 3.5 PEP 465 and using Matrix exponentation for faster computation of powers

class Mat(list) :

def __matmul__(self, B) :

print('\n%i:' % i)

printtable(power(m, i))

{{Output}}

<pre>

===Library function call===

{{libheader|Biodem}}

m <- matrix(c(3,2,2,1), nrow=2)

mtx.exp(m, 0)

# [,1] [,2]

# [1,] 1346269 832040

Note that non-integer powers are not supported with this function.

===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:

a^1

{{out}}

<pre>> a^1

[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?

{

is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol#See the docs for is.integer

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

nonSquareMatrix %^% 1

{{out}}

<pre>> a %^% 0

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:

And it will work just the same, except for being much faster and throwing an error on nonSquareMatrix %^% 1.

 

=={{header|Racket}}==

 

#lang racket

(require math)

(for ([i (in-range 1 11)])

(printf "a^~a = ~s\n" i (matrix-expt a i)))

 

=={{header|Raku}}==

(formerly Perl 6)

 

multi infix:<*>(SqMat $a, SqMat $b) {[

say "### Order $order";

show @m ** $order;

{{out}}

<pre>### Order 0

=={{header|Rust}}==

Rust (1.37.0) does not allow to overload the ** operator, instead ^ (bitwise xor) is used.

use std::ops;

const WIDTH: usize = 6;

println!("Power of {}:\n{:?}", i, sm.clone() ^ i);

}

{{out}}

<pre>

 

=={{header|Scala}}==

{

import n._

}

}

{{out}}

<pre>-- m --

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.

 

(define (dec x)

(- x 1))

(define (square-matrix mat)

(matrix-multiply mat mat))

 

 

*A ''for'' loop which loops over values listed in an array literal

 

include "float.s7i";

 

writeln(m ** exponent);

end for;

 

Original source of matrix exponentiation: [http://seed7.sourceforge.net/algorith/math.htm#matrix_exponentiation]

 

=={{header|Sidef}}==

method ** (Number n { .>= 0 }) {

var tmp = self

var t = (m ** order)

say (' ', t.join("\n "))

{{out}}

<pre>

{{works with|OpenAxiom}}

{{works with|Axiom}}

 

+0 - %i+

(4) | |

+%i 0 +

 

Domain:[http://fricas.github.io/api/Matrix.html?highlight=matrix Matrix(R)]

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).

 

real matrix p, x

real scalar i, s

}

return(s?luinv(p):p)

 

Here is an example to compute Fibonacci numbers:

 

[symmetric]

1 2

1 | 34 |

2 | 55 89 |

 

=={{header|Tcl}}==

Using code at [[Matrix multiplication#Tcl]] and [[Matrix Transpose#Tcl]]

namespace path {::tcl::mathop ::tcl::mathfunc}

 

for {set n 0} {$n < $size} {incr n} {lset i $n $n 1}

return $i

<pre>% print_matrix [matrix_exp {{3 2} {2 1}} 1]

3 2

=={{header|TI-89 BASIC}}==

Built-in exponentiation:

Output: <math>\begin{bmatrix}417 & 208 \\ 416 & 209\end{bmatrix}</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.

#import lin

 

id = @h ^|CzyCK33/1.! 0.!*

Alternatively, this version uses the fast binary algorithm.

This test program raises a 2 by 2 matrix to a selection of powers.

 

output:

<pre><

=={{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.

Private Function Identity(n As Integer) As Variant

Dim I() As Variant

Debug.Print

pp MatrixExponentiation(M3, 10)

<pre>-1 2

2 -3

 

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.

import "/fmt" for Fmt

 

Fmt.mprint(m, 2, 0)

System.print("\nRaised to power of 10:\n")

 

{{out}}