Theory Bernstein_01

section Bernstein Polynomials over the interval [0, 1]

theory Bernstein_01
  imports "HOL-Computational_Algebra.Computational_Algebra" 
    "Budan_Fourier.Budan_Fourier"
    "RRI_Misc"
begin

text 
The theorem of three circles is a statement about the Bernstein coefficients of a polynomial, the
coefficients when a polynomial is expressed as a sum of Bernstein polynomials. These coefficients
behave nicely under translations and rescaling and are the coefficients of a particular polynomial
in the [0, 1] case. We shall define the [0, 1] case now and consider the general case later,
deriving all the results by rescaling.


subsection Definition and basic results

definition Bernstein_Poly_01 :: "nat  nat  real poly" where
  "Bernstein_Poly_01 j p = (monom (p choose j) j) 
                              * (monom 1 (p-j) p [:1, -1:])"

lemma degree_Bernstein: 
  assumes hb: "j  p" 
  shows "degree (Bernstein_Poly_01 j p) = p"
proof -
  have ha: "monom (p choose j) j  (0::real poly)" using hb by force
  have hb: "monom 1 (p-j) p [:1, -1:]  (0::real poly)"
  proof
    assume "monom 1 (p-j) p [:1, -1:] = (0::real poly)"
    hence "lead_coeff (monom 1 (p - j) p [:1, -1:]) = (0::real)"
      apply (subst leading_coeff_0_iff)
      by simp
    moreover have "lead_coeff (monom (1::real) (p - j) 
        p [:1, -1:]) = (((- 1) ^ (p - j))::real)"
      by (subst lead_coeff_comp, auto simp: degree_monom_eq)
    ultimately show "False" by auto
  qed
  from ha hb show ?thesis
    by (auto simp add: Bernstein_Poly_01_def degree_mult_eq 
          degree_monom_eq degree_pcompose)
qed

lemma coeff_gt: 
  assumes hb: "j > p" 
  shows "Bernstein_Poly_01 j p = 0"
  by (simp add: hb Bernstein_Poly_01_def)

lemma degree_Bernstein_le: "degree (Bernstein_Poly_01 j p)  p"
  apply (cases "j  p")
  by (simp_all add: degree_Bernstein coeff_gt)

lemma poly_Bernstein_nonneg: 
  assumes "x  0" and "1  x" 
  shows "poly (Bernstein_Poly_01 j p) x  0"
  using assms by (simp add: poly_monom poly_pcompose Bernstein_Poly_01_def)

lemma Bernstein_symmetry: 
  assumes "j  p"
  shows "(Bernstein_Poly_01 j p) p [:1, -1:] = Bernstein_Poly_01 (p-j) p"
proof -
  have "(Bernstein_Poly_01 j p) p [:1, -1:]
         = ((monom (p choose j) j) * (monom 1 (p-j) p [:1, -1:])) p [:1, -1:]"
    by (simp add: Bernstein_Poly_01_def)
  also have "... = (monom (p choose (p-j)) j * 
                    (monom 1 (p-j) p [:1, -1:])) p [:1, -1:]" 
    by (fastforce simp: binomial_symmetric[OF assms])
  also have "... = monom (p choose (p-j)) j p [:1, -1:] * 
                   (monom 1 (p-j)) p ([:1, -1:] p [:1, -1:])"
    by (force simp: pcompose_mult pcompose_assoc)
  also have "... = (monom (p choose (p-j)) j p [:1, -1:]) * monom 1 (p-j)"
    by (force simp: pcompose_pCons)
  also have "... = smult (p choose (p-j)) (monom 1 j p [:1, -1:]) 
                    * monom 1 (p-j)" 
    by (simp add: assms smult_monom pcompose_smult[symmetric])
  also have "... = (monom 1 j p [:1, -1:]) * monom (p choose (p-j)) (p-j)"
    apply (subst mult_smult_left)
    apply (subst mult_smult_right[symmetric])
    apply (subst smult_monom)
    by force
  also have "... = Bernstein_Poly_01 (p-j) p" using assms
    by (auto simp: Bernstein_Poly_01_def)
  finally show ?thesis .
