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.




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