theory Gauss_Jordan_Imp
imports
"$ISABELLE_HOME/src/HOL/Main"
"$ISABELLE_HOME/src/HOL/Imperative_HOL/Imperative_HOL"
"$ISABELLE_HOME/src/HOL/Imperative_HOL/ex/Subarray"
"$ISABELLE_HOME/src/HOL/Multivariate_Analysis/Multivariate_Analysis"
"$ISABELLE_HOME/src/HOL/Library/Code_Target_Numeral"
"Code_Bit"
begin
instantiation bit::heap
begin
instance
proof
show " ∃to_nat::bit=>nat. inj to_nat"
by (rule exI[of _ "λn::bit. if n=0 then 0 else 1"], simp add: inj_on_def)
qed
end
fun swap :: "bit array => nat => nat => unit Heap" where
"swap a i j = do {
x \<leftarrow> Array.nth a i;
y \<leftarrow> Array.nth a j;
Array.upd i y a;
Array.upd j x a;
return ()
}"
definition "example = do {
a \<leftarrow> Array.of_list [42, 2, 3, 5, 0, 1705, 8, 3, 15];
swap a 0 1;
return a
}"
ML_val {* @{code example} ()*}
section{*Gauss-Jordan with mutable arrays*}
subsection{*Elementary operations*}
text{*Arrays must be proven to be instances of heap type class,
over any type being a heap itself.*}
instantiation array :: (heap) heap
begin
instance
proof
qed
end
text{*Interchange rows is equivalent to swap:*}
fun interchange_rows :: "(bit array) array => nat => nat => unit Heap" where
"interchange_rows M i j = do {
x \<leftarrow> Array.nth M i;
y \<leftarrow> Array.nth M j;
Array.upd i y M;
Array.upd j x M;
return ()
}"
definition "example2 = do {
a \<leftarrow> Array.of_list [42, 2, 3, 5/7, 0, 1705, 8];
b \<leftarrow> Array.of_list [24, 1, 4, 6, 9, 517, 98];
M \<leftarrow> Array.of_list [a, b, a, b];
interchange_rows M 0 3;
return M
}"
ML_val {* @{code example2} ()*}
text{*The following proof fails with "fun"; no argument decreases,
but the difference between "i" and "j":*}
text{*Another interesting point is that we use a second "aux" bit array;
this is to avoid later unwanted lateral effects:*}
text{*It may be more elegant to return a bit array, but then it should be built
inside of the recursive function, I guess...*}
function mult_array_rec :: "(bit array) => (bit array) => bit => nat => nat => unit Heap"
where
"mult_array_rec v w q i j = (if (i < j) then
do {
x \<leftarrow> Array.nth v i;
w \<leftarrow> Array.upd i (x * q) w;
mult_array_rec v w q (i + 1) j
}
else return ())"
by pat_completeness auto
termination by (relation "measure (λ(v,w,q,i,j). j - i)") simp+
text{*Using the previous function, imposing bounds 0 and "Array.len v",
and a newly allocated array, we recursively multiply the input array
by a constant:*}
fun mult_array :: "(bit array) => bit => (bit array) Heap" where
"mult_array v q =
do {
l \<leftarrow> Array.len v;
w \<leftarrow> Array.make l (%x. 0);
mult_array_rec v w q 0 l;
return w
}"
definition "example3 = do {
a \<leftarrow> Array.of_list [42, 2, 3, 5, 0, 1705, 8::bit];
w \<leftarrow> mult_array a 3;
return w
}"
ML_val {* @{code example3} ()*}
fun mult_row :: "(bit array) array => nat => bit => unit Heap" where
"mult_row M i q = do {
x \<leftarrow> Array.nth M i;
y \<leftarrow> mult_array x q;
Array.upd i y M;
return ()
}"
definition "example4 = do {
a \<leftarrow> Array.of_list [42, 2, 3, 5, 0, 1705, 8];
b \<leftarrow> Array.of_list [24, 1, 4, 6, 9, 517, 98];
M \<leftarrow> Array.of_list [a, b, a, b];
mult_row M 2 0.5;
return M
}"
ML_val {* @{code example4} ()*}
function add_array_rec :: "(bit array) => (bit array) => nat => nat => unit Heap" where
"add_array_rec v w i j = (if (i < j) then
do {
a \<leftarrow> Array.nth v i;
b \<leftarrow> Array.nth w i;
v' \<leftarrow> Array.upd i (a + b) v;
add_array_rec v' w (i + 1) j
}
else return ())"
by pat_completeness auto
termination by (relation "measure (λ(v,q,i,j). j - i)") simp+
text{*Using the previous function, imposing bounds 0 and "Array.len v",
we recursively add two arrays (they are supposed to be of equal length:)*}
fun add_array :: "(bit array) => (bit array) => unit Heap" where
"add_array v w =
do {
l \<leftarrow> Array.len v;
add_array_rec v w 0 l;
return()
}"
definition "example5 = do {
a \<leftarrow> Array.of_list [42, 2, 3, 5, 0, 1705, 8];
b \<leftarrow> Array.of_list [24, 1, 4, 6, 9, 517, 98];
add_array a b;
return a
}"
ML_val {* @{code example5} ()*}
text{*Row i is added row j times q:*}
fun row_add :: "(bit array) array => nat => nat => bit => unit Heap" where
"row_add M i j q = do {
a \<leftarrow> Array.nth M i;
b \<leftarrow> Array.nth M j;
c \<leftarrow> mult_array b q;
add_array a c;
Array.upd i a M;
return ()
}"
definition "example6 = do {
a \<leftarrow> Array.of_list [42, 2, 3, 5, 0, 1705, 8];
b \<leftarrow> Array.of_list [24, 1, 4, 6, 9, 517, 98];
c \<leftarrow> Array.of_list [24, 1, 4, 6, 9, 517, 98];
d \<leftarrow> Array.of_list [42::bit, 2, 3, 5, 0, 1705, 8];
M \<leftarrow> Array.of_list [a, b, c, d];
row_add M 1 2 0.5;
return M
}"
ML_val {* @{code example6} () *}
text{*In the following function the first condition might be removed by using a
tail recursive function, but I encountered some problems to prove it terminating.*}
function lnzero_position_f_ind :: "bit array => nat => nat => nat Heap" where
"lnzero_position_f_ind v i j = (if (i < j) then do {
x \<leftarrow> Array.nth v i;
(if (x ≠ 0) then (return i) else (lnzero_position_f_ind v (i + 1) j))
}
else (raise ''array lookup: index out of range''))"
by pat_completeness auto
termination by (relation "measure (λ(v,i,j). j - i)") simp+
fun lnzero_position_array_upt :: "(bit array) => nat => nat Heap" where
"lnzero_position_array_upt v i =
do {
l \<leftarrow> Array.len v;
j \<leftarrow> lnzero_position_f_ind v i l;
return j
}"
fun lnzero_position_array :: "(bit array) => nat Heap" where
"lnzero_position_array v =
do {
l \<leftarrow> Array.len v;
i \<leftarrow> lnzero_position_f_ind v 0 l;
return i
}"
definition "example7 = do {
v \<leftarrow> Array.of_list [0, 0, 0, 5, 0, 0, 9];
b \<leftarrow> lnzero_position_array_upt v 4;
return b
}"
ML_val {* @{code example7} ()*}
text{*A predicate to detect columns array where there is no element to
be pivoted.*}
function pivot_rec :: "(bit array) => nat => nat => bool Heap"
where
"pivot_rec v i j = (if (i < j) then do {
x \<leftarrow> Array.nth v i;
(if (x = 0) then pivot_rec v (i + 1) j else return True)
}
else return False)"
by pat_completeness auto
termination by (relation "measure (λ(v,i,j). j - i)") simp+
fun pivot :: "(bit array) => nat => bool Heap"
where
"pivot v i = do {
l \<leftarrow> Array.len v;
b \<leftarrow> pivot_rec v i l;
return b
}"
definition "example8 = do {
v \<leftarrow> Array.of_list [0, 0, 0, 5, 0, 0, 7];
b \<leftarrow> pivot v 4;
return b
}"
ML_val {* @{code example8} ()*}
text{*The next definition requires that the element in position M i j is equal to one;
therefore, we must first pivot an element into position M i j, and then invoke it.*}
text{*Using row i, every row in the interval k to l is reduced in column j.*}
function row_reduce_rec :: "(bit array) array => nat => nat => nat => nat => unit Heap"
where
"row_reduce_rec M i j k l = (if (k < l) then do
{
(if (k = i) then row_reduce_rec M i j (k + 1) l else do {
fk \<leftarrow> Array.nth M k;
mkj \<leftarrow> Array.nth fk j;
row_add M k i (- mkj);
row_reduce_rec M i j (k + 1) l
}
)
}
else return ())"
by (pat_completeness) auto
termination by (relation "measure (λ(M,i,j,k,l). l - k)") simp+
text{*Applying row reduction with row i in column j*}
fun row_reduce :: "(bit array) array => nat => nat => unit Heap"
where
"row_reduce M i j = do
{
nr \<leftarrow> Array.len M;
row_reduce_rec M i j 0 nr
}"
definition "example9 = do {
a \<leftarrow> Array.of_list [42::bit, 2, 3, 5];
b \<leftarrow> Array.of_list [24::bit, 1, 4, 6];
c \<leftarrow> Array.of_list [24::bit, 1, 4, 6];
d \<leftarrow> Array.of_list [42::bit, 2, 3, 5];
M \<leftarrow> Array.of_list [a, b, c, d];
row_reduce M 1 1;
return M
}"
ML_val {* @{code example9} ()*}
text{*We need an operation to get a column of a matrix; we use an auxiliary array
in the second parameter.*}
function column_rec :: "(bit array) array => bit array => nat => nat => nat => unit Heap"
where
"column_rec M w k i j = (if (i < j) then do {
fi \<leftarrow> Array.nth M i;
mij \<leftarrow> Array.nth fi k;
Array.upd i mij w;
column_rec M w k (i + 1) j
}
else return ())"
by (pat_completeness) auto
termination by (relation "measure (λ(M,w,k,i,j). j - i)") simp+
fun column :: "(bit array) array => nat => (bit array) Heap"
where
"column M k = do {
l \<leftarrow> Array.len M;
w \<leftarrow> Array.make l (%x. 0);
column_rec M w k 0 l;
return w
}"
definition "example10 = do {
a \<leftarrow> Array.of_list [42::bit, 2, 3, 5];
b \<leftarrow> Array.of_list [24::bit, 1, 4, 6];
c \<leftarrow> Array.of_list [24::bit, 1, 4, 6];
d \<leftarrow> Array.of_list [42::bit, 2, 3, 5];
M \<leftarrow> Array.of_list [a, b, c, d];
w \<leftarrow> column M 3;
return w
}"
ML_val {* @{code example10} ()*}
text{*Still need to make zeros above and below the pivot;*}
fun Gauss_Jordan_in_ij :: "(bit array) array => nat => nat => unit Heap"
where
"Gauss_Jordan_in_ij M i j = do {
cj \<leftarrow> column M j;
n \<leftarrow> lnzero_position_array_upt cj i;
pivot \<leftarrow> Array.nth cj n;
interchange_rows M i n;
mult_row M i (1 / pivot);
row_reduce M i j;
return ()
}"
definition "example11 = do {
a \<leftarrow> Array.of_list [42, 2, 3, 5];
b \<leftarrow> Array.of_list [24, 1, 4, 6];
c \<leftarrow> Array.of_list [24, 1, 0, 6];
d \<leftarrow> Array.of_list [42, 2, 3, 5];
M \<leftarrow> Array.of_list [a, b, c, d];
Gauss_Jordan_in_ij M 2 2;
return M
}"
ML_val {* @{code example11} ()*}
fun Gauss_Jordan_column_k :: "((bit array) array × nat) => nat => ((bit array) array × nat) Heap"
where
"Gauss_Jordan_column_k Mp k = do {
ck \<leftarrow> column (fst Mp) k;
b \<leftarrow> pivot ck (snd Mp);
nr \<leftarrow> Array.len (fst Mp);
(if ((¬ b) ∨ (nr = (snd Mp))) then return (((fst Mp), (snd Mp)))
else do {
Gauss_Jordan_in_ij (fst Mp) (snd Mp) k;
return (((fst Mp), (snd Mp) + 1))
})
}"
definition "example12 = do {
a \<leftarrow> Array.of_list [42::bit, 2, 3, 10];
b \<leftarrow> Array.of_list [24::bit, 1, 4, 6];
c \<leftarrow> Array.of_list [24::bit, 1, 0, 6];
d \<leftarrow> Array.of_list [42::bit, 2, 3, 5];
M \<leftarrow> Array.of_list [a, b,c, d];
Gauss_Jordan_column_k (M, 0) 0;
return M
}"
ML_val {* @{code example12} ()*}
text{*In the following function k denotes the last column to be processed,
i denotes the column being processed in each recursive call,
and j the row in which the pivot in row i is being looked for.*}
function Gauss_Jordan_upt_k_rec :: "(bit array) array => nat => nat => nat => unit Heap"
where
"Gauss_Jordan_upt_k_rec M k i j = (if (i < k) then do {
Mp \<leftarrow> Gauss_Jordan_column_k (M, j) i;
Gauss_Jordan_upt_k_rec (fst Mp) k (i + 1) (snd Mp)
} else return ())"
by (pat_completeness) auto
termination by (relation "measure (λ(M,k,i,j). k - i)") simp+
text{*Using the previous function, starting from column 0 and row 0, we define a new function
to produce Gauss-Jordan up to a column k.*}
fun Gauss_Jordan_upt_k :: "(bit array) array => nat => unit Heap"
where
"Gauss_Jordan_upt_k M k = Gauss_Jordan_upt_k_rec M k 0 0"
definition "example13 = do {
a \<leftarrow> Array.of_list [42, 2, 3, 10];
b \<leftarrow> Array.of_list [24, 1, 4, 6];
c \<leftarrow> Array.of_list [24, 1, 0, 6];
d \<leftarrow> Array.of_list [42, 2, 3, 5];
M \<leftarrow> Array.of_list [a, b, c, d];
Gauss_Jordan_upt_k M 3;
return M
}"
ML_val {* @{code example13} ()*}
text{*Note that we only count the elements in the first row to compute the number of columns;
input matrices have to be well formed.*}
fun cols :: "(bit array) array => nat Heap"
where "cols M = do {
v \<leftarrow> Array.nth M 0;
n \<leftarrow> Array.len v;
return n
}"
text{*I am not pretty sure if the column to stop is ncols or "ncols - 1"*}
fun Gauss_Jordan :: "(bit array) array => unit Heap"
where
"Gauss_Jordan M = do {
k \<leftarrow> cols M;
Gauss_Jordan_upt_k M k;
return ()
}"
definition "example14 = do {
a \<leftarrow> Array.of_list [42, 2, 3, 10];
b \<leftarrow> Array.of_list [24, 1, 4, 6];
c \<leftarrow> Array.of_list [24, 1, 0, 6];
d \<leftarrow> Array.of_list [42, 2, 3, 5];
M \<leftarrow> Array.of_list [a, b, c, d];
Gauss_Jordan M;
return M
}"
ML_val {* @{code example14} ()*}
definition "example15 = do {
a \<leftarrow> Array.of_list [1, 0, 0, 0];
b \<leftarrow> Array.of_list [0, 1, 0, 0];
c \<leftarrow> Array.of_list [0, 0, 1, 0];
d \<leftarrow> Array.of_list [1, 0, 1, 0];
M \<leftarrow> Array.of_list [a, b, c, d, a, c, d];
Gauss_Jordan M;
return M
}"
ML_val {* @{code example15} ()*}
export_code Gauss_Jordan example14 in SML
file "Gauss_Jordan_imp.sml"
definition "matrix_z2 = do {
a \<leftarrow> Array.of_list [0, 1, 0, 1::bit];
b \<leftarrow> Array.of_list [0, 1, 1, 1];
c \<leftarrow> Array.of_list [1, 1, 1, 1];
d \<leftarrow> Array.of_list [1, 0, 0, 1];
M \<leftarrow> Array.of_list [a, b, c, d];
return M
}"
definition "print_result_Gauss M = do {
Gauss_Jordan M;
return M
}"
term Gauss_Jordan
term print_result_Gauss
export_code matrix_z2 print_result_Gauss Gauss_Jordan in SML
module_name "Gauss_SML"
file "Gauss_Jordan_Imp.sml"
end