type linear_congruence = { A : int [] []; b : int []; m : int }

module LinearCongruence = struct
  type t = linear_congruence 

  let wf_coefficient m coefficient  = 0 <= coefficient && coefficient < m 

  let wf lc =
    let in_bounds coefficient = wf_coefficient lc.m coefficient in
    Array.length lc.A = Array.length lc.b &&
    Array.for_all (Array.for_all in_bounds) lc.A &&
    Array.for_all in_bounds lc.b

  let printf lc = 
    let mid = Array.length lc.A/2 in
    Array.iteri 
      (fun idx vec -> 
        printf "[";
        Array.iter (printf " %d") vec;
        printf " ]";
        if idx = mid then 
         printf " = [ %d ] mod %d\n" lc.b.[idx] lc.m 
	else 
	  printf "   [ %d ]\n" lc.b.[idx]
      )
      lc.A
      
end


// ---------- Calculate Join of Two Congruences -------------------


// Vectors
let zeroVec num = Array.init num (fun i -> 0)
let oneVec num = Array.init num (fun i -> 1)
let negateVector m vec = Array.map (fun v -> (m - v) % m) vec
let getDimension0 m = Array.length m
let printVector vec = printf "[ "; Array.iter (printf "%d ") vec; printf "]\n"


// Matrices
let printMatrix m = Array.iter printVector m
let ithPlace dim place value = Array.init dim (fun x -> if x <> place then 0 else value)
let diagMatrix dim v = Array.init dim (fun row -> ithPlace dim row v)
let idMatrix dim = diagMatrix dim 1
let zeroMatrix rows cols = Array.init rows (fun r -> Array.create cols 0)
let negateMatrix m matrix = Array.map (negateVector m) matrix
let transpose vec = Array.map (fun v -> [|v|] ) vec
let zeroCol dim = transpose (zeroVec dim)
let oneCol  dim = transpose (oneVec dim)

let getDimension1 m = if getDimension0 m = 0 then 0 else getDimension0 m.[0]
let getDimensions m = (getDimension0 m, getDimension1 m)

let rec power2 x = if x = 0 then 1 else 2*(power2 (x-1))

let rec log2_ceiling r =  
  match r with
  | 0 -> 0
  | 1 -> 0
  | _ -> (log2_ceiling (r / 2) ) + 1


let norm num m = 
  let n = num % m in 
  if n < 0 then n + m else n
	

let (@@) a b = Array.append a b
let (@@@) m1 m2 = 
  assert (getDimension0 m1 = getDimension0 m2);
  Array.map2 Array.append m1 m2

(*
 * Given A1 x = b1, A2 x = b2 create the matrix: 
 *
 * [ 1     1   0   0  0 1]
 * [ -b1   0  A1   0  0 0]
 * [ 0   -b2   0  A2  0 0]
 * [ 0     0  -I  -I  I 0]
 *)
let createJointMatrix lc1 lc2 =
  assert (getDimension1 lc1.A >= getDimension1 lc2.A);
  let dim = getDimension1 lc1.A in
  let m = lc1.m in
  let (rows1, cols1) = getDimensions lc1.A in
  let (rows2, cols2) = getDimensions lc2.A in
  let b1 = transpose (negateVector m lc1.b) in
  let b2 = transpose (negateVector m lc2.b) in
  let zV dim = zeroVec dim in
  let zC dim = zeroCol dim in 
  let dM dim v = diagMatrix dim v in
  let zM rows cols = zeroMatrix rows cols in
  let matrix = 
    [| [| 1;       1 |]     @@  zV cols1      @@  zV cols2       @@  zV dim       @@ [| 1|] |]  @@
    (b1       @@@ zC rows1 @@@ lc1.A          @@@ zM rows1 cols2 @@@ zM rows1 dim @@@ zC rows1) @@
    (zC rows2 @@@ b2       @@@ zM rows2 cols1 @@@ lc2.A          @@@ zM rows2 dim @@@ zC rows2) @@
    (zC dim   @@@ zC dim   @@@ dM dim (m-1)   @@@ dM dim (m-1)   @@@ dM dim 1     @@@ zC dim)
  in
  //printf "joint matrix\n";
  //printMatrix matrix;
  matrix

let leading vec = Array.find_index (fun x -> x <> 0) vec 
    
let vecDiff vec vec1 = Array.map2 (fun x y -> x - y) vec vec1
				    