qed

subsection @{term Bernstein_Poly_01} and @{term reciprocal_poly}

lemma Bernstein_reciprocal: 
  "reciprocal_poly p (Bernstein_Poly_01 i p) 
    = smult (p choose i) ([:-1, 1:]^(p-i))"
proof cases
  assume "i  p"
  hence "reciprocal_poly p (Bernstein_Poly_01 i p) = 
         reciprocal_poly (degree (Bernstein_Poly_01 i p)) (Bernstein_Poly_01 i p)"
    by (auto simp: degree_Bernstein)
  also have "... = reflect_poly (Bernstein_Poly_01 i p)"
    by (rule reciprocal_degree)
  also have "... = smult (p choose i) ([:-1, 1:]^(p-i))"
    by (auto simp: Bernstein_Poly_01_def reflect_poly_simps monom_altdef
         pcompose_pCons reflect_poly_pCons' hom_distribs)
  finally show ?thesis .
next
  assume h:"¬ i  p"
  hence "reciprocal_poly p (Bernstein_Poly_01 i p) = (0::real poly)"
    by (auto simp: coeff_gt reciprocal_poly_def)
  also have "... = smult (p choose i) ([:-1, 1:]^(p - i))" using h
    by fastforce
  finally show ?thesis .
qed

lemma Bernstein_reciprocal_translate: 
  "reciprocal_poly p (Bernstein_Poly_01 i p) p [:1, 1:] = 
   monom (p choose i) (p - i)"
  by (auto simp: Bernstein_reciprocal pcompose_smult pcompose_pCons monom_altdef hom_distribs)

lemma coeff_Bernstein_sum_01: fixes b::"nat  real" assumes hi: "p  i"
  shows 
    "coeff (reciprocal_poly p 
            (x = 0..p. smult (b x) (Bernstein_Poly_01 x p)) p [:1, 1:]) 
      (p - i) = (p choose i) * (b i)" (is "?L = ?R")
proof -
  define P where "P  (x = 0..p. (smult (b x) (Bernstein_Poly_01 x p)))"

  have "x. degree (smult (b x) (Bernstein_Poly_01 x p))  p"
  proof -
    fix x
    show "degree (smult (b x) (Bernstein_Poly_01 x p))  p"
      apply (cases "x  p")
      by (auto simp: degree_Bernstein coeff_gt)
  qed
  hence "reciprocal_poly p P = 
         (x = 0..p. reciprocal_poly p (smult (b x) (Bernstein_Poly_01 x p)))"
    apply (subst P_def)
    apply (rule reciprocal_sum)
    by presburger
  also have
    "... = (x = 0..p. (smult (b x * (p choose x)) ([:-1, 1:]^(p-x))))"
  proof (rule sum.cong)
    fix x assume "x  {0..p}"
    hence "x  p" by simp
    thus "reciprocal_poly p (smult (b x) (Bernstein_Poly_01 x p)) =
          smult ((b x) * (p choose x)) ([:-1, 1:]^(p-x))"
      by (auto simp add: reciprocal_smult degree_Bernstein Bernstein_reciprocal)
  qed (simp)
  finally have 
    "reciprocal_poly p P = 
     (x = 0..p. (smult ((b x) * (p choose x)) ([:-1, 1:]^(p-x))))" .
  hence 
    "(reciprocal_poly p P) p [:1, 1:] = 
     (x = 0..p. (smult ((b x) * (p choose x)) ([:-1, 1:]^(p-x))) p [:1, 1:])"
    by (simp add: pcompose_sum pcompose_add)
  also have "... = (x = 0..p. (monom ((b x) * (p choose x)) (p - x)))"
  proof (rule sum.cong)
    fix x assume "x  {0..p}"
    hence "x  p" by simp
    thus "smult (b x * (p choose x)) ([:- 1, 1:] ^ (p - x)) p [:1, 1:] =
          monom (b x * (p choose x)) (p - x)"
      by (simp add: hom_distribs pcompose_smult pcompose_pCons monom_altdef)
  qed (simp)
  finally have "(reciprocal_poly p P) p [:1, 1:] = 
                (x = 0..p. (monom ((b x) * (p choose x)) (p - x)))" .
  hence "?L = (x = 0..p. if p - x = p - i then b x * real (p choose x) else 0)"
    by (auto simp add: P_def coeff_sum)
  also have "... = (x = 0..p. if x = i then b x * real (p choose x) else 0)"
  proof (rule sum.cong)
    fix x assume "x  {0..p}"
    hence "x  p" by simp
    thus "(if p - x = p - i then b x * real (p choose x) else 0) =
          (if x = i then b x * real (p choose x) else 0)" using hi
      by (auto simp add: leI)
  qed (simp)
  also have "... = ?R" by simp
  finally show ?thesis .
