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, sequence b)
sequence c
    if length(a[1]) != length(b) then
        return 0
    else
        c = repeat(repeat(0,length(b[1])),length(a))
        for i=1 to length(a) do
            for j=1 to length(b[1]) do
                for k=1 to length(a[1]) do
                    c[i][j] += a[i][k]*b[k][j]
                end for
            end for
        end for
        return c
    end if
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}}
constant M2 = {{3, 2},
               {2, 1}}
constant 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))