let solve_deg_4 ?(imeps=0.000001) ?(compeps=0.00000000001) (t,a1,b1,c1,d1) =
  !print_check
    (Printf.sprintf "solving %fX^4 + %fX^3 + %fX^2 + %fX + %f = 0"
                       t a1 b1 c1 d1);
  let a1 = a1 /. t
  and b1 = b1 /. t
  and c1 = c1 /. t
  and d1 = d1 /. t in
  let b = (-. 3. /. 8.) *. (a1 ** 2.0) +. b1
  and c = ((a1 ** 3.0) /. 8.) -. (a1 *. b1 /. 2.) +. c1
  and d = (-. 3. /. 256.) *. (a1 ** 4.0) +.
      (b1 *. (a1 ** 2.0) /. 16.) -. (a1 *. c1 /. 4.) +. d1 in
  let real_or_complex c =
    if c.Complex.im >= -. imeps && c.Complex.im <= imeps then
      `Real c.Complex.re
    else
      `Complex c
  in
  let to_sol = function
      `Real z ->
        let x = z -. a1 /. 4. in
        !print_check
          (Printf.sprintf "%f^4 + a1 %f^3 + b1 %f^2 + c1 %f + d1 = %f"
             x x x x
             ((x**4.0) +. a1 *. (x**3.0) +. b1 *. (x**2.0) +. c1 *. x +. d1)
          );
        `Real x
    | `Complex z ->
        let c = { z with Complex.re = z.Complex.re -. a1 /. 4. } in
        let s = string_of_complex c in
        let sum =
          List.fold_left
            Complex.add
            {Complex.re = d1; Complex.im = 0.0}
            [ (c **? 4.0) ;
              a1 *? (c **? 3.0);
              b1 *? (c **? 2.0);
              c1 *? c ;
            ]
        in
        !print_check (Printf.sprintf " %s^4 + b/a %s^3 + c/a %s^2 + d/a %s + e/a = %s"
                         s s s s (string_of_complex sum));

        `Complex c
  in
  let to_sol = function
      `Real x -> to_sol (`Real x)
    | `Complex c -> to_sol (real_or_complex c)
  in
  if d >= -. compeps && d <= compeps then
    let (z1,z2,z3) = solve_deg_3 (1.0, 0.0, b, c) in
    (to_sol (`Real 0.0), to_sol z1, to_sol z2, to_sol z3)
  else
    if c >= -. compeps && c <= compeps then
      let (z1,z2) = solve_deg_2 (1.0, b, d) in
      let to_sol2 = function
          `Real x -> to_sol (`Real (sqrt x))
        | `Complex _ -> `Real nan (* A VOIR *)
      in
      let to_sol2' = function
          `Real x -> to_sol (`Real (-. (sqrt x)))
        | `Complex _ -> `Real nan (* A VOIR *)
      in
      (to_sol2 z1, to_sol2' z1, to_sol2 z2, to_sol2' z2)
    else
      begin
        let (g1,g2,g3) =
          match solve_deg_3 (1.0, 8. *. b, 16. *.((b**2.0) -. 4. *. d), -. (64.0 *. (c**2.0))) with
            (`Real g1, `Real g2, `Real g3) ->
              ({ Complex.re = g1 ; im = 0.0 },
               { Complex.re = g2 ; im = 0.0 },
               { Complex.re = g3 ; im = 0.0 })
          | (`Real g1, `Complex g2, `Complex g3) ->
              ({ Complex.re = g1 ; im = 0.0 }, g2, g3)
          | _ -> assert false
        in
        let ro1 = Complex.sqrt g1 in
        let ro2 = Complex.sqrt g2 in
        let ro3 =
          let r = { Complex.re = (-. 8.) *. c; im = 0.0 } in
          prerr_endline (Printf.sprintf "-8c = %s" (string_of_complex r));
          Complex.div r (Complex.mul ro1 ro2)
        in
        let g3' = ro3 **? 2.0 in
        !print_check (Printf.sprintf "g3'=%s =? %s=g3"
                        (string_of_complex g3')
                      (string_of_complex g3)
                     );
(*
   let ro3 =
   if c >= 0.0 then
   Complex.neg (g3 **? (1. /. 2.))
   else
   g3 **? (1. /. 2.)
   in
*)
        let foo = Complex.mul ro1 (Complex.mul ro2 ro3) in
        !print_check
          (Printf.sprintf "ro1.ro2.ro3 = %s\n-8c = %f"
           (string_of_complex foo)
             (-. 8. *. c)
          );
        let (+?) = Complex.add and (-?) = Complex.sub in
        let q = 1. /. 4. in
        let x1 = to_sol (`Complex (q *? (ro1 +? ro2 +? ro3)))
        and x2 = to_sol (`Complex (q *? ((ro1 -? ro2) -? ro3)))
        and x3 = to_sol (`Complex (q *? ((ro2 -? ro3) -? ro1)))
        and x4 = to_sol (`Complex (q *? ((ro3 -? ro1) -? ro2))) in
        (x1, x2, x3, x4)
      end