(** Corrigé du TP5 *)

(** Exercice 1 *)

let rec randlist maxval n = match n with
    | 0 -> []
    | _ -> random__int maxval :: randlist maxval (n-1) ;;

let mesurer_temps f x =
    let n = 20 in  (* nombre d'iterations *)
    let start = sys__time () in
    for i = 1 to n do
        let _ = f x in  (* ignorer le résultat *)
        ()
    done;
    let stop = sys__time () in
    (stop -. start) /. (float_of_int n);;

let rec insert x l = match l with
    | [] -> [x]
        | y::l' -> if x <= y then x :: l else y :: insert x l';;

let rec insert_sort l = match l with
    | [] -> []
        | x :: l' -> insert x (insert_sort l');;

(* Voir le cours pour plus de détails *)

let rec merge l1 l2 = match l1, l2 with
    | [], _ -> l2
    | _, [] -> l1
    | x::l1', y::l2' ->
        if x < y then x :: merge l1' l2
        else if x > y then y :: merge l1 l2'
        else x :: y :: merge l1' l2'
;;

let split l =
    let rec split_aux acc1 acc2 l = match l with
        | [] -> acc1, acc2
        | x::l' -> split_aux acc2 (x::acc1) l'
    in split_aux [] [] l;;

let rec merge_sort l = match l with
    | [] -> []
    | [x] -> [x]
    | _ ->
        let l1, l2 = split l in
        let l1' = merge_sort l1 in
        let l2' = merge_sort l2 in
        merge l1' l2'
;;

