let abs x =
  if x < 0 then -x else x
;;

let absf x =
  if x < 0. then -. x else x
;;

(** Somme des entiers de 1 à n *)
let somme n =
  let r = ref 0 in
  for i = 1 to n do
    (* invariant: !r + somme_{k=i...n} k = somme_{k=1...n} k *)
    r:= !r + i;
  done;
  !r
;;

somme 10;;

(** Factorielle *)
let fact n =
  let r = ref 1 in
  for i = 1 to n do
    (* invariant: !r * produit_{k=i...n} k = fact n
     * variant: i croit strictement et est inférieur à n *)
    r := !r * i;
  done;
  !r
;;

fact 5 ;;

(** Puissance entière *)
let puissance x n =
  let r = ref 1 in
  for i = 1 to n do
    (* invariant: !r = x^i *)
    r := x * !r;
  done;
  !r
;;

puissance 2 10;;  (* 1024 *)
puissance 2 0;; (* 1 *)

(** Racine carrée (approximation) *)
let racine x epsilon =
  let s = ref 1. in
  while absf (!s -. (x /. !s)) >= epsilon do
    (* invariant: difficile à exprimer, l'algorithme est approximé
        (ceci dit, on sait que !s > 0 comme il convient à une racine);
       variant: la distance entre !s et x décroît *)
    s := (!s +. (x /. !s)) /. 2. ;
  done;
  !s
;;

racine 2. 0.01;;
racine 2. 1e-5;;
racine 16. 1e-5;;

for i = 32 to 127 do
  print_char (char_of_int i);
  if (i-32) mod 32 = 0
    then print_newline();
done
;;


(** "Miroir" d'un nombre en base 10 *)
let miroir n0 =
  let n = ref n0 in
  let result = ref 0 in
  while !n <> 0 do
    (* invariant:
        let string = string_of_int in
        string(!result) ^ string(miroir !n) = string(miroir n0)
       variant: !n décroit strictement mais reste positif *)
    let q, r = !n/10, !n mod 10 in
    result := 10 * !result + r;
    n := q;
  done;
  !result
;;

miroir 12345678;;

(** Nombre parfait: somme des ses diviseurs stricts *)
let parfait n =
  let sum = ref 0 in
  for i = 1 to n-1 do
    (* invariant: !sum = somme des diviseurs de n qui sont < i.
        note: on pourrait s'arrêter à la racine entière de n plutôt qu'à n-1. *)
    if n mod i = 0
    then sum := !sum + i;
  done;
  !sum = n
;;

for i = 1 to 1000 do
  if parfait i then begin
    print_int i;
    print_newline();
  end;
done;;

let solve (a,b,c) =
  let epsilon = 1e-4 in
  let delta = b ** 2. -. 4. *. a *. c in
  if delta >= 0.
  then (((racine delta epsilon) -. b) /. (2. *. a),
		 (-. (racine delta epsilon) -. b) /. (2. *. a))
  else failwith "pas de solutions"
;;

solve (1., 4., 1.);;

(** Élément maximal d'un tableau d'entiers *)
let maximum a =
  if vect_length a = 0
  then failwith "tableau vide";
  let r = ref a.(0) in
  for i = 1 to vect_length a - 1 do
    (* invariant: !r = élément maximal parmis a.(0.... i-1) *)
    r := max !r a.(i);
  done;
  !r
;;

maximum [| 1; 3; 11; 0; -5; 10; 200; 0 |];;

let affiche v =
  let first = ref true in
  for i = vect_length v - 1 downto 0 do
    if v.(i) <> 0 then begin
      (* doit-on afficher un "+" ? Uniquement si on a déjà affiché un
          coefficient *)
      if not !first then print_string " + ";
      first := false;
      print_int v.(i);
      if i > 0 then begin
        print_string " * X^";
        print_int i;
      end
    end
  done
;;

affiche [| 1; 2; 0; 4; 0; -2 |];;

(** Polynomes. Notons p[0...i] le polynome de degré au plus i, dont les
   coefficients sont les mêmes que ceux de p pour X^k avec 0<=k<= i.
   En particulier p[i] est le coefficient de X^i dans p *)

(** Évaluation d'un polynôme *)
let eval p x =
  let r = ref 0. in
  for i = 0 to vect_length p - 1 do
    (* invariant: !r = eval p[0...i-1] *)
    r := !r +. (float_of_int p.(i)) *. (x ** (float_of_int i));
  done;
  !r
;;

(** Somme de deux polynômes *)
let somme p1 p2 =
  let p1, p2 = if vect_length p1 > vect_length p2 then p2, p1 else p1, p2 in
  (* p2: polynome de plus grand degre *)
  let n = vect_length p2 in
  let p = copy_vect p2 in
  for i = 0 to vect_length p1 - 1 do
    (* invariant: p = p1[0... i-1] + p2[0... i-1] (de degré i-1 donc) *)
    p.(i) <- p.(i) + p1.(i)
  done;
  p
;;

(** Produit de deux polynômes. *)
let produit p1 p2 =
  let n = vect_length p1 + vect_length p2 - 1 in
  let p = make_vect n 0 in
  for i = 0 to vect_length p1 - 1 do
    (* invariant: p le produit de p1[0... i-1] et p2 *)
    for j = 0 to vect_length p2 - 1 do
      (* invariant: p est   p1[0...i-1] * p2 + p1[i] * p2[0... j-1] *)
      p.(i+j) <- p.(i+j) + p1.(i) * p2.(j)
    done;
  done;
  p
;;