Matrix-exponentiation operator

Matrix-exponentiation operator
You are encouraged to solve this task according to the task description, using any language you may know.

Most programming languages have a built-in implementation of exponentiation for integers and reals only.

Task

Demonstrate how to implement matrix exponentiation as an operator.

11l

Translation of: Python

F matrix_mul(m1, m2)
   assert(m1[0].len == m2.len)
   V r = [[0] * m2[0].len] * m1.len
   L(j) 0 .< m1.len
      L(i) 0 .< m2[0].len
         V s = 0
         L(k) 0 .< m2.len
            s += m1[j][k] * m2[k][i]
         r[j][i] = s
   R r

F identity(size)
   V rsize = 0 .< size
   R rsize.map(j -> @rsize.map(i -> Int(i == @j)))

F matrixExp(m, pow)
   assert(pow >= 0 & Int(pow) == pow, ‘Only non-negative, integer powers allowed’)
   V accumulator = identity(m.len)
   L(i) 0 .< pow
      accumulator = matrix_mul(accumulator, m)
   R accumulator

F printtable(data)
   L(row) data
      print(row.map(cell -> ‘#<5’.format(cell)).join(‘ ’))

V m = [[3, 2], [2, 1]]
L(i) 5
   print("\n#.:".format(i))
   printtable(matrixExp(m, i))

print("\n10:")
printtable(matrixExp(m, 10))

Output:

0:
1     0    
0     1    

1:
3     2    
2     1    

2:
13    8    
8     5    

3:
55    34   
34    21   

4:
233   144  
144   89   

10:
1346269 832040
832040 514229

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²ⁿ:

with Ada.Text_IO;  use Ada.Text_IO;
 
procedure Test_Matrix is
   generic
      type Element is private;
      Zero : Element;
      One  : Element;
      with function "+" (A, B : Element) return Element is <>;
      with function "*" (A, B : Element) return Element is <>;
      with function Image (X : Element) return String is <>;
   package Matrices is
      type Matrix is array (Integer range <>, Integer range <>) of Element;
      function "*" (A, B : Matrix) return Matrix;
      function "**" (A : Matrix; Power : Natural) return Matrix;
      procedure Put (A : Matrix);
   end Matrices;

   package body Matrices is
      function "*" (A, B : Matrix) return Matrix is
         R   : Matrix (A'Range (1), B'Range (2));
         Sum : Element := Zero;
      begin
         for I in R'Range (1) loop
            for J in R'Range (2) loop
               Sum := Zero;
               for K in A'Range (2) loop
                  Sum := Sum + A (I, K) * B (K, J);
               end loop;
               R (I, J) := Sum;
            end loop;
         end loop;
         return R;
      end "*";

      function "**" (A : Matrix; Power : Natural) return Matrix is
      begin
         if Power = 1 then
            return A;
         end if;
         declare
            R : Matrix (A'Range (1), A'Range (2)) := (others => (others => Zero));
            P : Matrix  := A;
            E : Natural := Power;
         begin
            for I in P'Range (1) loop -- R is identity matrix
               R (I, I) := One;
            end loop;
            if E = 0 then
               return R;
            end if;
            loop
               if E mod 2 /= 0 then
                  R := R * P;
               end if;
               E := E / 2;
               exit when E = 0;
               P := P * P;
            end loop;
            return R;
         end;
      end "**";
      
      procedure Put (A : Matrix) is
      begin
         for I in A'Range (1) loop
            for J in A'Range (1) loop
               Put (Image (A (I, J)));
            end loop;
            New_Line;
         end loop;
      end Put;
   end Matrices;
   
   package Integer_Matrices is new Matrices (Integer, 0, 1, Image => Integer'Image);
   use Integer_Matrices;
   
   M : Matrix (1..2, 1..2) := ((3,2),(2,1));
begin
   Put_Line ("M =");       Put (M);
   Put_Line ("M**0 =");    Put (M**0);
   Put_Line ("M**1 =");    Put (M**1);
   Put_Line ("M**2 =");    Put (M**2);
   Put_Line ("M*M =");     Put (M*M);
   Put_Line ("M**3 =");    Put (M**3);
   Put_Line ("M*M*M =");   Put (M*M*M);
   Put_Line ("M**4 =");    Put (M**4);
   Put_Line ("M*M*M*M ="); Put (M*M*M*M);
   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);
end Test_Matrix;

Sample output:

M =
 3 2
 2 1
M**0 =
 1 0
 0 1
M**1 =
 3 2
 2 1
M**2 =
 13 8
 8 5
M*M =
 13 8
 8 5
M**3 =
 55 34
 34 21
M*M*M =
 55 34
 34 21
M**4 =
 233 144
 144 89
M*M*M*M =
 233 144
 144 89
M**10 =
 1346269 832040
 832040 514229
M*M*M*M*M*M*M*M*M*M =
 1346269 832040
 832040 514229

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.Text_IO;                  use Ada.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.Real_Arrays;     use Ada.Numerics.Real_Arrays;
with Ada.Numerics.Complex_Arrays;  use Ada.Numerics.Complex_Arrays;
with Ada.Numerics.Complex_Elementary_Functions; use Ada.Numerics.Complex_Elementary_Functions;

procedure Test_Matrix is
   function "**" (A : Complex_Matrix; Power : Complex) return Complex_Matrix is
      L  : Real_Vector (A'Range (1));
      X  : Complex_Matrix (A'Range (1), A'Range (2));
      R  : Complex_Matrix (A'Range (1), A'Range (2));
      RL : Complex_Vector (A'Range (1));
   begin
      Eigensystem (A, L, X);
      for I in L'Range loop
         RL (I) := (L (I), 0.0) ** Power;
      end loop;
      for I in R'Range (1) loop
         for J in R'Range (2) loop
            declare
               Sum : Complex := (0.0, 0.0);
            begin
               for K in RL'Range (1) loop
                  Sum := Sum + X (I, K) * RL (K) * X (J, K);
               end loop;
               R (I, J) := Sum;
            end;
         end loop;
      end loop;
      return R;
   end "**";
   procedure Put (A : Complex_Matrix) is
   begin
      for I in A'Range (1) loop
        for J in A'Range (2) loop
           Put (A (I, J));
        end loop;
        New_Line;
      end loop;
   end Put;
   M : Complex_Matrix (1..2, 1..2) := (((3.0,0.0),(2.0,1.0)),((2.0,-1.0),(1.0,0.0)));
begin
   Put_Line ("M =");      Put (M);
   Put_Line ("M**0 =");   Put (M**(0.0,0.0));
   Put_Line ("M**1 =");   Put (M**(1.0,0.0));
   Put_Line ("M**0.5 ="); Put (M**(0.5,0.0));
end Test_Matrix;

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.

(Another person is talking here:) I have made small corrections and tested this in 2023, and it did not work as I expected. However, I have questions about the mathematical libraries. I tried both GCC 12 and GCC 13. (I also tried the last GNAT Community Edition, but it no longer functions on my system.) What might be needed here is one's own eigensystem routine.

On the other hand, I did get a version working to raise a real matrix to a natural number power, thus demonstrating the correctness of the approach:

with Ada.Text_IO;                  use Ada.Text_IO;
with Ada.Float_Text_IO;            use Ada.Float_Text_IO;
with Ada.Numerics.Real_Arrays;     use Ada.Numerics.Real_Arrays;

procedure Test_Matrix is
   procedure Put (A : Real_Matrix) is
   begin
      for I in A'Range (1) loop
        for J in A'Range (2) loop
           Put (" ");
           Put (A (I, J));
        end loop;
        New_Line;
      end loop;
   end Put;
   function "**" (A : Real_Matrix; Power : Integer) return Real_Matrix is
      L  : Real_Vector (A'Range (1));
      X  : Real_Matrix (A'Range (1), A'Range (2));
      R  : Real_Matrix (A'Range (1), A'Range (2));
      RL : Real_Vector (A'Range (1));
   begin
      Eigensystem (A, L, X);
      for I in L'Range loop
         RL (I) := L (I) ** Power;
      end loop;
      for I in R'Range (1) loop
         for J in R'Range (2) loop
            declare
               Sum : Float := 0.0;
            begin
               for K in RL'Range loop
                  Sum := Sum + X (I, K) * RL (K) * X (J, K);
               end loop;
               R (I, J) := Sum;
            end;
         end loop;
      end loop;
      return R;
   end "**";
   M : Real_Matrix (1..2, 1..2) := ((3.0, 2.0), (2.0, 1.0));
begin
   Put_Line ("M =");      Put (M);
   Put_Line ("M**0 =");   Put (M**0);
   Put_Line ("M**1 =");   Put (M**1);
   Put_Line ("M**2 =");   Put (M**2);
   Put_Line ("M**3 =");   Put (M**3);
   Put_Line ("M**50 =");  Put (M**50);
end Test_Matrix;

Output:

M =
  3.00000E+00  2.00000E+00
  2.00000E+00  1.00000E+00
M**0 =
  1.00000E+00  0.00000E+00
  0.00000E+00  1.00000E+00
M**1 =
  3.00000E+00  2.00000E+00
  2.00000E+00  1.00000E+00
M**2 =
  1.30000E+01  8.00000E+00
  8.00000E+00  5.00000E+00
M**3 =
  5.50000E+01  3.40000E+01
  3.40000E+01  2.10000E+01
M**50 =
  1.61305E+31  9.96919E+30
  9.96919E+30  6.16130E+30

ALGOL 68

Works with: ALGOL 68 version Revision 1 - no extensions to language used.

Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny.

File: Matrix_algebra.a68

INT default upb=3;
MODE VEC = [default upb]COSCAL;
MODE MAT = [default upb,default upb]COSCAL;

OP * = (VEC a,b)COSCAL: (
    COSCAL result:=0;
    FOR i FROM LWB a TO UPB a DO result+:= a[i]*b[i] OD;
    result
  );

OP * = (VEC a, MAT b)VEC: ( # overload vec times matrix #
    [2 LWB b:2 UPB b]COSCAL result;
    FOR j FROM 2 LWB b TO 2 UPB b DO result[j]:=a*b[,j] OD;
    result
  );

OP * = (MAT a, b)MAT: ( # overload matrix times matrix #
    [LWB a:UPB a, 2 LWB b:2 UPB b]COSCAL result;
    FOR k FROM LWB result TO UPB result DO result[k,]:=a[k,]*b OD;
    result
  );

OP IDENTITY = (INT upb)MAT:(
  [upb,upb] COSCAL out;
  FOR i TO upb DO
    FOR j TO upb DO
      out[i,j]:= ( i=j |1|0)
    OD
  OD;
  out
);

File: Matrix-exponentiation_operator.a68

OP ** = (MAT base, INT exponent)MAT: (
  BITS binary exponent:=BIN exponent ;
  MAT out := IF bits width ELEM binary exponent THEN base ELSE IDENTITY UPB base FI;
  MAT sq:=base;

  WHILE
    binary exponent := binary exponent SHR 1;
    binary exponent /= BIN 0
  DO
    sq := sq * sq;
    IF bits width ELEM binary exponent THEN out := out * sq FI
  OD;
  out
);

File: test_Matrix-exponentiation_operator.a68

#!/usr/local/bin/a68g --script #

MODE COSCAL = COMPL;
PR READ "Matrix_algebra.a68" PR
PR READ "Matrix-exponentiation_operator.a68" PR

PROC compl mat printf= (FORMAT scal fmt, MAT m)VOID:(
  FORMAT
    vec math = $n(2 UPB m)(f(scal fmt)"&")$,
    mat math = $"<math>\begin{bmat}"ln(UPB m)(xxf(vec fmt)"\\"l)"\end{bmat}</math>"$,
    vec fmt = $"("n(2 UPB m-1)(f(scal fmt)",")f(scal fmt)")"$,
    mat fmt = $x"("n(UPB m-1)(f(vec fmt)","lxx)f(vec fmt)");"$;
  # finally print the result #
  printf((mat fmt,m))
);

FORMAT scal fmt = $-d.dddd,+d.dddd"i"$; # width of 4, with no leading '+' sign, 1 decimals #
MAT mat=((sqrt(0.5)I0         , sqrt(0.5)I0        , 0I0),
         (        0I-sqrt(0.5),         0Isqrt(0.5), 0I0),
         (        0I0         ,         0I0        , 0I1))

printf(($" mat ** "g(0)":"l$,24));
compl mat printf(scal fmt, mat**24);
print(newline)

Output:

 mat ** 24:
 (( 1.0000+0.0000i, 0.0000+0.0000i, 0.0000+0.0000i),
  ( 0.0000+0.0000i, 1.0000+0.0000i, 0.0000+0.0000i),
  ( 0.0000+0.0000i, 0.0000+0.0000i, 1.0000+0.0000i));

ATS

(* I will write a GENERAL template for raising something to a
   non-negative integer power, and then apply that template to matrix
   multiplication. *)

#include "share/atspre_staload.hats"

(*------------------------------------------------------------------*)
(* The interface. *)

extern fn {a : t@ype} nonnegative_integer_power : (a, intGte 0) -> a
extern fn {a : t@ype} zeroth_power : () -> a
extern fn {a : t@ype} product : (a, a) -> a

(*------------------------------------------------------------------*)
(* The implementation of "nonnegative_integer_power". *)

(* I use the squaring method. See
https://en.wikipedia.org/w/index.php?title=Exponentiation_by_squaring&oldid=1144956501
*)

implement {a}
nonnegative_integer_power (M, i) =
  let
    fun
    repeat {i : nat}     (* <-- This number consistently shrinks. *)
           .<i>.         (* <-- Proof the recursion will terminate. *)
           (Accum : a,   (* "Accumulator" *)
            Base  : a,
            i     : int i)
        : a =
      if i = 0 then
        Accum
      else
        let
          val i_halved = half i (* Integer division. *)
          and Base_squared = product<a> (Base, Base)
        in
          if i_halved + i_halved = i then
            repeat (Accum, Base_squared, i_halved)
          else
            repeat (product<a> (Base, Accum), Base_squared, i_halved)
        end
  in
    repeat (zeroth_power<a> (), M, i)
  end

(*------------------------------------------------------------------*)
(* Application of nonnegative_integer_power to mtrxszref. *)

fn {tk : tkind}
npow_mtrxszref (M : mtrxszref (g0float tk),
                p : intGte 0)
    : mtrxszref (g0float tk) =
  let
    typedef a = g0float tk

    val n = mtrxszref_get_nrow M
    val () =
      if mtrxszref_get_ncol M <> n then
        $raise IllegalArgExn ("npow_mtrxszref:matrix_not_square")

    implement
    zeroth_power<mtrxszref a> () =
      (* Return an n-by-n identity matrix. *)
      let
        val I = mtrxszref_make_elt<a> (n, n, g0i2f 0)
        var k : Size_t
      in
        for (k := i2sz 0; k <> n; k := succ k)
          I[k, k] := g0i2f 1;
        I
      end

    implement
    product<mtrxszref a> (A, B) =
      (* Return the matrix product of A and B. *)
      let
        val C = mtrxszref_make_elt<a> (n, n, g0i2f 0)
        var i : Size_t
      in
        for (i := i2sz 0; i <> n; i := succ i)
          let
            var j : Size_t
          in
            for (j := i2sz 0; j <> n; j := succ j)
              let
                var k : Size_t
              in
                for (k := i2sz 0; k <> n; k := succ k)
                  C[i, j] := C[i, j] + (A[i, k] * B[k, j])
              end
          end;
        C
      end
  in
    nonnegative_integer_power<mtrxszref a> (M, p)
  end

overload ** with npow_mtrxszref

(*------------------------------------------------------------------*)

implement
main0 () =
  let
    (* This matrix is borrowed from the entry for the programming
       language Chapel:
       
          1 2 0
          0 3 1
          1 0 0

    *)
    val A = mtrxszref_make_elt (i2sz 3, i2sz 3, 0.0)
    val () = A[0, 0] := 1.0
    val () = A[0, 1] := 2.0
    val () = A[1, 1] := 3.0
    val () = A[1, 2] := 1.0
    val () = A[2, 0] := 1.0

    var p : intGte 0
  in
    for (p := 0; p <> 11; p := succ p)
      let
        val B = A ** p
      in
        fprint_val<string> (stdout_ref, "power = ");
        fprint_val<int> (stdout_ref, p);
        fprint_val<string> (stdout_ref, "\n");
        fprint_mtrxszref_sep<double> (stdout_ref, B, "\t", "\n");
        fprint_val<string> (stdout_ref, "\n\n")
      end
  end

(*------------------------------------------------------------------*)

Output:

$ patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW matrix_exponentiation_task.dats -lgc && ./a.out
power = 0
1.000000	0.000000	0.000000
0.000000	1.000000	0.000000
0.000000	0.000000	1.000000

power = 1
1.000000	2.000000	0.000000
0.000000	3.000000	1.000000
1.000000	0.000000	0.000000

power = 2
1.000000	8.000000	2.000000
1.000000	9.000000	3.000000
1.000000	2.000000	0.000000

power = 3
3.000000	26.000000	8.000000
4.000000	29.000000	9.000000
1.000000	8.000000	2.000000

power = 4
11.000000	84.000000	26.000000
13.000000	95.000000	29.000000
3.000000	26.000000	8.000000

power = 5
37.000000	274.000000	84.000000
42.000000	311.000000	95.000000
11.000000	84.000000	26.000000

power = 6
121.000000	896.000000	274.000000
137.000000	1017.000000	311.000000
37.000000	274.000000	84.000000

power = 7
395.000000	2930.000000	896.000000
448.000000	3325.000000	1017.000000
121.000000	896.000000	274.000000

power = 8
1291.000000	9580.000000	2930.000000
1465.000000	10871.000000	3325.000000
395.000000	2930.000000	896.000000

power = 9
4221.000000	31322.000000	9580.000000
4790.000000	35543.000000	10871.000000
1291.000000	9580.000000	2930.000000

power = 10
13801.000000	102408.000000	31322.000000
15661.000000	116209.000000	35543.000000
4221.000000	31322.000000	9580.000000

BBC BASIC

      DIM matrix(1,1), output(1,1)
      matrix() = 3, 2, 2, 1
      
      FOR power% = 0 TO 9
        PROCmatrixpower(matrix(), output(), power%)
        PRINT "matrix()^" ; power% " = "
        FOR row% = 0 TO DIM(output(), 1)
          FOR col% = 0 TO DIM(output(), 2)
            PRINT output(row%,col%);
          NEXT
          PRINT
        NEXT row%
      NEXT power%
      END
      
      DEF PROCmatrixpower(src(), dst(), pow%)
      LOCAL i%
      dst() = 0
      FOR i% = 0 TO DIM(dst(), 1) : dst(i%,i%) = 1 : NEXT
      IF pow% THEN
        FOR i% = 1 TO pow%
          dst() = dst() . src()
        NEXT
      ENDIF
      ENDPROC

Output:

matrix()^0 =
         1         0
         0         1
matrix()^1 =
         3         2
         2         1
matrix()^2 =
        13         8
         8         5
matrix()^3 =
        55        34
        34        21
matrix()^4 =
       233       144
       144        89
matrix()^5 =
       987       610
       610       377
matrix()^6 =
      4181      2584
      2584      1597
matrix()^7 =
     17711     10946
     10946      6765
matrix()^8 =
     75025     46368
     46368     28657
matrix()^9 =
    317811    196418
    196418    121393

BQN

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

MatMul ← +˝∘×⎉1‿∞

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

(>⟨3‿2
   2‿1⟩) MatEx 1‿2‿3‿4‿10

┌─                                                           
· ┌─      ┌─       ┌─        ┌─          ┌─                  
  ╵ 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.

MatEx ← MatMul{𝔽´𝔽˜⍟(/2|⌊∘÷⟜2⍟(↕1+·⌊2⋆⁼⊢)𝕩)𝕨}

Burlesque

blsq ) {{1 1} {1 0}} 10 .*{mm}r[
{{89 55} {55 34}}

C

C doesn't support classes or allow operator overloading. The following is code that defines a function, SquareMtxPower that will raise a matrix to a positive integer power.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

typedef struct squareMtxStruct {
    int   dim;
    double *cells;
    double **m;
} *SquareMtx;

/* function for initializing row r of a new matrix */
typedef void (*FillFunc)( double *cells, int r, int dim, void *ff_data);

SquareMtx NewSquareMtx( int dim, FillFunc fillFunc, void *ff_data ) 
{
    SquareMtx sm = malloc(sizeof(struct squareMtxStruct));
    if (sm) {
        int rw;
        sm->dim = dim;
        sm->cells = malloc(dim*dim * sizeof(double));
        sm->m = malloc( dim * sizeof(double *));
        if ((sm->cells != NULL) && (sm->m != NULL)) {
            for (rw=0; rw<dim; rw++) {
                sm->m[rw] = sm->cells + dim*rw;
                fillFunc( sm->m[rw], rw, dim, ff_data );
            }
        }
        else {
            free(sm->m);
            free(sm->cells);
            free(sm);
            printf("Square Matrix allocation failure\n");
            return NULL;
        }
    }
    else {
        printf("Malloc failed for square matrix\n");
    }
    return sm;
}

void ffMatxSquare( double *cells, int rw, int dim, SquareMtx m0 )
{
    int col, ix;
    double sum;
    double *m0rw = m0->m[rw];
    
    for (col = 0; col < dim; col++) {
        sum = 0.0;
        for (ix=0; ix<dim; ix++)
            sum += m0rw[ix] * m0->m[ix][col];
        cells[col] = sum;
    }
}

void ffMatxMulply( double *cells, int rw, int dim, SquareMtx mplcnds[] )
{
    SquareMtx mleft = mplcnds[0];
    SquareMtx mrigt = mplcnds[1];
    double sum;
    double *m0rw = mleft->m[rw];
    int col, ix;

    for (col = 0; col < dim; col++) {
        sum = 0.0;
        for (ix=0; ix<dim; ix++)
            sum += m0rw[ix] * mrigt->m[ix][col];
        cells[col] = sum;
    }
}

void MatxMul( SquareMtx mr, SquareMtx left, SquareMtx rigt)
{
    int rw;
    SquareMtx mplcnds[2];
    mplcnds[0] = left; mplcnds[1] = rigt;

    for (rw = 0; rw < left->dim; rw++)  
        ffMatxMulply( mr->m[rw], rw, left->dim, mplcnds);
}

void ffIdentity( double *cells, int rw, int dim, void *v )
{
    int col;
    for (col=0; col<dim; col++) cells[col] = 0.0;
    cells[rw] = 1.0;
}
void ffCopy(double *cells, int rw, int dim, SquareMtx m1)
{
    int col;
    for (col=0; col<dim; col++) cells[col] = m1->m[rw][col];
}

void FreeSquareMtx( SquareMtx m ) 
{
    free(m->m);
    free(m->cells);
    free(m);
}

SquareMtx SquareMtxPow( SquareMtx m0, int exp )
{
    SquareMtx v0 = NewSquareMtx(m0->dim, ffIdentity, NULL);
    SquareMtx v1 = NULL;
    SquareMtx base0 = NewSquareMtx( m0->dim, ffCopy, m0);
    SquareMtx base1 = NULL;
    SquareMtx mplcnds[2], t;

    while (exp) {
        if (exp % 2) {
            if (v1)
                MatxMul( v1, v0, base0);
            else  {
                mplcnds[0] = v0; mplcnds[1] = base0;
                v1 = NewSquareMtx(m0->dim, ffMatxMulply, mplcnds); 
            }
            {t = v0; v0=v1; v1 = t;}
        }
        if (base1)
            MatxMul( base1, base0, base0);
        else 
            base1 = NewSquareMtx( m0->dim, ffMatxSquare, base0);
        t = base0; base0 = base1; base1 = t;
        exp = exp/2;
    }
    if (base0) FreeSquareMtx(base0);
    if (base1) FreeSquareMtx(base1);
    if (v1) FreeSquareMtx(v1);
    return v0;
}

FILE *fout;
void SquareMtxPrint( SquareMtx mtx, const char *mn ) 
{
    int rw, col;
    int d = mtx->dim;

    fprintf(fout, "%s dim:%d =\n", mn, mtx->dim);

    for (rw=0; rw<d; rw++) {
        fprintf(fout, " |");
        for(col=0; col<d; col++) 
            fprintf(fout, "%8.5f ",mtx->m[rw][col] );
        fprintf(fout, " |\n");
    }
    fprintf(fout, "\n");
}

void fillInit( double *cells, int rw, int dim, void *data)
{
    double theta = 3.1415926536/6.0;
    double c1 = cos( theta);
    double s1 = sin( theta);

    switch(rw) {
    case 0:
        cells[0]=c1; cells[1]=s1; cells[2]=0.0;
        break;
    case 1:
        cells[0]=-s1; cells[1]=c1; cells[2]=0;
        break;
    case 2:
        cells[0]=0.0; cells[1]=0.0; cells[2]=1.0;
        break;
    }
}

int main()
{
    SquareMtx m0 = NewSquareMtx( 3, fillInit, NULL);
    SquareMtx m1 = SquareMtxPow( m0, 5);
    SquareMtx m2 = SquareMtxPow( m0, 9);
    SquareMtx m3 = SquareMtxPow( m0, 2);

//  fout = stdout;
    fout = fopen("matrx_exp.txt", "w");
    SquareMtxPrint(m0, "m0"); FreeSquareMtx(m0);
    SquareMtxPrint(m1, "m0^5"); FreeSquareMtx(m1);
    SquareMtxPrint(m2, "m0^9"); FreeSquareMtx(m2);
    SquareMtxPrint(m3, "m0^2"); FreeSquareMtx(m3);
    fclose(fout);

    return 0;
}

Output:

m0 dim:3 =
 | 0.86603  0.50000  0.00000  |
 |-0.50000  0.86603  0.00000  |
 | 0.00000  0.00000  1.00000  |

m0^5 dim:3 =
 |-0.86603  0.50000  0.00000  |
 |-0.50000 -0.86603  0.00000  |
 | 0.00000  0.00000  1.00000  |

m0^9 dim:3 =
 | 0.00000 -1.00000  0.00000  |
 | 1.00000  0.00000  0.00000  |
 | 0.00000  0.00000  1.00000  |

m0^2 dim:3 =
 | 0.50000  0.86603  0.00000  |
 |-0.86603  0.50000  0.00000  |
 | 0.00000  0.00000  1.00000  |

C#

using System;
using System.Collections;
using System.Collections.Generic;
using static System.Linq.Enumerable;

public static class MatrixExponentation
{
    public static double[,] Identity(int size) {
        double[,] matrix = new double[size, size];
        for (int i = 0; i < size; i++) matrix[i, i] = 1;
        return matrix;
    }

    public static double[,] Multiply(this double[,] left, double[,] right) {
        if (left.ColumnCount() != right.RowCount()) throw new ArgumentException();
        double[,] m = new double[left.RowCount(), right.ColumnCount()];
        foreach (var (row, column) in from r in Range(0, m.RowCount()) from c in Range(0, m.ColumnCount()) select (r, c)) {
            m[row, column] = Range(0, m.RowCount()).Sum(i => left[row, i] * right[i, column]);
        }
        return m;
    }

    public static double[,] Pow(this double[,] matrix, int exp) {
        if (matrix.RowCount() != matrix.ColumnCount()) throw new ArgumentException("Matrix must be square.");
        double[,] accumulator = Identity(matrix.RowCount());
        for (int i = 0; i < exp; i++) {
            accumulator = accumulator.Multiply(matrix);
        }
        return accumulator;
    }

    private static int RowCount(this double[,] matrix) => matrix.GetLength(0);
    private static int ColumnCount(this double[,] matrix) => matrix.GetLength(1);

    private static void Print(this double[,] m) {
        foreach (var row in Rows()) {
            Console.WriteLine("[ " + string.Join("   ", row) + " ]");
        }
        Console.WriteLine();

        IEnumerable<IEnumerable<double>> Rows() =>
            Range(0, m.RowCount()).Select(row => Range(0, m.ColumnCount()).Select(column => m[row, column]));
    }

    public static void Main() {
        var matrix = new double[,] {
            { 3, 2 },
            { 2, 1 }
        };
        
        matrix.Pow(0).Print();
        matrix.Pow(1).Print();
        matrix.Pow(2).Print();
        matrix.Pow(3).Print();
        matrix.Pow(4).Print();
        matrix.Pow(50).Print();
    }

}

Output:

[ 1   0 ]
[ 0   1 ]

[ 3   2 ]
[ 2   1 ]

[ 13   8 ]
[ 8   5 ]

[ 55   34 ]
[ 34   21 ]

[ 233   144 ]
[ 144   89 ]

[ 1.61305314249046E+31   9.9692166771893E+30 ]
[ 9.9692166771893E+30   6.16131474771528E+30 ]

C++

This is an implementation in C++.

#include <complex>
#include <cmath>
#include <iostream>
using namespace std;

template<int MSize = 3, class T = complex<double> >
class SqMx {
  typedef T Ax[MSize][MSize];
  typedef SqMx<MSize, T> Mx;

private:
  Ax a;
  SqMx() { }

public:
  SqMx(const Ax &_a) { // constructor with pre-defined array
    for (int r = 0; r < MSize; r++)
      for (int c = 0; c < MSize; c++)
        a[r][c] = _a[r][c];
  }

  static Mx identity() {
    Mx m;
    for (int r = 0; r < MSize; r++)
      for (int c = 0; c < MSize; c++)
        m.a[r][c] = (r == c ? 1 : 0);
    return m;
  }

  friend ostream &operator<<(ostream& os, const Mx &p)
  { // ugly print
    for (int i = 0; i < MSize; i++) {
      for (int j = 0; j < MSize; j++)
        os << p.a[i][j] << ',';
      os << endl;
    }
    return os;
  }

  Mx operator*(const Mx &b) {
    Mx d;
    for (int r = 0; r < MSize; r++)
      for (int c = 0; c < MSize; c++) {
        d.a[r][c] = 0;
        for (int k = 0; k < MSize; k++)
          d.a[r][c] += a[r][k] * b.a[k][c];
      }
    return d;
  }

This is the task part.

  // C++ does not have a ** operator, instead, ^ (bitwise Xor) is used.
  Mx operator^(int n) {
    if (n < 0)
      throw "Negative exponent not implemented";

    Mx d = identity();
    for (Mx sq = *this; n > 0; sq = sq * sq, n /= 2)
      if (n % 2 != 0)
        d = d * sq;
    return d;
  } 
};

typedef SqMx<> M3;
typedef complex<double> creal;

int main() {
  double q = sqrt(0.5);
  creal array[3][3] = { { { q,  0 }, { q, 0 }, { 0, 0 } },
                        { { 0, -q }, { 0, q }, { 0, 0 } },
                        { { 0,  0 }, { 0, 0 }, { 0, 1 } } };
  M3 m(array);

  cout << "m ^ 23=" << endl
       << (m ^ 23) << endl;

  return 0;
}

Output:

m ^ 23=
(0.707107,0),(0,0.707107),(0,0),
(0.707107,0),(0,-0.707107),(0,0),
(0,0),(0,0),(0,-1),

An alternative way would be to implement operator*= and conversion from number (giving multiples of the identity matrix) for the matrix and use the generic code from Exponentiation operator#C++ with support for negative exponents removed (or alternatively, implement matrix inversion as well, implement /= in terms of it, and use the generic code unchanged). Note that the algorithm used there is much faster as well.

Chapel

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

proc **(a, e) {
    // create result matrix of same dimensions
    var r:[a.domain] a.eltType;
    // and initialize to identity matrix
    forall ij in r.domain do
        r(ij) = if ij(1) == ij(2) then 1 else 0;

    for 1..e do
        r *= a;

    return r;
}

Usage example (like Perl):

var m:[1..3, 1..3] int;
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(3,1) = 1; m(3,2) = 0; m(3,3) = 0;

config param n = 10;

for i in 0..n do {
    writeln("Order ", i);
    writeln(m ** i, "\n");
}

Output:

Order 0
1 0 0
0 1 0
0 0 1

Order 1
1 2 0
0 3 1
1 0 0

Order 2
1 8 2
1 9 3
1 2 0

Order 3
3 26 8
4 29 9
1 8 2

Order 4
11 84 26
13 95 29
3 26 8

Order 5
37 274 84
42 311 95
11 84 26

Order 6
121 896 274
137 1017 311
37 274 84

Order 7
395 2930 896
448 3325 1017
121 896 274

Order 8
1291 9580 2930
1465 10871 3325
395 2930 896

Order 9
4221 31322 9580
4790 35543 10871
1291 9580 2930

Order 10
13801 102408 31322
15661 116209 35543
4221 31322 9580

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.

(defun multiply-matrices (matrix-0 matrix-1)
  "Takes two 2D arrays and returns their product, or an error if they cannot be multiplied"
  (let* ((m0-dims (array-dimensions matrix-0))
         (m1-dims (array-dimensions matrix-1))
         (m0-dim (length m0-dims))
         (m1-dim (length m1-dims)))
    (if (or (/= 2 m0-dim) (/= 2 m1-dim))
        (error "Array given not a matrix")
        (let ((m0-rows (car m0-dims))
              (m0-cols (cadr m0-dims))
              (m1-rows (car m1-dims))
              (m1-cols (cadr m1-dims)))
          (if (/= m0-cols m1-rows)
              (error "Incompatible dimensions")
              (do ((rarr (make-array (list m0-rows m1-cols)
                                     :initial-element 0) rarr)
                   (n 0 (if (= n (1- m0-cols)) 0 (1+ n)))
                   (cc 0 (if (= n (1- m0-cols))
                             (if (/= cc (1- m1-cols)) 
                                 (1+ cc) 0) cc))
                   (cr 0 (if (and (= (1- m0-cols) n) 
                                  (= (1- m1-cols) cc))
                             (1+ cr)
                             cr)))
                  ((= cr m0-rows) rarr)
                (setf (aref rarr cr cc)
                      (+ (aref rarr cr cc)
                         (* (aref matrix-0 cr n)
                            (aref matrix-1 n cc))))))))))

(defun matrix-identity (dim)
  "Creates a new identity matrix of size dim*dim"
  (do ((rarr (make-array (list dim dim)
                         :initial-element 0) rarr)
       (n 0 (1+ n)))
      ((= n dim) rarr)
    (setf (aref rarr n n) 1)))

(defun matrix-expt (matrix exp)
  "Takes the first argument (a matrix) and multiplies it by itself exp times"
  (let* ((m-dims (array-dimensions matrix))
         (m-rows (car m-dims))
         (m-cols (cadr m-dims)))
    (cond
      ((/= m-rows m-cols) (error "Non-square matrix"))
      ((zerop exp) (matrix-identity m-rows))
      ((= 1 exp) (do ((rarr (make-array (list m-rows m-cols)) rarr)
                      (cc 0 (if (= cc (1- m-cols))
                                0
                                (1+ cc)))
                      (cr 0 (if (= cc (1- m-cols))
                                (1+ cr)
                                cr)))
                     ((= cr m-rows) rarr)
                   (setf (aref rarr cr cc) (aref matrix cr cc))))
      ((zerop (mod exp 2)) (let ((me2 (matrix-expt matrix (/ exp 2))))
                             (multiply-matrices me2 me2)))
      (t (let ((me2 (matrix-expt matrix (/ (1- exp) 2))))
           (multiply-matrices matrix (multiply-matrices me2 me2)))))))

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

CL-USER> (setf 5x5-matrix
               (make-array '(5 5)
                           :initial-contents
                           '((0    1 -1   -2    2)
                             (0.4  4  3.2 -3  -10)
                             (4.5 -2  0.5  1    7)
                             (10   1  0    1.5 -2)
                             (4    5 -3   -2    1))))
#2A((0 1 -1 -2 2)
    (0.4 4 3.2 -3 -10)
    (4.5 -2 0.5 1 7)
    (10 1 0 1.5 -2)
    (4 5 -3 -2 1))
CL-USER> (matrix-expt 5x5-matrix 3)
#2A((-163.25 -19.5 92.25 -7.5999985 -184.3)
    (156.6 -412.09998 0.7999954 331.45 597.4)
    (-129.82501 401.25 -66.975 -302.55 -390.15)
    (-148.9 39.25 -5.200001 -67.225006 -7.300003)
    (-495.05 -231.5 310.85 33.0 -328.5))
CL-USER> (setf 4x4-matrix
               (make-array '(4 4)
                           :initial-contents
                           '(( 1/2 -1/2  4    8)
                             (-3/4  7/3  8/5 -2)
                             (-5   17   20/3 -5/2)
                             ( 3/2 -1   -7/3  6))))                            
#2A((1/2 -1/2 4 8) (-3/4 7/3 8/5 -2) (-5 17 20/3 -5/2) (3/2 -1 -7/3 6))
CL-USER> (matrix-expt 4x4-matrix 3)
#2A((-233/8 182723/720 757/30 353/6)
    (-73517/480 838241/2160 77789/450 -67537/180)
    (-5315/9 66493/45 90883/135 -54445/36)
    (37033/144 -27374/45 -15515/54 12109/18))

D

import std.stdio, std.string, std.math, std.array, std.algorithm;

struct SquareMat(T = creal) {
    public static string fmt = "%8.3f";
    private alias TM = T[][];
    private TM a;

    public this(in size_t side) pure nothrow @safe
    in {
        assert(side > 0);
    } body {
        a = new TM(side, side);
    }

    public this(in TM m) pure nothrow @safe
    in {
        assert(!m.empty);
        assert(m.all!(row => row.length == m.length)); // Is square.
    } body {
        // 2D dup.
        a.length = m.length;
        foreach (immutable i, const row; m)
            a[i] = row.dup;
    }

    string toString() const @safe {
        return format("<%(%(" ~ fmt ~ ", %)\n %)>", a);
    }

    public static SquareMat identity(in size_t side) pure nothrow @safe {
        auto m = SquareMat(side);
        foreach (immutable r, ref row; m.a)
            foreach (immutable c; 0 .. side)
                row[c] = (r == c) ? 1+0i : 0+0i;
        return m;
    }

    public SquareMat opBinary(string op:"*")(in SquareMat other)
    const pure nothrow @safe in {
        assert (a.length == other.a.length);
    } body {
        immutable side = other.a.length;
        auto d = SquareMat(side);
        foreach (immutable r; 0 .. side)
            foreach (immutable c; 0 .. side) {
                d.a[r][c] = 0+0i;
                foreach (immutable k, immutable ark; a[r])
                    d.a[r][c] += ark * other.a[k][c];
            }
        return d;
    }

    public SquareMat opBinary(string op:"^^")(int n) // The task part.
    const pure nothrow @safe in {
        assert(n >= 0, "Negative exponent not implemented.");
    } body {
        auto sq = SquareMat(this.a);
        auto d = SquareMat.identity(a.length);
        for (; n > 0; sq = sq * sq, n >>= 1)
            if (n & 1)
                d = d * sq;
        return d;
    }
}

void main() {
    alias M = SquareMat!();
    enum real q = 0.5.sqrt;
    immutable m = M([[   q + 0*1.0Li,    q + 0*1.0Li, 0.0L + 0.0Li],
                     [0.0L - q*1.0Li, 0.0L + q*1.0Li, 0.0L + 0.0Li],
                     [0.0L +   0.0Li, 0.0L +   0.0Li, 0.0L + 1.0Li]]);
    M.fmt = "%5.2f";
    foreach (immutable p; [0, 1, 23, 24])
        writefln("m ^^ %d =\n%s", p, m ^^ p);
}

Output:

m ^^ 0 =
< 1.00+ 0.00i,  0.00+ 0.00i,  0.00+ 0.00i
  0.00+ 0.00i,  1.00+ 0.00i,  0.00+ 0.00i
  0.00+ 0.00i,  0.00+ 0.00i,  1.00+ 0.00i>
m ^^ 1 =
< 0.71+ 0.00i,  0.71+ 0.00i,  0.00+ 0.00i
  0.00+-0.71i,  0.00+ 0.71i,  0.00+ 0.00i
  0.00+ 0.00i,  0.00+ 0.00i,  0.00+ 1.00i>
m ^^ 23 =
< 0.71+ 0.00i,  0.00+ 0.71i,  0.00+ 0.00i
  0.71+ 0.00i,  0.00+-0.71i,  0.00+ 0.00i
  0.00+ 0.00i,  0.00+ 0.00i,  0.00+-1.00i>
m ^^ 24 =
< 1.00+ 0.00i,  0.00+ 0.00i,  0.00+ 0.00i
  0.00+ 0.00i,  1.00+ 0.00i,  0.00+ 0.00i
  0.00+ 0.00i,  0.00+ 0.00i,  1.00+ 0.00i>

Delphi

program Matrix_exponentiation_operator;

{$APPTYPE CONSOLE}

{$R *.res}

uses
  System.SysUtils;

type
  TCells = array of array of double;

  TMatrix = record
  private
    FCells: TCells;
    function GetCells(r, c: Integer): Double;
    procedure SetCells(r, c: Integer; const Value: Double);
    class operator Implicit(a: TMatrix): string;
    class operator BitwiseXor(a: TMatrix; e: Integer): TMatrix;
    class operator Multiply(a: TMatrix; b: TMatrix): TMatrix;
  public
    constructor Create(w, h: integer); overload;
    constructor Create(c: TCells); overload;
    constructor Ident(size: Integer);
    function Rows: Integer;
    function Columns: Integer;
    property Cells[r, c: Integer]: Double read GetCells write SetCells; default;
  end;

{ TMatrix }

constructor TMatrix.Create(c: TCells);
begin
  Create(Length(c), Length(c[0]));
  FCells := c;
end;

constructor TMatrix.Create(w, h: integer);
begin
  SetLength(FCells, w, h);
end;

class operator TMatrix.BitwiseXor(a: TMatrix; e: Integer): TMatrix;
begin
  if e < 0 then
    raise Exception.Create('Matrix inversion not implemented');

  Result.Ident(a.Rows);
  while e > 0 do
  begin
    Result := Result * a;
    dec(e);
  end;
end;

function TMatrix.Rows: Integer;
begin
  Result := Length(FCells);
end;

function TMatrix.Columns: Integer;
begin
  Result := 0;
  if Rows > 0 then
    Result := Length(FCells);
end;

function TMatrix.GetCells(r, c: Integer): Double;
begin
  Result := FCells[r, c];
end;

constructor TMatrix.Ident(size: Integer);
var
  i: Integer;
begin
  Create(size, size);

  for i := 0 to size - 1 do
    Cells[i, i] := 1;
end;

class operator TMatrix.Implicit(a: TMatrix): string;
var
  i, j: Integer;
begin
  Result := '[';
  if a.Rows > 0 then
    for i := 0 to a.Rows - 1 do
    begin
      if i > 0 then
        Result := Trim(Result) + ']'#10'[';
      for j := 0 to a.Columns - 1 do
      begin
        Result := Result + Format('%f', [a[i, j]]) + ' ';
      end;
    end;
  Result := trim(Result) + ']';
end;

class operator TMatrix.Multiply(a, b: TMatrix): TMatrix;
var
  size: Integer;
  r: Integer;
  c: Integer;
  k: Integer;
begin
  if (a.Rows <> b.Rows) or (a.Columns <> b.Columns) then
    raise Exception.Create('The matrix must have same size');

  size := a.Rows;
  Result.Create(size, size);

  for r := 0 to size - 1 do
    for c := 0 to size - 1 do
    begin
      Result[r, c] := 0;
      for k := 0 to size - 1 do
        Result[r, c] := Result[r, c] + a[r, k] * b[k, c];
    end;
end;

procedure TMatrix.SetCells(r, c: Integer; const Value: Double);
begin
  FCells[r, c] := Value;
end;

var
  M: TMatrix;

begin
  M.Create([[3, 2], [2, 1]]);
// Delphi don't have a ** and can't override ^ operator, then XOR operator was used 
  Writeln(string(M xor 0), #10);
  Writeln(string(M xor 1), #10);
  Writeln(string(M xor 2), #10);
  Writeln(string(M xor 3), #10);
  Writeln(string(M xor 4), #10);
  Writeln(string(M xor 50), #10);
  Readln;
end.

Output:

[1,00 0,00]
[0,00 1,00]

[3,00 2,00]
[2,00 1,00]

[13,00 8,00]
[8,00 5,00]

[55,00 34,00]
[34,00 21,00]

[233,00 144,00]
[144,00 89,00]

[1,61305314249045832E31 9,96921667718930453E30]
[9,96921667718930115E30 6,16131474771527643E30]

ERRE

                               10
This example calculates | 3 2 |
                        | 2 1 |

PROGRAM MAT_PROD

!$MATRIX

!-----------------
! calculate A[]^N
!-----------------

CONST ORDER=1

DIM A[1,1],B[1,1],ANS[1,1]

BEGIN

DATA(3,2,2,1)
DATA(10)  ! integer power only

FOR I=0 TO ORDER DO
   FOR J=0 TO ORDER DO
      READ(A[I,J])
   END FOR
END FOR

READ(M)  N=M-1

IF N=0 THEN   ! A[]^0=matrice identit…
   for I=0 TO ORDER DO
      B[I,I]=1
   END FOR
 ELSE
   B[]=A[]
   FOR Z=1 TO N DO
      ANS[]=0
      FOR I=0 TO ORDER DO
         FOR J=0 TO ORDER DO
            FOR K=0 TO ORDER DO
               ANS[I,J]=ANS[I,J]+(A[I,K]*B[K,J])
            END FOR
         END FOR
      END FOR
      B[]=ANS[]
  END FOR
END IF

! print answer
  FOR I=0 TO ORDER DO
     FOR J=0 TO ORDER DO
        PRINT(B[I,J],)
     END FOR
     PRINT
  END FOR

END PROGRAM

Sample output:

 1346269   832040
 832040    514229

Factor

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

USING: kernel math math.matrices sequences ;

: my-m^n ( m n -- m' )
    dup 0 < [ "no negative exponents" throw ] [
        [ drop length identity-matrix ]
        [ swap '[ _ m. ] times ] 2bi
    ] if ;

( scratchpad ) { { 3 2 } { 2 1 } } 0 my-m^n .
{ { 1 0 } { 0 1 } }
( scratchpad ) { { 3 2 } { 2 1 } } 4 my-m^n .
{ { 233 144 } { 144 89 } }

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]^-1;                       {matrix inverse}
[a]^0;                        {identity matrix}
[a]^2;
[a]^3;
[a]^10;

Output:

[[  -3 / 2, 6, `
    5 / 2, -5   ]]

[[  1,  0, `
    0,  1   ]]

[[  32 / 45, 52 / 75, `
    13 / 45, 23 / 75   ]]

[[  476 / 675, 796 / 1125, `
    199 / 675, 329 / 1125   ]]

[[  81409466972 / 115330078125, 135682444612 / 192216796875, `
    33920611153 / 115330078125, 56534352263 / 192216796875   ]]

Fortran

Works with: Fortran version 90 and later

module matmod
  implicit none
   
! Overloading the ** operator does not work because the compiler cannot
! differentiate between matrix exponentiation and the elementwise raising
! of an array to a power therefore we define a new operator
  interface operator (.matpow.)
    module procedure matrix_exp
  end interface

contains

function matrix_exp(m, n) result (res)
  real, intent(in)  :: m(:,:)
  integer, intent(in)  :: n
  real :: res(size(m,1),size(m,2))
  integer :: i
   
  if(n == 0) then
    res = 0
    do i = 1, size(m,1)
      res(i,i) = 1
    end do
    return
  end if

  res = m
  do i = 2, n
    res = matmul(res, m)
  end do
  
end function matrix_exp
end module matmod

program Matrix_exponentiation
  use matmod
  implicit none

  integer, parameter :: n = 3
  real, dimension(n,n) :: m1, m2
  integer :: i, j
  
  m1 = reshape((/ (i, i = 1, n*n) /), (/ n, n /), order = (/ 2, 1 /))
 
  do i = 0, 4
    m2 = m1 .matpow. i
    do j = 1, size(m2,1)
      write(*,*) m2(j,:)
    end do
    write(*,*)
  end do

end program Matrix_exponentiation

Output

      1.00000         0.00000         0.00000
      0.00000         1.00000         0.00000
      0.00000         0.00000         1.00000
 
      1.00000         2.00000         3.00000
      4.00000         5.00000         6.00000
      7.00000         8.00000         9.00000
 
      30.0000         36.0000         42.0000
      66.0000         81.0000         96.0000
      102.000         126.000         150.000
 
      468.000         576.000         684.000
      1062.00         1305.00         1548.00
      1656.00         2034.00         2412.00
 
      7560.00         9288.00         11016.0
      17118.0         21033.0         24948.0
      26676.0         32778.0         38880.0

FreeBASIC

The include statements incorporate the code from Matrix multiplication#FreeBASIC, which defines the Matrix type and the matrix multiplication operator, Reduced row echelon form#FreeBASIC which contains a function for getting a matrix into row-echelon form, and Gauss-Jordan matrix inversion#FreeBASIC which gives the inverse of a matrix. Make sure to remove all the print statements first though.

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

#include once "matmult.bas"
#include once "rowech.bas"
#include once "matinv.bas"

operator ^ (byval M as Matrix, byval n as integer ) as Matrix
    dim as uinteger i, j, k = ubound( M.m, 1 )
    if n < 0 then return matinv(M) ^ (-n)
    if n = 0 then return M * matinv(M)
    return (M ^ (n-1)) * M
end operator

dim as Matrix M = Matrix(2,2), Q
dim as integer i, j, n
M.m(0,0) = 1./3 : M.m(0,1) = 2./3
M.m(1,0) = 2./7 : M.m(1,1) = 5./7

for n = -2 to 4
    Q = (M ^ n)
    for i = 0 to 1
        for j = 0 to 1
            print Q.m(i, j),
        next j
        print
    next i
    print
next n

Output:

 308.9999999999998          -307.9999999999998          
-132           133          

 14.99999999999999          -13.99999999999999          
-6.000000000000003           7.000000000000004          

 1             0            
 0             1            

 0.3333333333333333          0.6666666666666666         
 0.2857142857142857          0.7142857142857143         

 0.3015873015873016          0.6984126984126984         
 0.2993197278911565          0.7006802721088435         

 0.3000755857898715          0.6999244142101284         
 0.299967606090055           0.7000323939099449         

 0.3000035993233272          0.6999964006766727         
 0.2999984574328597          0.7000015425671401

GAP

# Matrix exponentiation is built-in
A := [[0 , 1], [1, 1]];
PrintArray(A);
#   [ [  0,  1 ],
#     [  1,  1 ] ]
PrintArray(A^10);
#   [ [  34,  55 ],
#     [  55,  89 ] ]

Go

Translation of: Kotlin

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.

package main

import "fmt"

type vector = []float64
type matrix []vector

func (m1 matrix) mul(m2 matrix) matrix {
    rows1, cols1 := len(m1), len(m1[0])
    rows2, cols2 := len(m2), len(m2[0])
    if cols1 != rows2 {
        panic("Matrices cannot be multiplied.")
    }
    result := make(matrix, rows1)
    for i := 0; i < rows1; i++ {
        result[i] = make(vector, cols2)
        for j := 0; j < cols2; j++ {
            for k := 0; k < rows2; k++ {
                result[i][j] += m1[i][k] * m2[k][j]
            }
        }
    }
    return result
}

func identityMatrix(n int) matrix {
    if n < 1 {
        panic("Size of identity matrix can't be less than 1")
    }
    ident := make(matrix, n)
    for i := 0; i < n; i++ {
        ident[i] = make(vector, n)
        ident[i][i] = 1
    }
    return ident
}

func (m matrix) pow(n int) matrix {
    le := len(m)
    if le != len(m[0]) {
        panic("Not a square matrix")
    }
    switch {
    case n < 0:
        panic("Negative exponents not supported")
    case n == 0:
        return identityMatrix(le)
    case n == 1:
        return m
    }
    pow := identityMatrix(le)
    base := m
    e := n
    for e > 0 {
        if (e & 1) == 1 {
            pow = pow.mul(base)
        }
        e >>= 1
        base = base.mul(base)
    }
    return pow
}

func main() {
    m := matrix{{3, 2}, {2, 1}}
    for i := 0; i <= 10; i++ {
        fmt.Println("** Power of", i, "**")
        fmt.Println(m.pow(i))
        fmt.Println()
    }
}

Output:

** Power of 0 **
[[1 0] [0 1]]

** Power of 1 **
[[3 2] [2 1]]

** Power of 2 **
[[13 8] [8 5]]

** Power of 3 **
[[55 34] [34 21]]

** Power of 4 **
[[233 144] [144 89]]

** Power of 5 **
[[987 610] [610 377]]

** Power of 6 **
[[4181 2584] [2584 1597]]

** Power of 7 **
[[17711 10946] [10946 6765]]

** Power of 8 **
[[75025 46368] [46368 28657]]

** Power of 9 **
[[317811 196418] [196418 121393]]

** Power of 10 **
[[1.346269e+06 832040] [832040 514229]]

Haskell

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

import Data.List (transpose)

(<+>)
  :: Num a
  => [a] -> [a] -> [a]
(<+>) = zipWith (+)

(<*>)
  :: Num a
  => [a] -> [a] -> a
(<*>) = (sum .) . zipWith (*)

newtype Mat a =
  Mat [[a]]
  deriving (Eq, Show)

instance Num a =>
         Num (Mat a) where
  negate (Mat x) = Mat $ map (map negate) x
  Mat x + Mat y = Mat $ zipWith (<+>) x y
  Mat x * Mat y =
    Mat
      [ [ xs Main.<*> ys -- Main prefix to distinguish fron applicative operator
        | ys <- transpose y ]
      | xs <- x ]
  abs = undefined
  fromInteger _ = undefined -- don't know dimension of the desired matrix
  signum = undefined

-- TEST ----------------------------------------------------------------------
main :: IO ()
main = print $ Mat [[1, 2], [0, 1]] ^ 4

Output:

Mat [[1,8],[0,1]]

This will work for matrices over any numeric type, including complex numbers. The implementation of (^) uses the fast binary algorithm for exponentiation.

Note: this implementation does not work for a power of 0.

With Numeric.LinearAlgebra

import Numeric.LinearAlgebra

a :: Matrix I
a = (2><2) 
  [1,2
  ,0,1]

main = do
  print $ a^4
  putStrLn "power of zero: "
  print $ a^0

Output:

(2><2)
 [ 1, 16
 , 0,  1 ]
power of zero: 
(1><1)
 [ 1 ]

J

mp=: +/ .*   NB. Matrix multiplication 
pow=: pow0=: 4 : 'mp&x^:y =i.#x'

or, from the J wiki, and faster for large exponents:

pow=: pow1=: 4 : 'mp/ mp~^:(I.|.#:y) x'

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

Example use:

    (3 2,:2 1) pow 3
 55 34
 34 21

JavaScript

Works with: SpiderMonkey

for the print() and Array.forEach() functions.

Extends Matrix Transpose#JavaScript and Matrix multiplication#JavaScript

// IdentityMatrix is a "subclass" of Matrix
function IdentityMatrix(n) {
    this.height = n;
    this.width = n;
    this.mtx = [];
    for (var i = 0; i < n; i++) {
        this.mtx[i] = [];
        for (var j = 0; j < n; j++) {
            this.mtx[i][j] = (i == j ? 1 : 0);
        }
    }
}
IdentityMatrix.prototype = Matrix.prototype;

// the Matrix exponentiation function
// returns a new matrix
Matrix.prototype.exp = function(n) {
    var result = new IdentityMatrix(this.height);
    for (var i = 1; i <= n; i++) {
        result = result.mult(this);
    }
    return result;
}

var m = new Matrix([[3, 2], [2, 1]]);
[0,1,2,3,4,10].forEach(function(e){print(m.exp(e)); print()})

output

1,0
0,1

3,2
2,1

13,8
8,5

55,34
34,21

233,144
144,89

1346269,832040
832040,514229

jq

In this section we define matrix_exp(n) for computing the n-th power of the input matrix, where it is assumed that n is a non-negative integer.

The implementation here can be used with any matrix multiplication function, multiply(A;B), for example as defined at Matrix_multiplication#jq. Thus matrix_exp(n) could be used with complex-valued matrices.

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.

# produce an array of length n that is 1 at i and 0 elsewhere
def indicator(i;n): [range(0;n) | 0] | .[i] = 1;

# Identity matrix:
def identity(n): reduce range(0;n) as $i ([]; . + [indicator( $i; n )] );

def direct_matrix_exp(n):
  . as $in
  | if n == 0 then identity($in|length)
    else reduce range(1;n) as $i ($in; . as $m | multiply($m; $in))
    end;

def matrix_exp(n):
  if n < 4 then direct_matrix_exp(n)
  else . as $in
  | ((n|2)|floor) as $m
  | matrix_exp($m) as $ans
  | multiply($ans;$ans) as $ans
  | (n - (2 * $m) ) as $residue
  | if $residue == 0 then $ans
    else matrix_exp($residue) as $residue
    | multiply($ans; $residue )
    end
  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 pi: 4 * (1|atan);

def rotation_matrix(theta):
  [[(theta|cos), (theta|sin)], [-(theta|sin), (theta|cos)]];

def demo_matrix_exp(n):
  rotation_matrix( pi / 4 ) | matrix_exp(n) ;

def demo_direct_matrix_exp(n):
  rotation_matrix( pi / 4 ) | direct_matrix_exp(n) ;

Results:

# For demo_matrix_exp(10000)
$ time jq -n -c -f Matrix-exponentiation_operator.rc
[[1,-1.1102230246251565e-12],[1.1102230246251565e-12,1]]
user	0m0.490s
sys	0m0.008s

# For demo_direct_matrix_exp(10000)
$ time jq -n -c -f Matrix-exponentiation_operator.rc
[[1,-7.849831895612169e-13],[7.849831895612169e-13,1]]
user	0m0.625s
sys	0m0.006s

Jsish

Based on Javascript matrix entries.

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

/* Matrix exponentiation, in Jsish */
require('Matrix');

if (Interp.conf('unitTest')) {
    var m = new Matrix([[3, 2], [2, 1]]);
;    m;
;    m.exp(0);
;    m.exp(1);
;    m.exp(2);
;    m.exp(4);
;    m.exp(10);
}

/*
=!EXPECTSTART!=
m ==> { height:2, mtx:[ [ 3, 2 ], [ 2, 1 ] ], width:2 }
m.exp(0) ==> { height:2, mtx:[ [ 1, 0 ], [ 0, 1 ] ], width:2 }
m.exp(1) ==> { height:2, mtx:[ [ 3, 2 ], [ 2, 1 ] ], width:2 }
m.exp(2) ==> { height:2, mtx:[ [ 13, 8 ], [ 8, 5 ] ], width:2 }
m.exp(4) ==> { height:2, mtx:[ [ 233, 144 ], [ 144, 89 ] ], width:2 }
m.exp(10) ==> { height:2, mtx:[ [ 1346269, 832040 ], [ 832040, 514229 ] ], width:2 }
=!EXPECTEND!=
*/

Output:

prompt$ jsish -u matrixExponentiation.jsi
[PASS] matrixExponentiation.jsi

Julia

Matrix exponentiation is implemented by the built-in ^ operator.

julia> [1 1 ; 1 0]^10
2x2 Array{Int64,2}:
 89  55
 55  34

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:

Output:

K Console - Enter \ for help

  \l mpow

  a:(3 2;2 1)
(3 2
 2 1)
  pow[a;0]
(1 0
 0 1)
  pow[a;1]
(3 2
 2 1)
  pow[a;2]
(13 8
 8 5)
  pow[a;3]
(55 34
 34 21)
  pow[a;4]
(233 144
 144 89)
  pow[a;10]
(1346269 832040
 832040 514229)

Kotlin

// version 1.1.3

typealias Vector = DoubleArray
typealias Matrix = Array<Vector>

operator fun Matrix.times(other: Matrix): Matrix {
    val rows1 = this.size
    val cols1 = this[0].size
    val rows2 = other.size
    val cols2 = other[0].size
    require(cols1 == rows2)
    val result = Matrix(rows1) { Vector(cols2) }
    for (i in 0 until rows1) {
        for (j in 0 until cols2) {
            for (k in 0 until rows2) {
                result[i][j] += this[i][k] * other[k][j]
            }
        }
    }
    return result
}

fun identityMatrix(n: Int): Matrix {
    require(n >= 1) 
    val ident = Matrix(n) { Vector(n) }
    for (i in 0 until n) ident[i][i] = 1.0
    return ident
}

infix fun Matrix.pow(n : Int): Matrix {
    require (n >= 0 && this.size == this[0].size)
    if (n == 0) return identityMatrix(this.size)
    if (n == 1) return this
    var pow = identityMatrix(this.size)
    var base = this
    var e = n
    while (e > 0) {
        if ((e and 1) == 1) pow *= base
        e = e shr 1
        base *= base
    }
    return pow
}  

fun printMatrix(m: Matrix, n: Int) {
    println("** Power of $n **")
    for (i in 0 until m.size) println(m[i].contentToString())
    println()
}

fun main(args: Array<String>) {
    val m = arrayOf(
        doubleArrayOf(3.0, 2.0),
        doubleArrayOf(2.0, 1.0)
    )
    for (i in 0..10) printMatrix(m pow i, i)
}

Output:

** Power of 0 **
[1.0, 0.0]
[0.0, 1.0]

** Power of 1 **
[3.0, 2.0]
[2.0, 1.0]

** Power of 2 **
[13.0, 8.0]
[8.0, 5.0]

** Power of 3 **
[55.0, 34.0]
[34.0, 21.0]

** Power of 4 **
[233.0, 144.0]
[144.0, 89.0]

** Power of 5 **
[987.0, 610.0]
[610.0, 377.0]

** Power of 6 **
[4181.0, 2584.0]
[2584.0, 1597.0]

** Power of 7 **
[17711.0, 10946.0]
[10946.0, 6765.0]

** Power of 8 **
[75025.0, 46368.0]
[46368.0, 28657.0]

** Power of 9 **
[317811.0, 196418.0]
[196418.0, 121393.0]

** Power of 10 **
[1346269.0, 832040.0]
[832040.0, 514229.0]

Lambdatalk

{require lib_matrix}

{def M.exp
 {lambda {:m :n}
  {if {= :n 0}
   then {M.new [ [1,0],[0,1] ]}
   else {S.reduce M.multiply {S.map {{lambda {:m _} :m} :m} {S.serie 1 :n}}}}}}
-> M.exp

'{def M 
 {M.new [[3,2],
         [2,1]]}}
-> M

{S.map {lambda {:i} {br}M{sup :i} = {M.exp {M} :i}} 
       0 1 2 3 4 10} 
->
M^0 = [[1,0],[0,1]] 
M^1 = [[3,2],[2,1]] 
M^2 = [[13,8],[8,5]] 
M^3 = [[55,34],[34,21]] 
M^4 = [[233,144],[144,89]] 
M^10 = [[1346269,832040],[832040,514229]]

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"


print "Exponentiation of a matrix"
call DisplayMatrix MatrixD$
print "         Raised to power 5 ="
MatrixE$ =MatrixToPower$( MatrixD$, 5)
call DisplayMatrix MatrixE$
print "         Raised to power 9 ="
MatrixE$ =MatrixToPower$( MatrixD$, 9)
call DisplayMatrix MatrixE$

Output:

Exponentiation of a matrix
| 0.86603 0.50000 0.00000 |
| -0.50000 0.86603 0.00000 |
| 0.00000 0.00000 1.00000 |

Raised to power 5 =
| -0.86604 0.50002 0.00000 |
| -0.50002 -0.86604 0.00000 |
| 0.00000 0.00000 1.00000 |

Raised to power 9 =
| -0.00002 -1.00004 0.00000 |
| 1.00004 -0.00002 0.00000 |
| 0.00000 0.00000 1.00000 |

Lua

Matrix = {}

function Matrix.new( dim_y, dim_x )
    assert( dim_y and dim_x )
    
    local matrix = {}
    local metatab = {}
    setmetatable( matrix, metatab )
    metatab.__add = Matrix.Add
    metatab.__mul = Matrix.Mul
    metatab.__pow = Matrix.Pow

    matrix.dim_y = dim_y
    matrix.dim_x = dim_x 
    
    matrix.data = {}
    for i = 1, dim_y do
        matrix.data[i] = {}
    end
    return matrix
end

function Matrix.Show( m )
    for i = 1, m.dim_y do
        for j = 1, m.dim_x do
            io.write( tostring( m.data[i][j] ), " " )
        end
        io.write( "\n" )
    end
end

function Matrix.Add( m, n )
    assert( m.dim_x == n.dim_x and m.dim_y == n.dim_y )
 
    local r = Matrix.new( m.dim_y, m.dim_x )
    for i = 1, m.dim_y do
        for j = 1, m.dim_x do
            r.data[i][j] = m.data[i][j] + n.data[i][j]
        end
    end
    return r
end

function Matrix.Mul( m, n )
    assert( m.dim_x == n.dim_y )
  
    local r = Matrix.new( m.dim_y, n.dim_x )
    for i = 1, m.dim_y do
        for j = 1, n.dim_x do
            r.data[i][j] = 0
            for k = 1, m.dim_x do
                r.data[i][j] = r.data[i][j] + m.data[i][k] * n.data[k][j]
            end
        end
    end
    return r
end

function Matrix.Pow( m, p )
    assert( m.dim_x == m.dim_y )
    
    local r = Matrix.new( m.dim_y, m.dim_x )
    
    if p == 0 then 
        for i = 1, m.dim_y do
            for j = 1, m.dim_x do
                if i == j then
                    r.data[i][j] = 1
                else
                    r.data[i][j] = 0
                end
            end
        end
    elseif p == 1 then
        for i = 1, m.dim_y do
            for j = 1, m.dim_x do
                r.data[i][j] = m.data[i][j]
            end
        end        
    else
        r = m
        for i = 2, p do
            r = r * m
        end
    end
    
    return r
end


m = Matrix.new( 2, 2 )
m.data = { { 1, 2 }, { 3, 4 } }

n = m^4;

Matrix.Show( n )

M2000 Interpreter

Module CheckIt {
	Class cArray {
		a=(,)
		Function Power(n as integer){
			cArr=This     ' create a copy
			dim new()
			new()=cArr.a   ' get a pointer from a to new()
			Let cArr.a=new()    ' now new() return a copy
			cArr.a*=0  ' make zero all elements
			link cArr.a to v()
			for i=dimension(cArr.a,1,0) to dimension(cArr.a, 1,1) : v(i,i)=1: next i
			while n>0
				let cArr=cArr*this    ' * is the operator "*"
				n--
			end while
			=cArr
		}
		Operator "*"{
			Read cArr
			b=cArr.a
			if dimension(.a)<>2 or dimension(b)<>2 then Error "Need two 2D arrays "
			let a2=dimension(.a,2), b1=dimension(b,1)
			if a2<>b1 then Error "Need columns of first array equal to rows of second array"
			let a1=dimension(.a,1), b2=dimension(b,2)
			let aBase=dimension(.a,1,0)-1, bBase=dimension(b,1,0)-1
			let aBase1=dimension(.a,2,0)-1, bBase1=dimension(b,2,0)-1
			link .a,b to a(), b()  ' change interface for arrays
			dim base 1, c(a1, b2)
			for i=1 to a1 : let ia=i+abase : for j=1 to b2 : let jb=j+bBase1 : for k=1 to a2
			c(i,j)+=a(ia,k+aBase1)*b(k+bBase,jb)
			next k : next j : next i
			\\ redim to base 0
			dim base 0, c(a1, b2)
			.a<=c()
			}
		Module Print {
			link .a to v()
			for i=dimension(.a,1,0) to dimension(.a, 1,1) 
			for j=dimension(.a,2,0) to dimension(.a, 2,1) 
			print  v(i,j),: next j: print : next i
				
		}
	Class:
		\\ this module used as constructor, and not returned to final group (user object in M2000)
		Module cArray (r) {
			c=r
			Dim a(r,c)
			For i=0 to r-1 : For j=0 to c-1: Read a(i,j): Next j : Next i
			.a<=a()
		}
	}
	Print "matrix():"
	P=cArray(2,3,2,2,1)
	P.Print
	For i=0 to 9 
		Print "matrix()^"+str$(i,0)+"="
		K=P.Power(i)
		K.Print
	next i
}
Checkit

Output:

matrix():
      3      2
      2      1
matrix()^0=
      1      0
      0      1
matrix()^1=
      3      2
      2      1
matrix()^2=
     13      8
      8      5
matrix()^3=
     55     34
     34     21
matrix()^4=
    233    144
    144     89
matrix()^5=
    987    610
    610    377
matrix()^6=
   4181   2584
   2584   1597
matrix()^7=
  17711  10946
  10946   6765
matrix()^8=
  75025  46368
  46368  28657
matrix()^9=
 317811 196418
 196418 121393

Maple

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

> M := <<1,2>|<3,4>>;
> M ^ 2;

$\left[{\begin{array}{cc}7&15\\10&22\end{array}}\right]$

If you want elementwise powers, you can use the elementwise ^~ operator:

> M := <<1,2>|<3,4>>;
> M ^~ 2;

$\left[{\begin{array}{cc}1&9\\4&16\end{array}}\right]$

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

a = {{3, 2}, {4, 1}};
MatrixPower[a, 0]
MatrixPower[a, 1]
MatrixPower[a, -1]
MatrixPower[a, 4]
MatrixPower[a, 1/2]
MatrixPower[a, Pi]

gives back:

$\left({\begin{array}{cc}1&0\\0&1\end{array}}\right)$

$\left({\begin{array}{cc}3&2\\4&1\end{array}}\right)$

$\left({\begin{array}{cc}-{\frac {1}{5}}&{\frac {2}{5}}\\{\frac {4}{5}}&-{\frac {3}{5}}\end{array}}\right)$

$\left({\begin{array}{cc}417&208\\416&209\end{array}}\right)$

$\left({\begin{array}{cc}{\frac {2{\sqrt {5}}}{3}}+{\frac {i}{3}}&{\frac {\sqrt {5}}{3}}-{\frac {i}{3}}\\{\frac {2{\sqrt {5}}}{3}}-{\frac {2i}{3}}&{\frac {\sqrt {5}}{3}}+{\frac {2i}{3}}\end{array}}\right)$

$\left({\begin{array}{cc}{\frac {(-1)^{\pi }}{3}}+2{\frac {5^{\pi }}{3}}&{\frac {5^{\pi }}{3}}-{\frac {1}{3}}(-1)^{\pi }\\2{\frac {5^{\pi }}{3}}-{\frac {2}{3}}(-1)^{\pi }&{\frac {2(-1)^{\pi }}{3}}+{\frac {5^{\pi }}{3}}\end{array}}\right)$

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

MatrixPower[{{i, j}, {k, l}}, m] // Simplify

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

$\left({\begin{array}{cc}{\frac {2^{-m-1}\left(\left({\sqrt {i^{2}-2il+4jk+l^{2}}}-i+l\right)\left(-{\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}+\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i-l\right)\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}\right)}{\sqrt {i^{2}-2il+4jk+l^{2}}}}&{\frac {j2^{-m}\left(\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}-\left(-{\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}\right)}{\sqrt {i^{2}-2il+4jk+l^{2}}}}\\{\frac {k2^{-m}\left(\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}-\left(-{\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}\right)}{\sqrt {i^{2}-2il+4jk+l^{2}}}}&{\frac {2^{-m-1}\left(\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i-l\right)\left(-{\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}+\left({\sqrt {i^{2}-2il+4jk+l^{2}}}-i+l\right)\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}\right)}{\sqrt {i^{2}-2il+4jk+l^{2}}}}\end{array}}\right)$

Final note: Do not confuse MatrixPower with MatrixExp; the former is for matrix exponentiation, and the latter for the matrix exponential (E^m).

MATLAB

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

function [output] = matrixexponentiation(matrixA, exponent)
   output = matrixA^(exponent);

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

function [output] = matrixexponentiation(matrixA, exponent)
   output = matrixA.^(exponent);

Maxima

a: matrix([3, 2],
          [4, 1])$

a ^^ 4;
/* matrix([417, 208],
          [416, 209]) */

a ^^ -1;
/* matrix([-1/5, 2/5],
          [4/5, -3/5]) */

Nim

import sequtils, strutils

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

func `*`[N, T](a, b: Matrix[N, T]): Matrix[N, T] =
  for i in 1..N:
    for j in 1..N:
      for k in 1..N:
        result[i][j] += a[i][k] * b[k][j]


func identityMatrix[N; T](): Matrix[N, T] =
  for i in 1..N:
    result[i][i] = T(1)


func `^`[N, T](m: Matrix[N, T]; n: Natural): Matrix[N, T] =
  if n == 0: return identityMatrix[N, T]()
  if n == 1: return m
  var n = n
  var m = m
  result = identityMatrix[N, T]()
  while n > 0:
    if (n and 1) != 0:
      result = result * m
    n = n shr 1
    m = m * m


proc `$`(m: Matrix): string =
  var lg = 0
  for i in 1..m.N:
    for j in 1..m.N:
      lg = max(lg, len($m[i][j]))
  for i in 1..m.N:
    echo m[i].mapIt(align($it, lg)).join(" ")


when isMainModule:

  let m1: Matrix[3, int] = [[ 3, 2, -1],
                            [-1, 0,  5],
                            [ 2, -1, 3]]
  echo m1^10

  import math
  const
    C30 = sqrt(3.0) / 2
    S30 = 1 / 2
  let m2: Matrix[2, float] = [[C30, -S30], [S30,  C30]]  # 30° rotation matrix.
  echo m2^12    # Nearly the identity matrix.

Output:

572880 154352 321344
480752 261648 306176
473168 161936 413472

    0.9999999999999993 -3.885780586188048e-16
 3.885780586188048e-16     0.9999999999999993

OCaml

We will use some auxiliary functions

(* identity matrix *)
let eye n =
  let a = Array.make_matrix n n 0.0 in
  for i=0 to n-1 do
    a.(i).(i) <- 1.0
  done;
  (a)
;;

(* matrix dimensions *)
let dim a = Array.length a, Array.length a.(0);;

(* make matrix from list in row-major order *)
let matrix p q v =
  if (List.length v) <> (p * q)
  then failwith "bad dimensions"
  else
    let a = Array.make_matrix p q (List.hd v) in
    let rec g i j = function
    | [] -> a
    | x::v ->
        a.(i).(j) <- x;
        if j+1 < q
        then g i (j+1) v
        else g (i+1) 0 v
    in
    g 0 0 v
;;

(* matrix product *)
let matmul a b =
  let n, p = dim a
  and q, r = dim b in
  if p <> q then failwith "bad dimensions" else
  let c = Array.make_matrix n r 0.0 in
  for i=0 to n-1 do
    for j=0 to r-1 do
      for k=0 to p-1 do
        c.(i).(j) <- c.(i).(j) +. a.(i).(k) *. b.(k).(j)
      done
    done
  done;
  (c)
;;

(* generic exponentiation, usual algorithm *)
let pow one mul a n =
  let rec g p x = function
  | 0 -> x
  | i ->
      g (mul p p) (if i mod 2 = 1 then mul p x else x) (i/2)
  in
  g a one n
;;

(* example with integers *)
pow 1 ( * ) 2 16;;
(* - : int = 65536 *)

Now matrix power is simply a special case of pow :

let matpow a n =
  let p, q = dim a in
  if p <> q then failwith "bad dimensions" else
  pow (eye p) matmul a n;;

matpow (matrix 2 2 [ 1.0; 1.0; 1.0; 0.0 ]) 10;;
(* - : float array array = [|[|89.; 55.|]; [|55.; 34.|]|] *)

(* use as infix operator *)
let ( ^^ ) = matpow;;

[| [| 1.0; 1.0|]; [| 1.0; 0.0 |] |] ^^ 10;;
(* - : float array array = [|[|89.; 55.|]; [|55.; 34.|]|] *)

Octave

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

M = [ 3, 2; 2, 1 ];
M^0
M^1
M^2
M^(-1)
M^0.5

Output:

ans =

   1   0
   0   1

ans =

   3   2
   2   1

ans =

   13    8
    8    5

ans =

  -1.0000   2.0000
   2.0000  -3.0000

ans =

   1.48931 + 0.13429i   0.92044 - 0.21729i
   0.92044 - 0.21729i   0.56886 + 0.35158i

(Of course this is not an implementation, but it can be used as reference for the results)

PARI/GP

M^n

Perl

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

use overload (
        '""'    => \&_string,   # overload string operator so we can just print
        '*'     => \&_mult,     # multiplication, needed for expo
        '*='    => \&_mult,     # ditto, explicitly defined to trigger copy
        '**'    => \&_expo,     # overload exponentiation
        '='     => \&_copy,     # copy operator
);

sub make {
        my $cls = shift;
        my $n = @_;
        for (@_) {
                # verify each row given is the right length
                confess "Bad data @$_: matrix must be square "
                        if @$_ != $n;
        }

        bless [ map [@$_], @_ ] # important: actually copy all the rows
}

sub identity {
        my $self = shift;
        my $n = @$self - 1;
        my @rows = map [ (0) x $_, 1, (0) x ($n - $_) ], 0 .. $n;
        bless \@rows
}

sub zero {
        my $self = shift;
        my $n = @$self;
        bless [ map [ (0) x $n ], 1 .. $n ]
}

sub _string {
        "[ ".join("\n  " =>
                map join(" " => map(sprintf("%12.6g", $_), @$_)), @{+shift}
        )."  ]\n";
}

sub _mult {
        my ($a, $b) = @_;
        my $x = $a->zero;
        my @idx = (0 .. $#$x);
        for my $j (@idx) {
                my @col = map($a->[$_][$j], @idx);
                for my $i (@idx) {
                        my $row = $b->[$i];
                        $x->[$i][$j] += $row->[$_] * $col[$_] for @idx;
                }
        }
        $x
}

sub _expo {
        my ($self, $n) = @_;
        confess "matrix **: must be non-negative integer power"
                        unless $n >= 0 && $n == int($n);

        my ($tmp, $out) = ($self, $self->identity);
        do {
                $out *= $tmp    if $n & 1;
                $tmp *= $tmp;
        } while $n >>= 1;

        $out
}

sub _copy { bless [ map [ @$_ ], @{+shift} ] }

# now use our matrix class
package main;

my $m = SquareMatrix->make(
                [1, 2, 0],
                [0, 3, 1],
                [1, 0, 0] );
print "### Order $_\n", $m ** $_        for 0 .. 10;

$m = SquareMatrix->make(
        [ 1.0001, 0,      0, 1       ],
        [ 0,      1.001,  0, 0       ],
        [ 0,      0,      1, 0.99998 ],
        [ 1e-8,   0,      0, 1.0002  ]);

print "\n### Matrix is now\n",  $m;
print "\n### Big power:\n",     $m ** 100_000;
print "\n### Too big:\n",       $m ** 1_000_000;
print "\n### WAY too big:\n",   $m ** 1_000_000_000_000;
print "\n### But identity matrix can handle that\n",
                $m->identity ** 1_000_000_000_000;

Phix

Phix does not permit operator overloading, however here is a simple function to raise a square matrix to a non-negative integer power.
First two routines copied straight from the Identity_matrix and Matrix_multiplication tasks.

with javascript_semantics
function identity(integer n)
    sequence res = repeat(repeat(0,n),n)
    for i=1 to n do
        res[i][i] = 1
    end for
    return res
end function
 
function matrix_mul(sequence a, b)
    integer {ha,wa,hb,wb} = apply({a,a[1],b,b[1]},length)
    if wa!=hb then return 0 end if
    sequence c = repeat(repeat(0,wb),ha)
    for i=1 to ha do
        for j=1 to wb do
            for k=1 to wa do
                c[i][j] += a[i][k]*b[k][j]
            end for
        end for
    end for
    return c
end function
 
function matrix_exponent(sequence m, integer n)
    integer l = length(m)
    if n=0 then return identity(l) end if
    sequence res = m
    for i=2 to n do
        res = matrix_mul(res,m)
    end for
    return res
end function
 
constant M1 = {{5}},
         M2 = {{3, 2},
               {2, 1}},
         M3 = {{1, 2, 0},
               {0, 3, 1},
               {1, 0, 0}}
 
ppOpt({pp_Nest,1})
pp(matrix_exponent(M1,0))
pp(matrix_exponent(M1,1))
pp(matrix_exponent(M1,2))
puts(1,"==\n")
pp(matrix_exponent(M2,0))
pp(matrix_exponent(M2,1))
pp(matrix_exponent(M2,2))
pp(matrix_exponent(M2,10))
puts(1,"==\n")
pp(matrix_exponent(M3,10))
puts(1,"==\n")
pp(matrix_exponent(identity(4),5))

Output:

{{1}}
{{5}}
{{25}}
==
{{1,0},
 {0,1}}
{{3,2},
 {2,1}}
{{13,8},
 {8,5}}
{{1346269,832040},
 {832040,514229}}
==
{{13801,102408,31322},
 {15661,116209,35543},
 {4221,31322,9580}}
==
{{1,0,0,0},
 {0,1,0,0},
 {0,0,1,0},
 {0,0,0,1}}

PicoLisp

Uses the 'matMul' function from Matrix multiplication#PicoLisp

(de matIdent (N)
   (let L (need N (1) 0)
      (mapcar '(() (copy (rot L))) L) ) )

(de matExp (Mat N)
   (let M (matIdent (length Mat))
      (do N
         (setq M (matMul M Mat)) )
      M ) )

(matExp '((3 2) (2 1)) 3)

Output:

-> ((55 34) (34 21))

Python

Using matrixMul from Matrix multiplication#Python

>>> from operator import mul
>>> def matrixMul(m1, m2):
  return map(
    lambda row:
      map(
        lambda *column:
          sum(map(mul, row, column)),
        *m2),
    m1)

>>> def identity(size):
	size = range(size)
	return [[(i==j)*1 for i in size] for j in size]

>>> def matrixExp(m, pow):
	assert pow>=0 and int(pow)==pow, "Only non-negative, integer powers allowed"
	accumulator = identity(len(m))
	for i in range(pow):
		accumulator = matrixMul(accumulator, m)
	return accumulator

>>> def printtable(data):
	for row in data:
		print ' '.join('%-5s' % ('%s' % cell) for cell in row)

		
>>> m = [[3,2], [2,1]]
>>> for i in range(5):
	print '\n%i:' % i
	printtable( matrixExp(m, i) )

	

0:
1     0    
0     1    

1:
3     2    
2     1    

2:
13    8    
8     5    

3:
55    34   
34    21   

4:
233   144  
144   89   
>>> printtable( matrixExp(m, 10) )
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) :
        A = self
        return Mat([[sum(A[i][k]*B[k][j] for k in range(len(B)))
                    for j in range(len(B[0])) ] for i in range(len(A))])
    

def identity(size):
    size = range(size)
    return [[(i==j)*1 for i in size] for j in size]

def power(F, n):
    result = Mat(identity(len(F)))
    b = Mat(F)
    while n > 0:
        if (n%2) == 0:
            b = b @ b
            n //= 2
        else:
            result = b @ result
            b = b @ b
            n //= 2
    return result

def printtable(data):
    for row in data:
        print (' '.join('%-5s' % ('%s' % cell) for cell in row))

m = [[3,2], [2,1]]
for i in range(5):
    print('\n%i:' % i)
    printtable(power(m, i))

Output:

0:
[[1, 0], [0, 1]]

1:
[[3, 2], [2, 1]]

2:
[[13, 8], [8, 5]]

3:
[[55, 34], [34, 21]]

4:
[[233, 144], [144, 89]]

R

Library function call

Library: Biodem

library(Biodem)
m <- matrix(c(3,2,2,1), nrow=2)
mtx.exp(m, 0)
#      [,1] [,2]
# [1,]    1    0
# [2,]    0    1
mtx.exp(m, 1)
#      [,1] [,2]
# [1,]    3    2
# [2,]    2    1
mtx.exp(m, 2)
#      [,1] [,2]
# [1,]   13    8
# [2,]    8    5
mtx.exp(m, 3)
#      [,1] [,2]
# [1,]   55   34
# [2,]   34   21
mtx.exp(m, 10)
#         [,1]   [,2]
# [1,] 1346269 832040
# [2,]  832040 514229

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 <- matrix(c(1, 2, 3, 4), 2, 2)
a^1
a^2

Output:

> a^1
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> a^2
     [,1] [,2]
[1,]    1    9
[2,]    4   16

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?

`%^%` <- function(mat, n)
{
  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(n==0) diag(nrow = nrow(mat))#Identity matrix of mat's dimensions
    else if(n == 1) mat
    else if(n > 1) mat %*% (mat %^% (n - 1))
    else stop("Invalid n.")
  }
  else stop("Invalid input type.")
}
#For output:
a %^% 0
a %^% 1
a %^% 2
a %*% a %*% a#Base R's equivalent of a %^% 3
a %^% 3
nonSquareMatrix <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3)
nonSquareMatrix %^% 1
nonSquareMatrix %^% 2#R's %*% will throw the error for us

Output:

> a %^% 0
     [,1] [,2]
[1,]    1    0
[2,]    0    1

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

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

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

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

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

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

> nonSquareMatrix %^% 2#R's %*% will throw the error for us
Error in mat %*% (mat %^% (n - 1)) : non-conformable arguments

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:

library(Biodem)
`%^%` <- function(mat, n) Biodem::mtx.exp(mat, n)

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

Racket

#lang racket
(require math)

(define a (matrix ((3 2) (2 1))))

;; Using the builtin matrix exponentiation
(for ([i 11])
  (printf "a^~a = ~s\n" i (matrix-expt a i)))

;; Output:
;; a^0 = (array #[#[1 0] #[0 1]])
;; a^1 = (array #[#[3 2] #[2 1]])
;; a^2 = (array #[#[13 8] #[8 5]])
;; a^3 = (array #[#[55 34] #[34 21]])
;; a^4 = (array #[#[233 144] #[144 89]])
;; a^5 = (array #[#[987 610] #[610 377]])
;; a^6 = (array #[#[4181 2584] #[2584 1597]])
;; a^7 = (array #[#[17711 10946] #[10946 6765]])
;; a^8 = (array #[#[75025 46368] #[46368 28657]])
;; a^9 = (array #[#[317811 196418] #[196418 121393]])
;; a^10 = (array #[#[1346269 832040] #[832040 514229]])

;; But it could be implemented manually, using matrix multiplication
(define (mpower M p)
  (cond [(= p 1) M]
        [(even? p) (mpower (matrix* M M) (/ p 2))]
        [else (matrix* M (mpower M (sub1 p)))]))
(for ([i (in-range 1 11)])
  (printf "a^~a = ~s\n" i (matrix-expt a i)))

Raku

(formerly Perl 6)

subset SqMat of Array where { .elems == all(.[]».elems) }

multi infix:<*>(SqMat $a, SqMat $b) {[
    for ^$a -> $r {[
        for ^$b[0] -> $c {
            [+] ($a[$r][] Z* $b[].map: *[$c])
        }
    ]}
]}

multi infix:<**> (SqMat $m, Int $n is copy where { $_ >= 0 }) {
    my $tmp = $m;
    my $out = [for ^$m -> $i { [ for ^$m -> $j { +($i == $j) } ] } ];
    loop {
        $out = $out * $tmp if $n +& 1;
        last unless $n +>= 1;
        $tmp = $tmp * $tmp;
    }

    $out;
}

multi show (SqMat $m) {
    my $size = $m.map( *.list».chars ).flat.max;
    say .fmt("%{$size}s", ' ') for $m.list;
}

my @m = [1, 2, 0],
        [0, 3, 1],
        [1, 0, 0];

for 0 .. 10 -> $order {
    say "### Order $order";
    show @m ** $order;
}

Output:

### Order 0
1 0 0
0 1 0
0 0 1
### Order 1
1 2 0
0 3 1
1 0 0
### Order 2
1 8 2
1 9 3
1 2 0
### Order 3
 3 26  8
 4 29  9
 1  8  2
### Order 4
11 84 26
13 95 29
 3 26  8
### Order 5
 37 274  84
 42 311  95
 11  84  26
### Order 6
 121  896  274
 137 1017  311
  37  274   84
### Order 7
 395 2930  896
 448 3325 1017
 121  896  274
### Order 8
 1291  9580  2930
 1465 10871  3325
  395  2930   896
### Order 9
 4221 31322  9580
 4790 35543 10871
 1291  9580  2930
### Order 10
 13801 102408  31322
 15661 116209  35543
  4221  31322   9580

RPL

Operators can not be overloaded, but we can easily create a new word, with same syntax as the classical exponentiation operator. the power must be a signed integer.

RPL code	Comment
≪ SWAP IF OVER 0 < THEN INV END DUP IDN → m id ≪ ABS id WHILE OVER REPEAT m * SWAP 1 - SWAP END SWAP DROP ≫ ≫ 'MATXP' STO	MATXP ( m n -- m^n ) inverse matrix if n<0 store matrix and identity initialize stack with abs(n) and identity multiply n times clean stack return m^n

[[3 2][2 1]] 0 MATXP
[[3 2][2 1]] 1 MATXP
[[3 2][2 1]] 2 MATXP
[[3 2][2 1]] 5 MATXP
[[3 2][2 1]] -5 MATXP
{{out}

5:            [[ 1 0 ]
              [ 0 1 ]]
4:            [[ 3 2 ]
              [ 2 1 ]]
3:           [[ 13 8 ]
              [ 8 5 ]]
2:        [[ 987 610 ]
          [ 610 377 ]]
1:       [[ -377 610 ]
          [ 610 -987]]

Ruby

Ruby's standard library already provides the matrix-exponentiation operator. It is Matrix#** from package 'matrix' of the standard library. MRI 1.9.x implements the matrix-exponentiation operator in file matrix.rb, def ** (around line 961).

$ irb
irb(main):001:0> require 'matrix'
=> true
irb(main):002:0> m=Matrix[[3,2],[2,1]]
=> Matrix[[3, 2], [2, 1]]
irb(main):003:0> m**0
=> Matrix[[1, 0], [0, 1]]
irb(main):004:0> m ** 1
=> Matrix[[3, 2], [2, 1]]
irb(main):005:0> m ** 2
=> Matrix[[13, 8], [8, 5]]
irb(main):006:0> m ** 5
=> Matrix[[987, 610], [610, 377]]
irb(main):007:0> m ** 10
=> Matrix[[1346269, 832040], [832040, 514229]]

Starting with Ruby 1.9.3, it can also calculate Matrix ** Float.

Works with: Ruby version 1.9.3

irb(main):008:0> m ** 1.5
=> Matrix[[(6.308803769316981-0.03170173099577213i), (3.8990551577913446+0.05129
4478253365354i)], [(3.899055157791345+0.05129447825336536i), (2.4097486115256355
-0.0829962092491375i)]]

Rust

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

use std::fmt;
use std::ops;
const WIDTH: usize = 6;

#[derive(Clone)]
struct SqMat {
    data: Vec<Vec<i64>>,
}

impl fmt::Debug for SqMat {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let mut row = "".to_string();
        for i in &self.data {
            for j in i {
                row += &format!("{:>w$} ", j, w = WIDTH);
            }
            row += &"\n";
        }
        write!(f, "{}", row)
    }
}

impl ops::BitXor<u32> for SqMat {
    type Output = Self;

    fn bitxor(self, n: u32) -> Self::Output {
        let mut aux = self.data.clone();
        let mut ans: SqMat = SqMat {
            data: vec![vec![0; aux.len()]; aux.len()],
        };
        for i in 0..aux.len() {
            ans.data[i][i] = 1;
        }
        let mut b = n;
        while b > 0 {
            if b & 1 > 0 {
                // ans = ans * aux
                let mut tmp = aux.clone();
                for i in 0..aux.len() {
                    for j in 0..aux.len() {
                        tmp[i][j] = 0;
                        for k in 0..aux.len() {
                            tmp[i][j] += ans.data[i][k] * aux[k][j];
                        }
                    }
                }
                ans.data = tmp;
            }
            b >>= 1;
            if b > 0 {
                // aux = aux * aux
                let mut tmp = aux.clone();
                for i in 0..aux.len() {
                    for j in 0..aux.len() {
                        tmp[i][j] = 0;
                        for k in 0..aux.len() {
                            tmp[i][j] += aux[i][k] * aux[k][j];
                        }
                    }
                }
                aux = tmp;
            }
        }
        ans
    }
}

fn main() {
    let sm: SqMat = SqMat {
        data: vec![vec![1, 2, 0], vec![0, 3, 1], vec![1, 0, 0]],
    };
    for i in 0..11 {
        println!("Power of {}:\n{:?}", i, sm.clone() ^ i);
    }
}

Output:

Power of 0:
     1      0      0
     0      1      0
     0      0      1

Power of 1:
     1      2      0
     0      3      1
     1      0      0

Power of 2:
     1      8      2
     1      9      3
     1      2      0

Power of 3:
     3     26      8
     4     29      9
     1      8      2

Power of 4:
    11     84     26
    13     95     29
     3     26      8

Power of 5:
    37    274     84
    42    311     95
    11     84     26

Power of 6:
   121    896    274
   137   1017    311
    37    274     84

Power of 7:
   395   2930    896
   448   3325   1017
   121    896    274

Power of 8:
  1291   9580   2930
  1465  10871   3325
   395   2930    896

Power of 9:
  4221  31322   9580
  4790  35543  10871
  1291   9580   2930

Power of 10:
 13801 102408  31322
 15661 116209  35543
  4221  31322   9580

Scala

class Matrix[T](matrix:Array[Array[T]])(implicit n: Numeric[T], m: ClassManifest[T])
{
  import n._
  val rows=matrix.size
  val cols=matrix(0).size
  def row(i:Int)=matrix(i)
  def col(i:Int)=matrix map (_(i))

  def *(other: Matrix[T]):Matrix[T] = new Matrix(
    Array.tabulate(rows, other.cols)((row, col) =>
      (this.row(row), other.col(col)).zipped.map(_*_) reduceLeft (_+_)
  ))

  def **(x: Int)=x match {
    case 0 => createIdentityMatrix
    case 1 => this
    case 2 => this * this
    case _ => List.fill(x)(this) reduceLeft (_*_)
  }
	
  def createIdentityMatrix=new Matrix(Array.tabulate(rows, cols)((row,col) => 
    if (row == col) one else zero)
  )

  override def toString = matrix map (_.mkString("[", ", ", "]")) mkString "\n"
}

object MatrixTest {
  def main(args:Array[String])={
    val m=new Matrix[BigInt](Array(Array(3,2), Array(2,1)))
    println("-- m --\n"+m)
				
    Seq(0,1,2,3,4,10,20,50) foreach {x =>
      println("-- m**"+x+" --")
      println(m**x)
    }
  }
}

Output:

-- m --
[3, 2]
[2, 1]
-- m**0 --
[1, 0]
[0, 1]
-- m**1 --
[3, 2]
[2, 1]
-- m**2 --
[13, 8]
[8, 5]
-- m**3 --
[55, 34]
[34, 21]
-- m**4 --
[233, 144]
[144, 89]
-- m**10 --
[1346269, 832040]
[832040, 514229]
-- m**20 --
[2504730781961, 1548008755920]
[1548008755920, 956722026041]
-- m**50 --
[16130531424904581415797907386349, 9969216677189303386214405760200]
[9969216677189303386214405760200, 6161314747715278029583501626149]

Scheme

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 (halve x)
  (/ x 2))

(define (row*col row col)
  (apply + (map * row col)))

(define (matrix-multiply m1 m2)
  (map 
    (lambda (row) 
      (apply map (lambda col (row*col row col)) 
        m2)) 
    m1))

(define (matrix-exp mat exp)
  (cond ((= exp 1) mat)
        ((even? exp) (square-matrix (matrix-exp mat (halve exp))))
        (else (matrix-multiply mat (matrix-exp mat (dec exp))))))

(define (square-matrix mat)
  (matrix-multiply mat mat))

Output:

> (matrix-exp '((3 2) (2 1)) 50)
((16130531424904581415797907386349 9969216677189303386214405760200)
 (9969216677189303386214405760200 6161314747715278029583501626149))

Seed7

The example below uses several features of Seed7:

• Overloading of the operators * and ** .
• The template enable_output, which allows writing a matrix with write (the function str must be defined before calling enable_output).
• A for loop which loops over values listed in an array literal

$ include "seed7_05.s7i";
  include "float.s7i";

const type: matrix is array array float;
 
const func string: str (in matrix: mat) is func
  result
    var string: stri is "";
  local
    var integer: row is 0;
    var integer: column is 0;
  begin
    for row range 1 to length(mat) do
      for column range 1 to length(mat[row]) do
        stri &:= str(mat[row][column]);
        if column < length(mat[row]) then
          stri &:= ", ";
        end if;
      end for;
      if row < length(mat) then
        stri &:= "\n";
      end if;
    end for;
  end func;

enable_output(matrix);

const func matrix: (in matrix: mat1) * (in matrix: mat2) is func
  result
    var matrix: product is matrix.value;
  local
    var integer: row is 0;
    var integer: column is 0;
    var integer: k is 0;
  begin
    product := length(mat1) times length(mat1) times 0.0;
    for row range 1 to length(mat1) do
      for column range 1 to length(mat1) do
        product[row][column] := 0.0;
        for k range 1 to length(mat1) do
          product[row][column] +:= mat1[row][k] * mat2[k][column];
        end for;
      end for;
    end for;
  end func;

const func matrix: (in var matrix: base) ** (in var integer: exponent) is func
  result
    var matrix: power is matrix.value;
  local
    var integer: row is 0;
    var integer: column is 0;
  begin
    if exponent < 0 then
      raise NUMERIC_ERROR;
    else
      if odd(exponent) then
        power := base;
      else
        # Create identity matrix
        power := length(base) times length(base) times 0.0;
        for row range 1 to length(base) do
          for column range 1 to length(base) do
            if row = column then
              power[row][column] := 1.0;
            end if;
          end for;
        end for;
      end if;
      exponent := exponent div 2;
      while exponent > 0 do
        base := base * base;
        if odd(exponent) then
          power := power * base;
        end if;
        exponent := exponent div 2;
      end while;
    end if;
  end func;

const proc: main is func
  local
    var matrix: m is [] (
      [] (4.0, 3.0),
      [] (2.0, 1.0));
    var integer: exponent is 0;
  begin
    for exponent range [] (0, 1, 2, 3, 5, 7, 11, 13, 17, 19, 23) do
      writeln("m ** " <& exponent <& " =");
      writeln(m ** exponent);
    end for;
  end func;

Original source of matrix exponentiation: [1]

Output:

m ** 0 =
1.0, 0.0
0.0, 1.0
m ** 1 =
4.0, 3.0
2.0, 1.0
m ** 2 =
22.0, 15.0
10.0, 7.0
m ** 3 =
118.0, 81.0
54.0, 37.0
m ** 5 =
3406.0, 2337.0
1558.0, 1069.0
m ** 7 =
98302.0, 67449.0
44966.0, 30853.0
m ** 11 =
81883680.0, 56183720.0
37455816.0, 25699956.0
m ** 13 =
2363278336.0, 1621541248.0
1081027456.0, 741736960.0
m ** 17 =
1968565387264.0, 1350712688640.0
900475125760.0, 617852567552.0
m ** 19 =
56815568027648.0, 38983467794432.0
25988979228672.0, 17832093941760.0
m ** 23 =
47326274699395072.0, 32472478198530048.0
21648320946503680.0, 14853792205897728.0

Sidef

class Array {
    method ** (Number n { .>= 0 }) {
        var tmp = self
        var out = self.len.of {|i| self.len.of {|j| i == j ? 1 : 0 }}
        loop {
            out = (out `mmul` tmp) if n.is_odd
            n >>= 1 || break
            tmp = (tmp `mmul` tmp)
        }
        return out
    }
}

var m = [[1, 2, 0],
         [0, 3, 1],
         [1, 0, 0]]

for order in (0..5) {
    say "### Order #{order}"
    var t = (m ** order)
    say ('  ', t.join("\n  "))
}

Output:

### Order 0
  [1, 0, 0]
  [0, 1, 0]
  [0, 0, 1]
### Order 1
  [1, 2, 0]
  [0, 3, 1]
  [1, 0, 0]
### Order 2
  [1, 8, 2]
  [1, 9, 3]
  [1, 2, 0]
### Order 3
  [3, 26, 8]
  [4, 29, 9]
  [1, 8, 2]
### Order 4
  [11, 84, 26]
  [13, 95, 29]
  [3, 26, 8]
### Order 5
  [37, 274, 84]
  [42, 311, 95]
  [11, 84, 26]

SPAD

Works with: FriCAS

Works with: OpenAxiom

Works with: Axiom

(1) -> A:=matrix [[0,-%i],[%i,0]]

        +0   - %i+
   (1)  |        |
        +%i   0  +
                                               Type: Matrix(Complex(Integer))
(2) -> A^4

        +1  0+
   (2)  |    |
        +0  1+
                                               Type: Matrix(Complex(Integer))
(3) -> A^(-1)

        +0   - %i+
   (3)  |        |
        +%i   0  +
                                     Type: Matrix(Fraction(Complex(Integer)))
(4) -> inverse A

        +0   - %i+
   (4)  |        |
        +%i   0  +
                          Type: Union(Matrix(Fraction(Complex(Integer))),...)

Domain:Matrix(R)

Stata

This implementation uses Exponentiation by squaring to compute a^n for a matrix a and an integer n (which may be positive, negative or zero).

real matrix matpow(real matrix a, real scalar n) {
	real matrix p, x
	real scalar i, s
	s = n<0
	n = abs(n)
	x = a
	p = I(rows(a))
	for (i=n; i>0; i=floor(i/2)) {
		if (mod(i,2)==1) p = p*x
		x = x*x
	}
	return(s?luinv(p):p)
}

Here is an example to compute Fibonacci numbers:

: matpow((0,1\1,1),10)
[symmetric]
        1    2
    +-----------+
  1 |  34       |
  2 |  55   89  |
    +-----------+

Tcl

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

package require Tcl 8.5
namespace path {::tcl::mathop ::tcl::mathfunc}

proc matrix_exp {m pow} {
    if { ! [string is int -strict $pow]} {
        error "non-integer exponents not implemented"
    }
    if {$pow < 0} {
        error "negative exponents not implemented"
    }
    lassign [size $m] rows cols
    # assume square matrix
    set temp [identity $rows]
    for {set n 1} {$n <= $pow} {incr n} {
        set temp [matrix_multiply $temp $m]
    }
    return $temp
}

proc identity {size} {
    set i [lrepeat $size [lrepeat $size 0]]
    for {set n 0} {$n < $size} {incr n} {lset i $n $n 1}
    return $i
}

% print_matrix [matrix_exp {{3 2} {2 1}} 1]
3 2 
2 1 
% print_matrix [matrix_exp {{3 2} {2 1}} 0]
1 0 
0 1 
% print_matrix [matrix_exp {{3 2} {2 1}} 2]
13 8 
 8 5 
% print_matrix [matrix_exp {{3 2} {2 1}} 3]
55 34 
34 21 
% print_matrix [matrix_exp {{3 2} {2 1}} 4]
233 144 
144  89 
% print_matrix [matrix_exp {{3 2} {2 1}} 10]
1346269 832040 
 832040 514229 

TI-89 BASIC

Built-in exponentiation:

[3,2;4,1]^4

Output: ${\begin{bmatrix}417&208\\416&209\end{bmatrix}}$

Ursala

For matrices of floating point numbers, the library function mmult can be used as shown. The user-defined id function takes a square matrix to the identity matrix of the same dimensions. The mex function takes a pair $(A,n)$ representing a real matrix $A$ and a natural exponent $n$ to the exponentiation $A^{n}$ using the naive algorithm.

#import nat
#import lin

id  = @h ^|CzyCK33/1.! 0.!*
mex = ||id@l mmult:-0^|DlS/~& iota

Alternatively, this version uses the fast binary algorithm.

mex = ~&ar^?\id@al (~&lr?/mmult@llPrX ~&r)^/~&alrhPX mmult@falrtPXPRiiX

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

#cast %eLLL

test = mex/*<<3.,2.>,<2.,1.>> <0,1,2,3,4,10>

output:

<
   <
      <1.000000e+00,0.000000e+00>,
      <0.000000e+00,1.000000e+00>>,
   <
      <3.000000e+00,2.000000e+00>,
      <2.000000e+00,1.000000e+00>>,
   <
      <1.300000e+01,8.000000e+00>,
      <8.000000e+00,5.000000e+00>>,
   <
      <5.500000e+01,3.400000e+01>,
      <3.400000e+01,2.100000e+01>>,
   <
      <2.330000e+02,1.440000e+02>,
      <1.440000e+02,8.900000e+01>>,
   <
      <1.346269e+06,8.320400e+05>,
      <8.320400e+05,5.142290e+05>>>

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.

Option Base 1
Private Function Identity(n As Integer) As Variant
    Dim I() As Variant
    ReDim I(n, n)
    For j = 1 To n
        For k = 1 To n
            I(j, k) = 0
        Next k
    Next j
    For j = 1 To n
        I(j, j) = 1
    Next j
    Identity = I
End Function
 Function MatrixExponentiation(ByVal x As Variant, ByVal n As Integer) As Variant
    If n < 0 Then
        x = WorksheetFunction.MInverse(x)
        n = -n
    End If
    If n = 0 Then
        MatrixExponentiation = Identity(UBound(x))
        Exit Function
    End If
    Dim y() As Variant
    y = Identity(UBound(x))
    Do While n > 1
        If n Mod 2 = 0 Then
            x = WorksheetFunction.MMult(x, x)
            n = n / 2
        Else
            y = WorksheetFunction.MMult(x, y)
            x = WorksheetFunction.MMult(x, x)
            n = (n - 1) / 2
        End If
    Loop
    MatrixExponentiation = WorksheetFunction.MMult(x, y)
End Function
Public Sub pp(x As Variant)
    For i_ = 1 To UBound(x)
        For j_ = 1 To UBound(x)
            Debug.Print x(i_, j_),
        Next j_
        Debug.Print
    Next i_
End Sub
Public Sub main()
    M2 = [{3,2;2,1}]
    M3 = [{1,2,0;0,3,1;1,0,0}]
    pp MatrixExponentiation(M2, -1)
    Debug.Print
    pp MatrixExponentiation(M2, 0)
    Debug.Print
    pp MatrixExponentiation(M2, 10)
    Debug.Print
    pp MatrixExponentiation(M3, 10)
End Sub

Output:

-1             2            
 2            -3            

 1             0            
 0             1            

 1346269       832040       
 832040        514229       

 13801         102408        31322        
 15661         116209        35543        
 4221          31322         9580 

Wren

Library: Wren-fmt

Library: Wren-matrix

Wren's Num class uses a method (pow) rather than an operator for exponentiation.

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 "./matrix" for Matrix
import "./fmt" for Fmt

var m = Matrix.new([[0, 1], [1, 1]])
System.print("Original:\n")
Fmt.mprint(m, 2, 0)
System.print("\nRaised to power of 10:\n")
Fmt.mprint(m ^ 10, 3, 0)

Output:

Original:

| 0  1|
| 1  1|

Raised to power of 10:

| 34  55|
| 55  89|