(** 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 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; ];;