(* param 1 = size of the matrices at the beginning
   param 2 = incrementation
   param 3 = size of the matrices at the end
   param 4 = number of iteration 

ocamlopt -unsafe -o ocamlmult mult.ml

*)

let multiplication_polymorphic n zero add mul a b c =
  for i=0 to n-1 do
   for j=0 to n-1 do
    let rec loop tmp k =
     if k=n then tmp else loop (add tmp (mul b.(k*n+j) a.(i*n+k))) (k+1)
    in 
     c.(i*n+j)<-loop zero 0
   done
  done

let multiplication_monomorphic n a b c =
 let tmp=ref Complex.zero in
  for i=0 to n-1 do
   for j=0 to n-1 do
    for k=0 to n-1 do
      tmp:=Complex.add !tmp (Complex.mul b.(k*n+j) a.(i*n+k))
    done;
    c.(i*n+j)<-(!tmp)
   done
  done
(* same efficienty with the loop function *)

let random_complex n = {Complex.re=Random.float n; im=Random.float n}

let matrixRandom n = Array.init (n*n) (fun _ -> random_complex 1000.)

let _ = 
	Random.self_init();
        let i = ref (int_of_string (Sys.argv).(1))
	and pas=(int_of_string (Sys.argv).(2))
	and maxi=(int_of_string (Sys.argv).(3)) 
	and iter=(int_of_string (Sys.argv).(4)) in
	while !i<maxi do
         let a=matrixRandom (!i)
         and b=matrixRandom (!i)
         and c=Array.make ((!i)*(!i)) Complex.zero in
	  for j=1 to iter do
           (* ignore(multiplication_polymorphic (!i) Complex.zero Complex.add Complex.mul a b c); *)
	   ignore(multiplication_monomorphic (!i) a b c);
          done;
	 i:=!i+pas;
	done