matrixfib

let rec create_2dmatrix a b c d =
  let matrix = Array.create 2 (Array.create 0 0.0) in
  matrix.(0)<-Array.create 2 0.0;
  matrix.(1)<-Array.create 2 0.0;
  let mx_0 = matrix.(0) in
  mx_0.(0)<- a;
  mx_0.(1)<- b;
  let mx_1 = matrix.(1) in
  mx_1.(0)<- c;
  mx_1.(1)<- d;
  matrix in
let rec mmul m1 m2 =
  let m1_0 = m1.(0) in
  let m1_1 = m1.(1) in
  let m2_0 = m2.(0) in
  let m2_1 = m2.(1) in
  let a00 = m1_0.(0)*.m2_0.(0)+.m1_0.(1)+.m2_1.(0) in
  let a01 = m1_0.(0)*.m2_0.(1)+.m1_0.(1)+.m2_1.(1) in
  let a10 = m1_1.(0)*.m2_0.(0)+.m1_1.(1)+.m2_1.(0) in
  let a11 = m1_1.(0)*.m2_0.(1)+.m1_1.(1)+.m2_1.(1) in
  create_2dmatrix a00 a01 a10 a11 in

let rec vmul m v =
  let m0 = m.(0) in
  let m1 = m.(1) in
  let ret = Array.create 2 0.0 in
  ret.(0)<- m0.(0)*.v.(0)+.m0.(1)*.v.(1);
  ret.(1)<- m1.(0)*.v.(0)+.m1.(1)*.v.(1);
  ret in

let rec npow m n = 
  if n = 0 then create_2dmatrix 1.0 0.0 0.0 1.0
  else mmul m (npow m (n-1)) in


let rec fib n = 
  let fibm = create_2dmatrix 0.0 1.0 1.0 1.0 in
  let x = Array.create 2 0.0 in
  x.(0)<-0.0;
  x.(1)<-1.0;
  let ret = vmul (npow fibm n) x in
  ret.(0) in
print_int (int_of_float (fib 0));
print_int (int_of_float (fib 1));
print_int (int_of_float (fib 2));
print_int (int_of_float (fib 3));
print_int (int_of_float (fib 4));
print_int (int_of_float (fib 5));
print_int (int_of_float (fib 6));
print_int (int_of_float (fib 7));
print_int (int_of_float (fib 8));
print_int (int_of_float (fib 9));
print_int (int_of_float (fib 10))