let multScalVec num vec = Array.map (fun x -> num * x) vec

let normVec m vec = Array.map (fun x -> norm x m) vec

let incr1 (a, b) = (a+1, b)

let rec rank num = 
  match num with 
  | 0 -> (0, 0)
  | _ when num % 2 <> 0 -> (0, num)
  | _ -> incr1 (rank (num / 2))

let leadingIndexEntryRankNonInvert vec =
  let index = leading vec     in
  let entry = vec.[index]     in
  let (rnk, inv) = rank entry in
  (index, entry, rnk, inv)

let is_zero vec = Array.for_all (fun x -> x = 0) vec

let is_non_zero vec =  not (is_zero vec)

let addInSet set inx vec = 
  if 0 <= inx && inx < Array.length set then 
    set.[inx] <- vec


let coeff r r1 i = (power2 (r - r1)) * i

let updateMatrix mat = Array.filter is_non_zero mat 

let rec triangulateSet (set : int [] []) vec m = 
  if is_non_zero vec then 
    let (inx, ent, rnk, inv) = leadingIndexEntryRankNonInvert vec in
    assert (0 <= inx && inx < Array.length set);
    if is_zero set.[inx] then 
      addInSet set inx vec
    else      
      let vec1 = set.[inx] in
      let (inx1, ent1, rnk1, inv1) = leadingIndexEntryRankNonInvert vec1 in
      if rnk1 <= rnk then 
	let vec2 = normVec m (vecDiff (multScalVec inv1 vec) (multScalVec (coeff rnk rnk1 inv) vec1)) in
	triangulateSet set vec2 m
      else 
	let vec2 = normVec m (vecDiff (multScalVec inv vec1) (multScalVec (coeff rnk1 rnk inv1) vec)) in
	(addInSet set inx vec);
        triangulateSet set vec2 m
            
 
let triangulate ( mat : int [] []) m = 
  let num_rows = getDimension0 mat in
  let num_cols = getDimension1 mat in
  let set =  Array.init num_cols (fun _ -> [||]) in
  for i = 0 to num_rows - 1 do
    if is_non_zero mat.[i] then       
      let l = leading mat.[i] in
      if set.[l] = [||] then 
	addInSet set l mat.[i] 
      else
	triangulateSet set mat.[i] m
  done;
  updateMatrix set

let findDimForProjection (mat: int [] []) dim = 
  let n = Array.length mat.[0] in
  let mutable i = Array.length mat - 1 in
  let mutable j = 0 in
  while (0 <= i) do
    if (is_zero (Array.sub mat.[i] 0 (n - 1 - dim))) then i <- i - 1
    else begin
      j <- i; i <- -2
    end
  done;
  //
  // (j = 0, i = - 2 ) = return 1 if there are at least two rows
  //      if is_zero (Array.sub mat.[i] 0 (n - 1 - dim) holds, it is covered
  //      if is_zero (Array.sub mat.[i] 0 (n - 1 - dim) does not holds, it is degeneric case, return -1
  // (j = 0, i = - 1 ) = return 0
  // (j = k > 0, i = - 2 ) = return k + 1 if not last, if last return number of rows
  //
  if (i = -1) then 
    0
  else if (j = (Array.length mat - 1)) then 
    if (is_zero (Array.sub mat.[j] 0 (n - 1 - dim))) then 
      j 
    else -1
  else j + 1


let projection (mat: int [] []) m n =
  let lM = getDimension0 mat - 1 in
  let lN = getDimension1 mat - 1 in
  let mat1 = Array.init ( lM - m + 1) (fun _ -> [||]) in
  let b1 = Array.init ( lM - m  + 1) (fun _ -> 0 ) in
  for i = m to lM  do
    let vec = Array.sub mat.[i] (lN - n ) n in
    mat1.[i - m] <- vec
  done;
  for i = m to lM do
    b1.[i - m] <- mat.[i].[lN]
  done;
  (mat1, b1)

let constructTrue dim m1 = { A=[| zeroVec dim |]; b=[| 0 |]; m=m1 }

let constructFalse dim m1 = { A=[| zeroVec dim |]; b=[| 1 |]; m=m1 }

let isZeroMatrix mat = Array.for_all (fun x -> is_zero x) mat

let isTrue lc = isZeroMatrix lc.A && is_zero lc.b

