55 lines
1.2 KiB
Text
55 lines
1.2 KiB
Text
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))
|