(* comparaison des deux fonctions de tri sur une liste de n valeurs,
  ce qui simplifie l'appel suivant et permet de tracer un graphe des
  temps d'éxecution *)
let comparer n =
  let l = randlist n n in
  n, [ "insertion", mesurer_temps insert_sort l;
        "fusion", mesurer_temps merge_sort l ]
;;

(* graphe d'éxecution. Voir le fichier
  'http://who.rocq.inria.fr/Simon.Cruanes/enseignement/tris.png' *)
let comparaison_tris =
  map comparer [10; 20; 50; 100; 1000; 1500; 2000; 2500; 3000; ] ;;

(** Exercice 2 *)

(** 2.1 *)

let puissance_naive x n =
  let rec aux acc n = match n with
    | 0 -> acc
    | 1 -> x*acc
    | _ -> aux (x * acc) (n-1)
  in
  aux 1 n;;

puissance_naive 2 10 = 1024;;

let rec puissance x n = match n with
  | 0 -> 1
  | 1 -> x
  | n when n mod 2 = 0 ->
      let x2 = x*x in
      puissance x2 (n/2)
  | _ ->
      let x2 = x*x in
      x * (puissance x2 (n/2))
;;

puissance 2 10 = 1024;;
puissance 2 15 = puissance_naive 2 15;;

(** 2.2 *)

(* Ici nous définissons un type abstrait de "monoide", une structure
  simple munie d'un produit associatif et d'un élément neutre *)

type 'a monoide = {
    neutre : 'a;
    produit : ('a -> 'a -> 'a);
};;

let monoide_int = {
    neutre = 1;
    produit = (fun x y -> x*y);
};;

(* Note: l'addition forme également un monoïde sur les entiers,
   d'élément neutre 0. En ce cas, l'exponentiation rapide sert
   à calculer le produit ! *)

let monoide_string = {
    neutre = "";
    produit = (fun x y -> x^y);  (* concaténation *)
};;


let puissance_naive2 monoide x n =
  let rec aux acc n = match n with
    | 0 -> acc
    | 1 -> monoide.produit x acc
    | _ -> aux (monoide.produit x acc) (n-1)
  in
  aux monoide.neutre n;;

(* "a" puissance 1000 est la chaîne composée de 1000 caractères `a` *)

puissance_naive2 monoide_string "a" 1000;;

let rec puissance2 monoide x n = match n with
  | 0 -> monoide.neutre
  | 1 -> x
  | _ when n mod 2 = 0 ->
      let x2 = monoide.produit x x in
      puissance2 monoide x2 (n/2)
  | _ ->
      let x2 = monoide.produit x x in
      monoide.produit x (puissance2 monoide x2 (n/2));;

(* version récursive terminale *)

let puissance3 monoide x n =
  (* acc: le produit des x *)
  let rec aux acc x n = match n with
    | 0 -> acc
    | 1 -> monoide.produit x acc
    | _ ->
        let acc' = (if n mod 2 = 0 then acc else monoide.produit x acc) in
        let x2 = monoide.produit x x in
        aux acc' x2 (n/2)
  in
  aux monoide.neutre x n;;

(** 2.3 *)

type matrix == int vect vect ;;

let produit_matrice m1 m2 =
  let n_lignes = vect_length m1 in
  let n_colonnes = vect_length m2.(0) in
  (* le nbr de colonnes de m1 doit correspondre
   * au nombre de lignes de m2 *)
  if vect_length m1.(0) <> vect_length m2
    then failwith "mauvaises dimensions";
  (* calcule la case (i,j) du résultat *)
  let case i j =
    let n = ref 0 in
    for k=0 to vect_length m2 - 1 do
      n := !n +  m1.(i).(k) * m2.(k).(j)
    done;
    !n
  in
  (* créer le résultat à la volée. *)
  let m3 = init_vect n_lignes
    (fun i ->
        init_vect n_colonnes
          (fun j -> case i j))
  in
  m3
;;

(* monoide pour les matrices 2x2 (impossible hélas de définir
  ce monoïde génériquement pour les matrices nxn car l'élément
  neutre devrait être arbitrairement grand ! *)
let monoide_matrix = {
  produit = produit_matrice;
  neutre = [| [| 1; 0 |]; [| 0; 1 |] |];
} ;;

(* exponentiation rapide sur cette matrice, ce qui donne
   une fonction efficace. *)
let fib n =
  let m = [| [| 1; 1 |]; [| 1; 0 |] |] in
  let m_n = puissance3 monoide_matrix m n in
  m_n.(0).(1)
;;

(** Exercice 3 *)

let base = 1024 ;;

let rec produit_int i j = match () with
  | _ when i < base && j < base -> i*j
  | _ ->
      let x1, x2 = i/base, i mod base in
      let y1, y2 = j/base, j mod base in
      let z1 = produit_int x1 y1 in
      let z2 = produit_int x2 y2 in
      let z3 = (produit_int (x1 + x2) (y1 + y2)) - z1 - z2 in
      z1 * (base*base) + z3 * base + z2;;

let degre p = vect_length p - 1 ;;

let somme_poly p q =
  let m = vect_length p in
  let n = vect_length q in
  (* lequel, de p et q, a le plus haut degré ? *)
  let bigger = if m > n then p else q in
  let r = make_vect (max m n) 0 in
  for i = 0 to (min m n)-1 do
    r.(i) <- p.(i) + q.(i)
  done;
  (* compléter le résultat si besoin, en ajoutant
    les coefficients du polynôme de degré maximal *)
  for i = min m n to (max m n) - 1 do
    r.(i) <- bigger.(i)
  done;
  r ;;

(* notation préfixe pour l'addition de polynômes *)
let prefix ++ p q = somme_poly p q;;

(* opposé *)
let opp_poly p =
  init_vect (vect_length p) (fun i -> - p.(i));;

let moins_poly p q = somme_poly p (opp_poly q);;

let prefix -- p q = moins_poly p q ;;

(* produit de polynômes simple. En pratique il est assez
  efficace et il est difficile (et compliqué) de faire mieux. *)
let prod_poly p q =
  let r = make_vect (degre p + degre q + 1) 0 in
  for i = 0 to degre p do
    for j = 0 to degre q do
      r.(i+j) <- r.(i+j) + p.(i) * q.(j);
    done;
  done;
  r;;

let prefix *** p q = prod_poly p q;;

(* quotient de la division euclidienne de p par X^m *)
let quotient_poly p m =
  let n = degre p in
  if n<m
  then [| 0 |]
  else init_vect (n-m+1) (fun i -> p.(i+m)) ;;

(* reste de la division euclidienne de p par X^m *)
let reste_poly p m =
  let degre_resultat = ref (min (degre p) (m-1)) in
  while !degre_resultat >= 0 &&  p.(!degre_resultat) = 0 do
    decr degre_resultat;
  done;
  let r = init_vect (!degre_resultat + 1) (fun i -> p.(i)) in
  r ;;

(* produit de p par X^m (décalage des coefficients
   de m cases vers la droite) *)
let shift_poly p m =
  let r = make_vect (m + vect_length p) 0 in
  for i = 0 to degre p do
    r.(i+m) <- p.(i);
  done;
  r;;

(* le produit se fait sur des "blocs" de degré 50 *)
let base_poly = 50;;

let rec karatsuba_poly p q =
  if degre p < base_poly && degre q < base_poly
  then p *** q
  else
    let p1 = quotient_poly p base_poly
    and p2 = reste_poly p base_poly
    and q1 = quotient_poly q base_poly
    and q2 = reste_poly q base_poly in
    (* recursion *)
    let z1 = karatsuba_poly p1 q1 in
    let z2 = p2 *** q2 in
    let z3 =
      ((karatsuba_poly (p1 ++ p2) (q1 ++ q2)) -- z1) -- z2 in
    (* reconstruction *)
    (shift_poly z1 (2*base_poly)) ++ (shift_poly z3 base_poly ++ z2) ;;

(* les polynômes aussi forment un monoïde !
   l'élément neutre est le polynôme constant "1" *)

let monoide_poly = {
  produit = prod_poly;
  neutre = [| 1 |];
} ;;

let exp_poly p n =
  puissance3 monoide_poly p n;;

exp_poly [| 1; 1 |] 10;;
exp_poly [| 1; 1 |] 20;;

(* admirer le triangle de Pascal dans toute sa splendeur. *)
let pascal_20 = init_vect 20 (exp_poly [| 1; 1 |]);;

(** Test *)

(* créer une string représentant le polynome *)
let affiche_poly v =
  let s = ref "[|" in
  let first = ref true in
  for i = 0 to vect_length v - 1 do
    (if !first then first := false else s := !s ^ "; ");
    s := !s ^ string_of_int v.(i)
  done;
  !s ^ "|]"
;;

(* afficher le triangle de Pascal plus lisiblement.
  "do_vect f a" appelle "f" sur chaque élément de "a". *)
do_vect
  (fun p -> print_endline (affiche_poly p))
  pascal_20;;

(* vérifier si p=q, en ignorant les 0 du début si l'un est plus long
 * que l'autre. *)
let egal_poly p q =
  let res = ref true in
  let m = degre p in
  let n = degre q in
  for i = 0 to min m n do
    if p.(i) <> q.(i) then res := false
  done;
  let bigger = if m > n then p else q in
  for i = min m n + 1 to max m n do
    if bigger.(i) <> 0 then res := false
  done;
  !res ;;

(* polynome aleatoire dont les coeffs sont entre -20 et 20 *)
let random_poly degre =
    let p = init_vect degre (fun _ -> random__int 41 - 20) in
    (* il faut etre sur que le dernier coefficient n'est pas 0 *)
    p.(vect_length p - 1) <- random__int 12 + 1;
    p;;

(* liste de n paires de polynomes aleatoires *)
let rec random_paires n = match n with
    | 0 -> []
    | _ ->
    (random_poly 30, random_poly 30) :: random_paires (n-1);;

(* tester que f1 et f2 se comportent pareillement sur des polynomes
   aleatoires *)
let tester_similaire f1 f2 =
    let l = random_paires 100 in
    (* vérifier, sur chaque paire de l, que f1 et f2 renvoient
     * le meme résultat *)
    list__for_all
        (fun (p1,p2) ->
          let ok = egal_poly (f1 p1 p2) (f2 p1 p2) in
          if not ok
            then print_endline
              ("probleme pour: " ^ affiche_poly p1 ^ " et " ^ affiche_poly p2);
          ok
        ) l;;

tester_similaire prod_poly karatsuba_poly;;
(* tester_similaire prod_poly karatsuba_bis;; *)

(* produit de polynômes: mesurer le temps. Ici on
voit que le produit de Karatsuba n'est pas si bon en
pratique à cause de l'implémentation (nombreuses copies
et sommes) *)
let uncurry f (x,y) = f x y;;

let temps_produit = map
  (fun n ->
    let p = random_poly n in
    let q = random_poly n in
    n, ["prod_poly", mesurer_temps (uncurry prod_poly) (p,q);
        "karatsuba", mesurer_temps (uncurry karatsuba_poly) (p,q);
      ])
  [10; 50; 80; 100; 150; 200; 300;  ];;