let isFalse lc = isZeroMatrix lc.A && is_non_zero lc.b

let constructNonTrivialJoin lc1 lc2 =
  assert (lc1.m = lc2.m);
  let dim = getDimension1 lc1.A in
  let m1 = triangulate (createJointMatrix lc1 lc2) lc1.m in 
  //printf "triangulated matrix\n";
  //printMatrix m1;
  let rawNum = findDimForProjection m1 dim in
  if (rawNum = -1) then 
    constructTrue dim lc1.m 
  else
    let (matJ, vecJ) = projection m1 rawNum dim in
    { 
      A = matJ; 
      b = vecJ; 
      m = lc1.m
    }

let constructJoin lc1 lc2 =
  if isTrue lc1 || isTrue lc2 then 
    constructTrue (getDimension1 lc1.A) lc1.m
  else if isFalse lc1 && isFalse lc2 then 
    constructFalse (getDimension1 lc1.A) lc1.m
  else if isFalse lc1 then 
    lc2
  else if isFalse lc2 then 
    lc1 
  else 
    constructNonTrivialJoin lc1 lc2


// END OF:   Calculate Join of Two Congruences -------------------


// ------  From Control Flow Graph to Formulas


open Microsoft.Z3
open System.Collections.Generic
open System

let par = new Config()
do  par.SetParamValue("MODEL", "true") 
let z3 = new TypeSafeContext(par)


type location = int
type var = int
type term = | Var of var 
            | Add of term * term
            | Sub of term * term
            | BvAnd of term * term
            | True
            | False
            | Num of int
            | Eq of term * term
            | Not of term
type guard      = term
type assignment = var * term
type command    = guard * assignment list
type edge       = { l_pre : location; l_post : location; command : command }

// location 0 is the initial location.
type program = { num_locations : int; num_vars : int; edges : edge list; m : int }


// mk_term takes a term and returns z3 term
let rec mk_term (z3:TypeSafeContext) ty term = 
  match term with 
  | Var v -> z3.MkConst(z3.MkSymbol(v), ty)
  | Add (t1,t2) -> z3.MkBvAdd ((mk_term z3 ty t1), (mk_term z3 ty t2))
  | Sub (t1, t2) -> z3.MkBvSub ((mk_term z3 ty t1), (mk_term z3 ty t2))
  | BvAnd (t1, t2) -> z3.MkBvAnd((mk_term z3 ty t1), (mk_term z3 ty t2))
  | True -> z3.MkTrue()
  | False -> z3.MkFalse() 
  | Num n -> z3.MkNumeral(n, ty)
  | Eq (t1, t2) -> z3.MkEq((mk_term z3 ty t1), (mk_term z3 ty t2))
  | Not t1 -> z3.MkNot (mk_term z3 ty t1)


let mk_formula (z3:TypeSafeContext) (program:program) ty (edge:edge) = 
  let guard, assignment = edge.command in
  let mk_var v = z3.MkConst(z3.MkSymbol(v + program.num_vars),ty) in
  let assignment = List.map (fun (v,t) -> z3.MkEq(mk_var v, mk_term z3 ty t)) assignment in
  let assignment = Array.of_list assignment in
  z3.MkAnd(mk_term z3 ty guard, z3.MkAnd assignment)


type todo_list = class
    val mutable list : int list;
    val active : bool array;
    new (num_items) = { list = []; active = Array.create num_items false }

    member this.add item = 
      if not this.active.[item] then 
        (this.active.[item] <- true; this.list <- item::this.list)

    member this.pop = 
      match this.list with 
      | [] -> None 
      | item::rest -> (this.list <- rest; this.active.[item] <- false; Some item)

    member this.print =
      List.iter (printf "%d ") this.list; printf "\n"
  end

     
// mapping from program locations to linear congruences.
// if we have multiple contexts, then we also need: contexts : TypeSafeContext array

type abstract_interpretation_state = 
  { program     : program;
    congruences : LinearCongruence.t array;  // one linear congruence per program location.
    todo        : todo_list;
  }


// mk_initial_congruence (program:program) -> int(location) -> LinearCongruence.t
let mk_initial_congruence (program:program) location = 
  if location = 0 then 
    constructTrue (program.num_vars * (log2_ceiling program.m)) program.m
  else 
    constructFalse program.num_vars program.m


