let rec root_bounding_interval = function
  (p, nbrp, prec, i) -> if gt_rat (prec, width (i))
                          then i
                          else
                            if eq_rat (eval_poly (p, midpoint (i)), cons_rat ("0", "1"))
                              then cons_single (midpoint (i))
                              else
                                let
                                  (i1, i2) = cut_in_half (i)
                                  (* i1 and i2 are non-empty open intervals *)
                                  (* since the midpoint of i is not a root of p,
                                     none of the bounds of i1 and i2 is a root of p
                                  *)
                                in
                                  if nbrp (i1) = 1
                                    then (* the number of real roots of p in i1 is 1 *)
                                         root_bounding_interval (p, nbrp, prec, i1)
                                    else (* the number of real roots of p in i2 is 1 *)
                                          root_bounding_interval (p, nbrp, prec, i2)