(* Matrices creuses *)

module M = Map.Make(struct type t = int let compare = compare end)

type t = { cols: int; rows: int M.t array }

let id n =
  { cols = n; rows = Array.init n (fun i -> M.singleton i 1) }

let init n m f =
  let rec row acc i j =
    if j = m then acc else
    row (let v = f i j in if v = 0 then acc else M.add j v acc) i (j + 1) in
  { cols = m; rows = Array.init n (fun i -> row M.empty i 0) }

let size a =
  (Array.length a.rows, a.cols)

let add a b =
  if size a <> size b then invalid_arg "add";
  let add j aij bij = match aij, bij with
    | Some x, Some y -> let z = x + y in if z = 0 then None else Some z
    | (Some _ as x), None | None, (Some _ as x) -> x
    | None, None -> None
  in
  { cols = a.cols;
    rows = Array.mapi (fun i r -> M.merge add r b.rows.(i)) a.rows }

let mul a b =
  let n, p = size a in
  let q, m = size b in
  if q <> p then invalid_arg "mul";
  let rows = Array.make n M.empty in
  let add i j v =
    let cij = v + try M.find j rows.(i) with Not_found -> 0 in
    rows.(i) <- if cij = 0 then M.remove j rows.(i) else M.add j cij rows.(i)
  in
  for i = 0 to n - 1 do
    M.iter
      (fun k aik ->
        M.iter
          (fun j bkj -> add i j (aik * bkj))
          b.rows.(k))
      a.rows.(i)
  done;
  { cols = m; rows = rows }

let rec power x n =
  if n = 0 then
    id x.cols
  else
    let y = power x (n/2) in
    if n mod 2 = 1 then mul x (mul y y) else mul y y

let rec print fmt m =
  Array.iter
    (fun rowi ->
      for j = 0 to m.cols - 1 do
        Format.fprintf fmt "%2d " (try M.find j rowi with Not_found -> 0)
      done;
      Format.fprintf fmt "@\n")
    m.rows

(* tests *)

let rec forall i j p = i >= j || p i && forall (i + 1) j p

let equal m1 m2 =
  m1.cols = m2.cols &&
  Array.length m1.rows = Array.length m2.rows &&
  forall 0 (Array.length m1.rows) (fun i -> M.equal (=) m1.rows.(i) m2.rows.(i))

let () =
  let m = id 17 in
  assert (equal (mul m m) m);
  let m' = init 17 17 (fun i j -> Random.int 100) in
  assert (equal (mul m m') m');
  assert (equal (mul m' m) m')

(* on calcule Fib(n) (voir exercices 10.11 et 10.12) avec une matrice
   inutilement grande (1000x1000) *)

let fib n =
  let size = 1000 in
  let mfib = init size size
    (fun i j -> if i = 0 && j <= 1 || i = 1 && j = 0 then 1 else 0) in
  try M.find 1 (power mfib n).rows.(0) with Not_found -> 0

let () =
  assert (fib 0 = 0);
  assert (fib 10 = 55);
  assert (fib 14 = 377)

This document was generated using caml2html