let solve_deg_3 (t,a2,a1,a0) =
  if t = 0.0 then raise Not_degre_3;
  let a2 = a2 /. t
  and a1 = a1 /. t
  and a0 = a0 /. t in
  let p = (-. 1. /. 3.) *. (a2 ** 2.0) +. a1
  and q = (2. /. 27.) *. (a2 ** 3.0) -. (1. /. 3.) *. a1 *. a2 +. a0 in
  let delta = 4. *. (p ** 3.0) +. 27. *. (q ** 2.0) in
  let j = Complex.polar 1.0 (2. *. pi /. 3.) in
  let j' = Complex.conj j in
  let check = function
      `Real x ->
        !print_check (Printf.sprintf " %f^3 + b/a %f^2 + c/a %f + d/a = %f"
                         x x x
                         ((x**3.0) +. a2 *. (x**2.0) +. a1 *. x +. a0))
    | `Complex c ->
        let x = string_of_complex c in
        let sum =
          List.fold_left
            Complex.add
            {Complex.re = a0; Complex.im = 0.0}
            [ (c **? 3.0) ;
              a2 *? (c **? 2.0);
              a1 *? c ;
            ]
        in
        !print_check (Printf.sprintf " %s^3 + b/a %s^2 + c/a %s + d/a = %s"
                         x x x (string_of_complex sum))
  in
  if delta >= 0.0 then
    (
     let u = ((-. 27. *. q +. 3. *. (sqrt 3.) *. (sqrt delta)) /. 2.0) ** (1. /. 3.) in
     let v0 = (-. 27. *. q -. 3. *. (sqrt 3.) *. (sqrt delta)) /. 2.0 in
     let v = if v0 >= 0.0 then v0 ** (1. /. 3.) else -. ((-. v0) ** (1. /. 3.)) in
     !print_check (Printf.sprintf "u=%f, v=%f" u v);
     let x1 = (1. /. 3.) *. (u +. v) in
     let x2 = (1. /. 3.) *? (Complex.add (u *? j) (v *? j'))
     and x3 = (1. /. 3.) *? (Complex.add (u *? j') (v *? j)) in
     let fc x = `Complex { x with Complex.re = x.Complex.re -. a2 /. 3. } in
     let x1 = `Real (x1 -. a2 /. 3.) in
     check x1;
     check (fc x2);
     check (fc x3);
     (x1, fc x2, fc x3)
    )
  else
    (
     let u0 =
       { Complex.re = (-. 27. *. q /. 2.0) ;
         Complex.im = (3. *. (sqrt 3.) *. (sqrt (-. delta)) /. 2.) ;
       }
     in
     let u = u0 **? (1. /. 3.) in
     !print_check (Printf.sprintf "u0=%f+%fi, u=%f+%fi"
                      u0.Complex.re u0.Complex.im
                      u.Complex.re u.Complex.im);
     let u' = Complex.conj u in
     let x1 = (Complex.add u u').Complex.re /. 3.
     and x2 = (Complex.add (Complex.mul j u) (Complex.mul j' u')).Complex.re /. 3.0
     and x3 =
       (Complex.add
          (Complex.mul (Complex.mul j j) u)
          (Complex.mul (Complex.mul j' j') u')).Complex.re /. 3.0
     in
     let f x = `Real (x -. a2 /. 3.) in
     check (f x1);
     check (f x2);
     check (f x3);
     (f x1, f x2, f x3)
    )