qed

lemma Bernstein_sum_01: assumes hP: "degree P  p"
  shows 
  "P = (j = 0..p. smult 
     (inverse (real (p choose j)) * 
      coeff (reciprocal_poly p P p [:1, 1:]) (p-j))
   (Bernstein_Poly_01 j p))"
proof -
  define Q where "Q  reciprocal_poly p P p [:1, 1:]"
  from hP Q_def have hQ: "degree Q  p" 
    by (auto simp: degree_reciprocal degree_pcompose)
  have "reciprocal_poly p (j = 0..p. 
        smult (inverse (real (p choose j)) * coeff Q (p-j)) 
        (Bernstein_Poly_01 j p)) p [:1, 1:] = Q"
  proof (rule poly_eqI)
    fix n
    show "coeff (reciprocal_poly p (j = 0..p. 
          smult (inverse (real (p choose j)) * coeff Q (p-j))
          (Bernstein_Poly_01 j p)) p [:1, 1:]) n = coeff Q n" 
      (is "?L = ?R")
    proof cases
      assume hn: "n  p"
      hence "?L = coeff (reciprocal_poly p (j = 0..p. 
             smult (inverse (real (p choose j)) * coeff Q (p-j)) 
             (Bernstein_Poly_01 j p)) p [:1, 1:]) (p - (p - n))"
        by force
      also have "... = (p choose (p-n)) * 
                       (inverse (real (p choose (p-n))) * 
                        coeff Q (p-(p-n)))"
        apply (subst coeff_Bernstein_sum_01)
        by auto
      also have "... = ?R" using hn
        by fastforce
      finally show "?L = ?R" .
    next
      assume hn: "¬ n  p"
      have "degree (j = 0..p.
            smult (inverse (real (p choose j)) * coeff Q (p - j))
            (Bernstein_Poly_01 j p))  p"
      proof (rule degree_sum_le)
        fix q assume "q  {0..p}"
        hence "q  p" by fastforce
        thus "degree (smult (inverse (real (p choose q)) * 
              coeff Q (p - q)) (Bernstein_Poly_01 q p))  p"
          by (auto simp add: degree_Bernstein degree_smult_le)
      qed simp
      hence "degree (reciprocal_poly p (j = 0..p.
            smult (inverse (real (p choose j)) * coeff Q (p - j)) 
            (Bernstein_Poly_01 j p)) p [:1, 1:])  p"
        by (auto simp add: degree_pcompose degree_reciprocal)
      hence "?L = 0" using hn by (auto simp add: coeff_eq_0)
      thus "?L = ?R" using hQ hn by (simp add: coeff_eq_0)
    qed
  qed
  hence "reciprocal_poly p P p [:1, 1:] = 
         reciprocal_poly p (j = 0..p. 
         smult (inverse (real (p choose j)) *
         coeff (reciprocal_poly p P p [:1, 1:]) (p-j))
         (Bernstein_Poly_01 j p)) p [:1, 1:]" 
    by (auto simp: degree_reciprocal degree_pcompose Q_def)
  hence "reciprocal_poly p P p ([:1, 1:] p [:-1, 1:]) =
         reciprocal_poly p (j = 0..p. smult (inverse (real (p choose j)) * 
         coeff (reciprocal_poly p P p [:1, 1:]) (p-j)) 
         (Bernstein_Poly_01 j p)) p ([:1, 1:] p [:-1, 1:])"
    by (auto simp: pcompose_assoc)
  hence "reciprocal_poly p P = reciprocal_poly p (j = 0..p. 
         smult (inverse (real (p choose j)) *
         coeff (reciprocal_poly p P p [:1, 1:]) (p-j)) (Bernstein_Poly_01 j p))" 
    by (auto simp: pcompose_pCons)
  hence "reciprocal_poly p (reciprocal_poly p P) = 
         reciprocal_poly p (reciprocal_poly p (j = 0..p. 
         smult (inverse (real (p choose j)) *
         coeff (reciprocal_poly p P p [:1, 1:]) (p-j)) (Bernstein_Poly_01 j p)))"
    by argo
  thus "P = (j = 0..p. smult (inverse (real (p choose j)) * 
        coeff (reciprocal_poly p P p [:1, 1:]) (p-j)) (Bernstein_Poly_01 j p))"
    using hP by (auto simp: reciprocal_reciprocal degree_sum_le degree_smult_le 
                 degree_Bernstein degree_add_le)