// mk_initial_abstract_interpretation_state -> program -> abstract_interpretation_stat
let mk_initial_abstract_interpretation_state program = 
  { 
    program = program; 
    congruences = Array.init program.num_locations (fun idx -> mk_initial_congruence program idx);
    todo = let todo = new todo_list(program.num_locations) in todo.add(0); todo
  }

let get_primed_vars (program:program) = 
  let m = program.m in
  let ty = z3.MkBvType(uint32 (log2_ceiling m)) in
  List.init program.num_vars (fun idx -> z3.MkConst(z3.MkSymbol(program.num_vars + idx), ty)) 

let mk_formulaFromCongruence offset (lc: LinearCongruence.t) = 
  assert (LinearCongruence.wf lc);
  let ty = z3.MkBvType(uint32 (log2_ceiling lc.m)) in
  let mk_var (var:int) = z3.MkConst(z3.MkSymbol(offset + var), ty) in
  let add a b = match a with None -> b | Some a -> z3.MkBvAdd(a,b) in
  let getOption a = match a with None -> z3.MkNumeral(0, ty) | Some a -> a in
  let mk_eq idx row = 
    let mutable a = None in
    for i = 0 to Array.length row - 1 do
      match row.[i] with
      | 0 -> ()
      | 1 -> a <- Some (add a (mk_var i))
      | c -> a <- Some (add a (z3.MkBvMul(z3.MkNumeral(c, ty), mk_var i)))
    done;
    z3.MkEq(getOption a, z3.MkNumeral(lc.b.[idx],ty))
  in
  z3.MkAnd (Array.mapi mk_eq lc.A)
  

(*
 * Convert a model to a linear congruence.
 *
 * create diagonal array
 * create column vector of values for the arrays.
 * remember modulus.
 *)
let extract_linear_congruence (model:TypeSafeModel) vars m =  
  let dimension = List.length vars               in
  let extract_val var = model.GetNumeralValueInt(model.Eval(var)) in 
  let vals = List.map extract_val vars           in
  assert (List.for_all (LinearCongruence.wf_coefficient m) vals);
  { 
    A = idMatrix dimension;
    b = Array.of_list vals;
    m = m
   }

let mk_formulaFromCongruence_bitwise offset (lc: LinearCongruence.t) = 
  assert (LinearCongruence.wf lc);
  let bv_size = log2_ceiling lc.m                      in
  assert ((getDimension1 lc.A) % bv_size = 0);
  let ty = z3.MkBvType(uint32 bv_size)                 in
  let mk_var (id:int) =
    let var = id / bv_size                             in
    let idx = id % bv_size                             in
    let v  = z3.MkConst(z3.MkSymbol(offset + var), ty) in
    z3.MkBvZeroExt(uint32 (bv_size-1), z3.MkBvExtract(uint32 idx, uint32 idx, v))
  in
  let add a b = match a with None -> b | Some a -> z3.MkBvAdd(a,b) in
  let getOption a = match a with None -> z3.MkNumeral(0, ty) | Some a -> a in
  let mk_eq idx row = 
    let mutable a = None in
    for i = 0 to Array.length row - 1 do
      match row.[i] with
      | 0 -> ()
      | 1 -> a <- Some (add a (mk_var i))
      | c -> a <- Some (add a (z3.MkBvMul(z3.MkNumeral(c, ty), mk_var i)))
    done;
    z3.MkEq(getOption a, z3.MkNumeral(lc.b.[idx],ty))
  in
  z3.MkAnd (Array.mapi mk_eq lc.A)

let extract_linear_congruence_bitwise (model:TypeSafeModel) vars m =  
  let bv_size = log2_ceiling m                                    in
  let dimension = List.length vars                                in
  let extract_val var = model.GetNumeralValueInt(model.Eval(var)) in 
  let expand_val v = List.init bv_size (fun idx -> if 0 <> ((1 <<< idx) &&& v) then 1 else 0) in
  let vals = List.map (extract_val >> expand_val) vars            in
  let vals = List.concat vals                                     in
  assert (List.for_all (LinearCongruence.wf_coefficient m) vals);
  { 
    A = idMatrix (dimension*bv_size);
    b = Array.of_list vals;
    m = m
   }
val map : ('a -> 'b) -> 'a list -> 'b list

let new_congruence lc = 
  let linear_congruence = mk_formulaFromCongruence 0 lc in // TBD: for incrementality, probably we want the primed ones
  z3.AssertCnstr(z3.MkNot(linear_congruence))
              
