Matrix-exponentiation operator: Difference between revisions
Content added Content deleted
(add RPL) |
|||
Line 332: | Line 332: | ||
( 0.0000+0.0000i, 1.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)); |
( 0.0000+0.0000i, 0.0000+0.0000i, 1.0000+0.0000i)); |
||
</pre> |
|||
=={{header|ATS}}== |
|||
<syntaxhighlight lang="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 |
|||
(*------------------------------------------------------------------*) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>$ 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 |
|||
</pre> |
</pre> |
||