qed

lemma Bernstein_Poly_01_span1: 
  assumes hP: "degree P  p"
  shows "P  poly_vs.span {Bernstein_Poly_01 x p | x. x  p}"
proof -
  have "Bernstein_Poly_01 x p
          poly_vs.span {Bernstein_Poly_01 x p |x. x  p}"
    if "x  {0..p}" for x
  proof -
    have "n. Bernstein_Poly_01 x p = Bernstein_Poly_01 n p  n  p"
      using that by force
    then show 
      "Bernstein_Poly_01 x p  poly_vs.span {Bernstein_Poly_01 n p |n. n  p}"
      by (simp add: poly_vs.span_base)
  qed
  thus ?thesis
    apply (subst Bernstein_sum_01[OF hP])
    apply (rule poly_vs.span_sum)
    apply (rule poly_vs.span_scale)
    by blast
qed

lemma Bernstein_Poly_01_span:
  "poly_vs.span {Bernstein_Poly_01 x p | x. x  p} 
      = {x. degree x  p}"
  apply (subst monom_span[symmetric])
  apply (subst poly_vs.span_eq)
  by (auto simp: monom_span degree_Bernstein_le
      Bernstein_Poly_01_span1 degree_monom_eq)

subsection Bernstein coefficients and changes

definition Bernstein_coeffs_01 :: "nat  real poly  real list" where 
  "Bernstein_coeffs_01 p P = 
   [(inverse (real (p choose j)) * 
    coeff (reciprocal_poly p P p [:1, 1:]) (p-j)). j  [0..<(p+1)]]"

lemma length_Bernstein_coeffs_01: "length (Bernstein_coeffs_01 p P) = p + 1"
  by (auto simp: Bernstein_coeffs_01_def)

lemma nth_default_Bernstein_coeffs_01: assumes "degree P  p"
  shows "nth_default 0 (Bernstein_coeffs_01 p P) i = 
         inverse (p choose i) * coeff (reciprocal_poly p P p [:1, 1:]) (p-i)"
  apply (cases "p = i")
  using assms by (auto simp: Bernstein_coeffs_01_def nth_default_append
                  nth_default_Cons Nitpick.case_nat_unfold binomial_eq_0)