let printAbstractInterpretationState s = 
  let congs = s.congruences in
  for j = 0 to congs.Length - 1 do
    printf "Location %d:\n" (j);
    LinearCongruence.printf (congs.[j])
  done;
  printf "todo list: "; s.todo.print
    
    
let propagate_edge location (state:abstract_interpretation_state) (edge:edge) = 
  if location = edge.l_pre then 
    let formulaPre =  mk_formulaFromCongruence_bitwise 0 state.congruences.[location] in
    let ty = z3.MkBvType(uint32 (log2_ceiling state.program.m)) in
    let formulaEdge =  mk_formula z3 state.program ty edge in
    let formulaPost =  mk_formulaFromCongruence_bitwise state.program.num_vars state.congruences.[edge.l_post] in
    let formula = z3.MkAnd [| formulaPre ; formulaEdge; z3.MkNot(formulaPost) |] in
    let model = ref (null : TypeSafeModel) in
    Printf.printf "%O\n" formula;
    z3.Push();
    z3.AssertCnstr(formula);
    (match z3.CheckAndGetModel(model) with
    | LBool.True -> 
	(*
	 * extract the linear congruence implied by the current model.
	 * join the extracted linear congruence with the current congruence of a state
	 *)
        let model = !model in
	printf "Model Found\n"; // current context is satisfiable.
	assert (model <> null);
	model.Display(Console.Out);
        let p_vars = get_primed_vars state.program in
        let lc1 = extract_linear_congruence_bitwise model p_vars state.program.m  in
        let lc2 = constructJoin lc1 state.congruences.[edge.l_post] in

        printf "congruence from model:\n"; 
	LinearCongruence.printf lc1;
        printf "congruence from state:\n"; 
	LinearCongruence.printf state.congruences.[edge.l_post];
        printf "their join: \n";  
	LinearCongruence.printf lc2;

        state.congruences.[edge.l_post] <- lc2;
        state.todo.add edge.l_post;
    	state.todo.add edge.l_pre;
        model.Dispose();

	(*
         * For incremental version:
	 * block the additional constraints from future models.
         * new_congruence lc1
	 *)
    | LBool.False -> 
	printf "fixed point\n" //Nothing, fixed point reached
    | _ -> 
	assert false 
    );
    z3.Pop()

    
let propagate_location location (state:abstract_interpretation_state) = 
  List.iter (propagate_edge location state) state.program.edges

let rec propagate (state:abstract_interpretation_state) i = 
  match state.todo.pop with
  | None -> ()
  | Some location -> 
      propagate_location location state;
      printf "end of iteration %d \n" i;
      printAbstractInterpretationState state;
      propagate state (i + 1)

let main_call program = 
  propagate (mk_initial_abstract_interpretation_state  program) 0

// -------------------------------------------------------------------
// Test zone


let sample_program = 
  let varX = 0 in
  let varY = 1 in
  let varC = 2 in
  { num_locations = 3;
    num_vars = 3;
    edges = 
    [
     {l_pre = 0; l_post = 1; command = (True, [(varY, Var varX); (varC, Num 0); (varX, Var varX)]) };
     {l_pre = 1; l_post = 1; command = (Not(Eq(Var varY, Num 0)), [(varY, BvAnd(Var varY, Sub(Var varY, Num 1))); (varC, Add(Var varC, Num 1)); (varX, Var varX)]) };
     {l_pre = 1; l_post = 2; command = (Eq(Var varY, Num 0), [(varX, Var varX); (varY, Var varY); (varC, Var varC)]) }
   ];
    m = 256;
  }
  
let test1() = main_call sample_program

let test_join lc1 lc2 = 
  let lc3 = constructJoin lc1 lc2 in
  printf "lc1:\n";
  LinearCongruence.printf lc1;
  printf "lc2:\n";
  LinearCongruence.printf lc2;
  printf "join:\n";
  LinearCongruence.printf lc3

let test2() = 
  let a1 = [| [| 1;0;0 |]; [| 0;1;0 |]; [| 0;0;1 |] |] in
  let a2 = [| [| 1;0;0 |]; [| 0;1;0 |]; [| 0;0;1 |] |] in
  let b1 = [| 1;1;0 |]                  in
  let b2 = [| 0;0;0|]                   in
  let lc1 = { A = a1; b = b1; m = 256 } in
  let lc2 = { A = a2; b = b2; m = 256 } in
  test_join lc1 lc2

