Matrix-exponentiation operator: Difference between revisions

=={{header|Phix}}==

Phix does not permit operator overloading, however here is a simple function to raise a square matrix to a non-negative integer power.<br>

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

<lang Phix>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))</lang>

{{out}}

<pre>

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

</pre>