lemma Bernstein_coeffs_01_sum: assumes "degree P  p"
  shows "P = (j = 0..p. smult (nth_default 0 (Bernstein_coeffs_01 p P) j) 
             (Bernstein_Poly_01 j p))"
  apply (subst nth_default_Bernstein_coeffs_01[OF assms])
  apply (subst Bernstein_sum_01[OF assms])
  by argo

definition Bernstein_changes_01 :: "nat  real poly  int" where
  "Bernstein_changes_01 p P = nat (changes (Bernstein_coeffs_01 p P))"

lemma Bernstein_changes_01_def': 
  "Bernstein_changes_01 p P = nat (changes [(inverse (real (p choose j)) * 
     coeff (reciprocal_poly p P p [:1, 1:]) (p-j)). j  [0..<p + 1]])"
  by (simp add: Bernstein_changes_01_def Bernstein_coeffs_01_def)

lemma Bernstein_changes_01_eq_changes: 
  assumes hP: "degree P  p"
  shows "Bernstein_changes_01 p P = 
         changes (coeffs ((reciprocal_poly p P) p [:1, 1:]))"
proof (subst Bernstein_changes_01_def')
  have h: 
    "map (λj. inverse (real (p choose j)) * 
     coeff (reciprocal_poly p P p [:1, 1:]) (p - j)) [0..<p + 1] = 
     map (λj. inverse (real (p choose j)) * 
     nth_default 0 [nth_default 0 (coeffs (reciprocal_poly p P p [:1, 1:]))
                    (p - j). j  [0..<p + 1]] j) [0..<p + 1]"
  proof (rule map_cong)
    fix x
    assume "x  set [0..<p+1]"
    hence hx: "x  p" by fastforce
    moreover have 1:
      "length (map (λj. nth_default 0 
       (coeffs (reciprocal_poly p P p [:1, 1:])) (p - j)) [0..<p + 1])  Suc p"
      by force
    moreover have "length (coeffs (reciprocal_poly p P p [:1, 1:]))  Suc p"
    proof (cases "P=0")
      case False
      then have "reciprocal_poly p P p [:1, 1:]  0" 
        using hP by (simp add: Missing_Polynomial.pcompose_eq_0 reciprocal_0_iff)
      moreover have "Suc (degree (reciprocal_poly p P p [:1, 1:]))  Suc p"
        using hP by (auto simp: degree_pcompose degree_reciprocal)
      ultimately show ?thesis 
        using length_coeffs_degree by force
    qed (auto simp: reciprocal_0)
    ultimately have h: 
      "nth_default 0 (map (λj. nth_default 0 (coeffs 
       (reciprocal_poly p P p [:1, 1:])) (p - j)) [0..<p + 1]) x =
       nth_default 0 (coeffs (reciprocal_poly p P p [:1, 1:])) (p - x)"
      (is "?L = ?R")
    proof -
      have "?L = (map (λj. nth_default 0 (coeffs
            (reciprocal_poly p P p [:1, 1:])) (p - j)) [0..<p + 1]) ! x"
        using hx by (auto simp: nth_default_nth)
      also have "... =  nth_default 0 
          (coeffs (reciprocal_poly p P p [:1, 1:])) (p - [0..<p + 1] ! x)"
        apply (subst nth_map)
        using hx by auto
      also have "... = ?R"
        apply (subst nth_upt)
        using hx by auto
      finally show ?thesis .
    qed
    show "inverse (real (p choose x)) *
          coeff (reciprocal_poly p P p [:1, 1:]) (p - x) =
          inverse (real (p choose x)) *
          nth_default 0 (map (λj. nth_default 0 
          (coeffs (reciprocal_poly p P p [:1, 1:])) (p - j)) [0..<p + 1]) x"
      apply (subst h)
      apply (subst nth_default_coeffs_eq)
      by blast
  qed auto

  have 1: 
    "rev (map (λj. nth_default 0 (coeffs (reciprocal_poly p P p [:1, 1:])) 
     (p - j)) [0..<p + 1]) = map (λj. nth_default 0 (coeffs 
     (reciprocal_poly p P p [:1, 1:])) j) [0..<p + 1]"
  proof (subst rev_map, rule map_cong')
    have "q. (q  p  rev [q-p..<q+1] = map ((-) q) [0..<p+1])"
    proof (induction p)
      case 0
      then show ?case by simp
    next
      case (Suc p)
      have IH: "q. (q  p  rev [q-p..<q+1] = map ((-) q) [0..<p+1])"
        using Suc.IH by blast
      show ?case
      proof
        assume hq: "Suc p  q"
        then have h: "rev [q - p..<q + 1] = map ((-) (q)) [0..<p + 1]"
          apply (subst IH)
          using hq by auto
        have "[q - Suc p..<q + 1] = (q - Suc p) # [q - p..<q + 1]"
          by (simp add: Suc_diff_Suc Suc_le_lessD hq upt_conv_Cons)
        hence "rev [q - Suc p..<q + 1] = rev [q - p..<q + 1] @ [q - Suc p]"
          by force
        also have "... = map ((-) (q)) [0..<p + 1] @ [q - Suc p]"
          using h by blast
        also have "... = map ((-) q) [0..<Suc p + 1]"
          by force
        finally show "rev [q - Suc p..<q + 1] = map ((-) q) [0..<Suc p + 1]" .
      qed
    qed
    thus "rev [0..<p + 1] = map ((-) p) [0..<p + 1]"
      by force
  next
    fix y
    assume "y  set [0..<p + 1]"
    hence "y  p" by fastforce
    thus "nth_default 0 (coeffs (reciprocal_poly p P p [:1, 1:])) (p - (p - y)) =
          nth_default 0 (coeffs (reciprocal_poly p P p [:1, 1:])) y"
      by fastforce
  qed

  have 2: " f. f  0  degree f  p 
           map (nth_default 0 (coeffs f)) [0..<p + 1] = 
           coeffs f @ replicate (p - degree f) 0"
  proof (induction p)
    case 0
    then show ?case by (auto simp: degree_0_iff)
  next
    fix f
    case (Suc p)
    hence IH: "(f  0 
                degree f  p 
                map (nth_default 0 (coeffs f)) [0..<p + 1] =
                coeffs f @ replicate (p - degree f) 0)" by blast
    then show ?case
    proof (cases)
      assume h': "Suc p = degree f"
      hence h: "[0..<Suc p + 1] = [0..<length (coeffs f)]" 
        by (metis add_is_0 degree_0 length_coeffs plus_1_eq_Suc zero_neq_one)
      thus ?thesis
        apply (subst h)
        apply (subst map_nth_default)
        using h' by fastforce
    next
      assume h': "Suc p  degree f"
      show ?thesis
      proof
        assume hf: "f  0"
        show "degree f  Suc p 
            map (nth_default 0 (coeffs f)) [0..<Suc p + 1] =
            coeffs f @ replicate (Suc p - degree f) 0"
        proof
          assume "degree f  Suc p"
          hence 1: "degree f  p" using h' by fastforce
          hence 2: "map (nth_default 0 (coeffs f)) [0..<p + 1] =
                  coeffs f @ replicate (p - degree f) 0" using IH hf by blast
          have "map (nth_default 0 (coeffs f)) [0..<Suc p + 1] = 
                map (nth_default 0 (coeffs f)) [0..<p + 1] @
                     [nth_default 0 (coeffs f) (Suc p)]"
            by fastforce
          also have
            "... = coeffs f @ replicate (p - degree f) 0 @ [coeff f (Suc p)]"
            using 2 
            by (auto simp: nth_default_coeffs_eq)
          also have "... = coeffs f @ replicate (p - degree f) 0 @ [0]"
            using degree f  Suc p h' le_antisym le_degree by blast
          also have "... = coeffs f @ replicate (Suc p - degree f) 0" using 1
            by (simp add: Suc_diff_le replicate_app_Cons_same)
          finally show "map (nth_default 0 (coeffs f)) [0..<Suc p + 1] =
                coeffs f @ replicate (Suc p - degree f) 0" .
        qed
      qed
    qed
  qed
  
  thus "int (nat (changes (map (λj. inverse (real (p choose j)) *
        coeff (reciprocal_poly p P p [:1, 1:]) (p - j)) [0..<p + 1]))) =
        changes (coeffs (reciprocal_poly p P p [:1, 1:]))"
  proof cases
    assume hP: "P = 0"
    show "int (nat (changes (map (λj. inverse (real (p choose j)) *
          coeff (reciprocal_poly p P p [:1, 1:]) (p - j)) [0..<p + 1]))) =
          changes (coeffs (reciprocal_poly p P p [:1, 1:]))" (is "?L = ?R")
    proof -
      have "?L = int (nat (changes (map (λj. 0::real) [0..<p+1])))"
        using hP by (auto simp: reciprocal_0 changes_nonneg)
      also have "... = 0"
        apply (induction p)
        by (auto simp: map_replicate_trivial changes_nonneg
            replicate_app_Cons_same)
      also have "0 = changes ([]::real list)" by simp
      also have "... = ?R" using hP by (auto simp: reciprocal_0)
      finally show ?thesis .
    qed
  next
    assume hP': "P  0"
    thus ?thesis
      apply (subst h)
      apply (subst changes_scale)
        apply auto[2]
      apply (subst changes_rev[symmetric])
      apply (subst 1)
      apply (subst 2)
      apply (simp add: pcompose_eq_0 hP reciprocal_0_iff)
      using assms apply (auto simp: degree_reciprocal)[1]
      by (auto simp: changes_append_replicate_0 changes_nonneg)
  qed