let test3() = 
  let matA = [| [| 1; 0; 0 |]; [| 0; 1; 0 |]; [| 0; 0; 1 |] |] in
  let matB = [| [| 0; 0; 0 |]  |] in
  let vecA = [| 0; 255; 0 |] in
  let vecB = [| 0 |] in
  let lc1 = { A=matA; b=vecA; m=256} in
  let lc2 = { A=matB; b=vecB; m=256} in
  test_join lc1 lc2

let test4() = 
  let fml =  mk_formulaFromCongruence 0 (constructTrue 3 4) in
  Printf.printf "%O\n" fml

let test5() = 
  let a1 = [| [| 1;0;0 |]; [| 0;1;0 |]; |] in
  let a2 = [| [| 1;0;0 |]; [| 0;1;0 |]; [| 0;0;1 |] |] in
  let b1 = [| 0;1 |]                  in
  let b2 = [| 1;0;1|]                 in
  let lc1 = { A = a1; b = b1; m = 2 } in
  let lc2 = { A = a2; b = b2; m = 2 } in
  test_join lc1 lc2

let test6() = 
  let a1 = [| [| 255;1;0 |]; [| 0;0;1 |] |] in
  let a2 = [| [| 1;0;0 |]; [| 0;1;0 |]; [| 0;0;1 |] |] in
  let b1 = [| 0; 0 |]                  in
  let b2 = [| 128;0;1|]                 in
  let lc1 = { A = a1; b = b1; m = 256 } in
  let lc2 = { A = a2; b = b2; m = 256 } in
  test_join lc1 lc2

(*

   a1: 

   xi = yi 
   ci = 0

   1; 0 ..., 255;  0; ... = 0
   0; 1; ...; 0; 255; ... = 0
   
   0               1 0 0  = 0
   0               0 1 0  = 0
   0; 0; ... ; 0   0 0 1  = 0

   x = 128
   y = 0
   c = 1
*)



let test7() = 
  let a1 = 
   [| 
   [| 1;0;0;0; 0;0;0;0; 255;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0 |]; 
   [| 0;1;0;0; 0;0;0;0; 0;255;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0 |]; 
   [| 0;0;1;0; 0;0;0;0; 0;0;255;0; 0;0;0;0; 0;0;0;0; 0;0;0;0 |]; 
   [| 0;0;0;1; 0;0;0;0; 0;0;0;255; 0;0;0;0; 0;0;0;0; 0;0;0;0 |]; 
   [| 0;0;0;0; 1;0;0;0; 0;0;0;0; 255;0;0;0; 0;0;0;0; 0;0;0;0 |]; 
   [| 0;0;0;0; 0;1;0;0; 0;0;0;0; 0;255;0;0; 0;0;0;0; 0;0;0;0 |]; 
   [| 0;0;0;0; 0;0;1;0; 0;0;0;0; 0;0;255;0; 0;0;0;0; 0;0;0;0 |]; 
   [| 0;0;0;0; 0;0;0;1; 0;0;0;0; 0;0;0;255; 0;0;0;0; 0;0;0;0 |]; 
   [| 0;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0;   1;0;0;0; 0;0;0;0 |];
   [| 0;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0;   0;1;0;0; 0;0;0;0 |];
   [| 0;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0;   0;0;1;0; 0;0;0;0 |];
   [| 0;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0;   0;0;0;1; 0;0;0;0 |];
   [| 0;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0;   0;0;0;0; 1;0;0;0 |];
   [| 0;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0;   0;0;0;0; 0;1;0;0 |];
   [| 0;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0;   0;0;0;0; 0;0;1;0 |];
   [| 0;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0;   0;0;0;0; 0;0;0;1 |]
   |]
  in
  let a2 = diagMatrix 24 1 in   
  let b1 = [| 0;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0 |]                   in
  let b2 = [| 1;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;0; 0;0;0;1 |]            in
  let lc1 = { A = a1; b = b1; m = 256 } in
  let lc2 = { A = a2; b = b2; m = 256 } in
  test_join lc1 lc2

(*
do test2()

do test3()

do test4()

do test5()

do test6()

do test7()

*)

do test1()


do System.Console.ReadLine() |> ignore