qed

lemma Bernstein_changes_01_test: fixes P::"real poly"
  assumes hP: "degree P  p" and h0: "P  0"
  shows "proots_count P {x. 0 < x  x < 1}  Bernstein_changes_01 p P 
        even (Bernstein_changes_01 p P - proots_count P {x. 0 < x  x < 1})"
proof -
  let ?Q = "(reciprocal_poly p P) p [:1, 1:]"

  have 1: "changes (coeffs ?Q)  proots_count ?Q {x. 0 < x}  
        even (changes (coeffs ?Q) - proots_count ?Q {x. 0 < x})"
    apply (rule descartes_sign)
    by (simp add: Missing_Polynomial.pcompose_eq_0 h0 hP reciprocal_0_iff)
  
  have "((+) (1::real) ` Collect ((<) (0::real))) = {x. (1::real)<x}"
  proof
    show "{x::real. 1 < x}  (+) 1 ` Collect ((<) 0)"
    proof
      fix x::real assume "x  {x. 1 < x}"
      hence "1 < x" by simp
      hence "-1 + x  Collect ((<) 0)" by auto
      hence "1 + (-1 + x)  (+) 1 ` Collect ((<) 0)" by blast
      thus "x  (+) 1 ` Collect ((<) 0)" by argo
    qed
  qed auto
  hence 2:  "proots_count P {x. 0 < x  x < 1} = proots_count ?Q {x. 0 < x}"
    using assms
    by (auto simp: proots_pcompose reciprocal_0_iff proots_count_reciprocal')
  
  show ?thesis
    apply (subst Bernstein_changes_01_eq_changes[OF hP])
    apply (subst Bernstein_changes_01_eq_changes[OF hP])
    apply (subst 2)
    apply (subst 2)
    by (rule 1)
qed

subsection Expression as a Bernstein sum

lemma Bernstein_coeffs_01_0: "Bernstein_coeffs_01 p 0 = replicate (p+1) 0"
  by (auto simp: Bernstein_coeffs_01_def reciprocal_0 map_replicate_trivial
      replicate_append_same)

lemma Bernstein_coeffs_01_1: "Bernstein_coeffs_01 p 1 = replicate (p+1) 1"
proof -
  have "Bernstein_coeffs_01 p 1 =
     map (λj. inverse (real (p choose j)) *
     coeff (kp. smult (real (p choose k)) ([:0, 1:] ^ k)) (p - j)) [0..<(p+1)]"
    by (auto simp: Bernstein_coeffs_01_def reciprocal_1 monom_altdef
        hom_distribs pcompose_pCons poly_0_coeff_0[symmetric] poly_binomial)
  also have "... = map (λj. inverse (real (p choose j)) * 
             real (p choose (p - j))) [0..<(p+1)]"
    by (auto simp: monom_altdef[symmetric] coeff_sum binomial)
  also have "... = map (λj. 1) [0..<(p+1)]"
    apply (rule map_cong)
    subgoal by argo
    subgoal apply (subst binomial_symmetric)
      by auto
    done
  also have "... = replicate (p+1) 1"
    by (auto simp: map_replicate_trivial replicate_append_same)
  finally show ?thesis .
qed

lemma Bernstein_coeffs_01_x: assumes "p  0"
  shows "Bernstein_coeffs_01 p (monom 1 1) = [i/p. i  [0..<(p+1)]]"
proof -
  have 
    "Bernstein_coeffs_01 p (monom 1 1) = map (λj. inverse (real (p choose j)) *
     coeff (monom 1 (p - Suc 0) p [:1, 1:]) (p - j)) [0..<(p+1)]"
    using assms by (auto simp: Bernstein_coeffs_01_def reciprocal_monom)
  also have 
    "... = map (λj. inverse (real (p choose j)) *
     (kp - Suc 0. coeff (monom (real (p -  1 choose k)) k) (p - j))) [0..<(p+1)]"
    by (auto simp: monom_altdef hom_distribs pcompose_pCons poly_binomial coeff_sum)
  also have"... = map (λj. inverse (real (p choose j)) *
            real (p -  1 choose (p - j))) [0..<(p+1)]"
    by auto
  also have "... = map (λj. j/p) [0..<(p+1)]"
  proof (rule map_cong)
    fix x assume "x  set [0..<(p+1)]"
    hence "x  p" by force
    thus "inverse (real (p choose x)) * real (p - 1 choose (p - x)) =
          real x / real p"
    proof (cases "x = 0")
      show "x = 0  ?thesis"
        using assms by fastforce
      assume 1: "x  p" and 2: "x  0"
      hence "p - x  p - 1" by force
      hence "(p - 1 choose (p - x)) = (p - 1 choose (x - 1))"
        apply (subst binomial_symmetric)
        using 1 2 by auto
      hence "x * (p choose x) = p * (p - 1 choose (p - x))"
         using 2 times_binomial_minus1_eq by simp
       hence "real x * real (p choose x) = real p * real (p - 1 choose (p - x))"
         by (metis of_nat_mult)
       thus ?thesis using 1 2
         by (auto simp: divide_simps)
    qed
  qed blast
  finally show ?thesis .
qed

lemma Bernstein_coeffs_01_add: 
  assumes "degree P  p" and "degree Q  p"
  shows "nth_default 0 (Bernstein_coeffs_01 p (P + Q)) i = 
    nth_default 0 (Bernstein_coeffs_01 p P) i +
    nth_default 0 (Bernstein_coeffs_01 p Q) i"
  using assms by (auto simp: nth_default_Bernstein_coeffs_01 degree_add_le
                    reciprocal_add pcompose_add algebra_simps)

lemma Bernstein_coeffs_01_smult: 
  assumes "degree P  p"
  shows "nth_default 0 (Bernstein_coeffs_01 p (smult a P)) i =
          a * nth_default 0 (Bernstein_coeffs_01 p P) i"
  using assms
  by (auto simp: nth_default_Bernstein_coeffs_01 reciprocal_smult
      pcompose_smult)

end