Theory Bernstein

section Bernstein Polynomials over any finite interval

theory Bernstein
  imports "Bernstein_01"
begin

subsection Definition and relation to Bernstein Polynomials over [0, 1]

definition Bernstein_Poly :: "nat  nat  real  real  real poly" where
  "Bernstein_Poly j p c d = smult ((p choose j)/(d - c)^p)
      (((monom 1 j) p [:-c, 1:]) * (monom 1 (p-j) p [:d, -1:]))"

lemma Bernstein_Poly_altdef: 
  assumes "c  d" and "j  p"
  shows "Bernstein_Poly j p c d = smult (p choose j) 
            ([:-c/(d-c), 1/(d-c):]^j * [:d/(d-c), -1/(d-c):]^(p-j))" 
    (is "?L = ?R")
proof -
  have "?L = smult (p choose j) (smult ((1/(d - c))^j)
        (smult ((1/(d - c))^(p-j)) ([:-c, 1:]^j * [:d, -1:]^(p-j))))"
    using assms by (auto simp: Bernstein_Poly_def monom_altdef hom_distribs
                    pcompose_pCons smult_eq_iff field_simps power_add[symmetric])
  also have "... = ?R"
    apply (subst mult_smult_right[symmetric])
    apply (subst mult_smult_left[symmetric])
    apply (subst smult_power)
    apply (subst smult_power)
    by auto
  finally show ?thesis .
qed

lemma Bernstein_Poly_nonneg: 
  assumes "c  x" and "x  d"
  shows "poly (Bernstein_Poly j p c d) x  0"
  using assms by (auto simp: Bernstein_Poly_def poly_pcompose poly_monom)

lemma Bernstein_Poly_01: "Bernstein_Poly j p 0 1 = Bernstein_Poly_01 j p"
  by (auto simp: Bernstein_Poly_def Bernstein_Poly_01_def monom_altdef)

lemma Bernstein_Poly_rescale: 
  assumes "a  b"
  shows "Bernstein_Poly j p c d p [:a, 1:] p [:0, b-a:] 
            = Bernstein_Poly j p ((c-a)/(b-a)) ((d-a)/(b-a))"
  (is "?L = ?R")
proof -
  have "?R = smult (real (p choose j) 
      / ((d - a) / (b - a) - (c - a) / (b - a)) ^ p)
      ([:- ((c - a) / (b - a)), 1:] ^ j 
      * [:(d - a) / (b - a), - 1:] ^ (p - j))"
    by (auto simp: Bernstein_Poly_def monom_altdef hom_distribs 
        pcompose_pCons) 
  also have "... = smult (real (p choose j) / ((d - c) / (b - a)) ^ p)
              ([:- ((c - a) / (b - a)), 1:] ^ j * [:(d - a) / (b - a), - 1:] 
            ^ (p - j))"
    by argo
  also have "... = smult (real (p choose j) / (d - c) ^ p) 
      (smult ((b - a) ^ (p - j)) (smult ((b - a) ^ j)
      ([:- ((c - a) / (b - a)), 1:] ^ j * [:(d - a) / (b - a), - 1:] 
      ^ (p - j))))"
    by (auto simp: power_add[symmetric] power_divide)
  also have "... = smult (real (p choose j) / (d - c) ^ p)
              ([:- (c - a), b - a:] ^ j * [:d - a, -(b - a):] ^ (p - j))"
    apply (subst mult_smult_left[symmetric])
    apply (subst mult_smult_right[symmetric])
    using assms by (auto simp: smult_power)
  also have "... = ?L"
    using assms 
    by (auto simp: Bernstein_Poly_def monom_altdef pcompose_mult 
        pcompose_smult hom_distribs pcompose_pCons)
  finally show ?thesis by presburger
qed

lemma Bernstein_Poly_rescale_01: 
  assumes "c  d"
  shows "Bernstein_Poly j p c d p [:c, 1:] p [:0, d-c:] 
          = Bernstein_Poly_01 j p"
  apply (subst Bernstein_Poly_rescale)
  using assms by (auto simp: Bernstein_Poly_01)

lemma Bernstein_Poly_eq_rescale_01: 
  assumes "c  d"
  shows "Bernstein_Poly j p c d = Bernstein_Poly_01 j p 
            p [:0, 1/(d-c):] p [:-c, 1:]"
  apply (subst Bernstein_Poly_rescale_01[symmetric])
  using assms by (auto simp: pcompose_pCons pcompose_assoc[symmetric])

lemma coeff_Bernstein_sum: 
  fixes b::"nat  real" and p::nat and c d::real
  defines "P  (j = 0..p. (smult (b j) (Bernstein_Poly j p c d)))"
  assumes "i  p" and "c  d"
  shows "coeff ((reciprocal_poly p (P p [:c, 1:] 
      p [:0, d-c:])) p [:1, 1:]) (p - i) = (p choose i) * (b i)"
proof -
  have h: "P p [:c, 1:] p [:0, d-c:] 
      = (j = 0..p. (smult (b j) (Bernstein_Poly_01 j p)))"
    using assms 
    by (auto simp: P_def pcompose_sum pcompose_smult 
          pcompose_add Bernstein_Poly_rescale_01)
  then show ?thesis
    using coeff_Bernstein_sum_01 assms by simp
qed

lemma Bernstein_sum: 
  assumes "c  d" and "degree P  p"
  shows "P = (j = 0..p. smult (inverse (real (p choose j))
     * coeff (reciprocal_poly p (P p [:c, 1:] p [:0, d-c:]) 
      p [:1, 1:]) (p-j)) (Bernstein_Poly j p c d))"
  apply (subst Bernstein_Poly_eq_rescale_01)
  subgoal using assms by blast
  subgoal 
    apply (subst pcompose_smult[symmetric])
    apply (subst pcompose_sum[symmetric])
    apply (subst pcompose_smult[symmetric])
    apply (subst pcompose_sum[symmetric])
    apply (subst Bernstein_sum_01[symmetric])
    using assms by (auto simp: degree_pcompose pcompose_assoc[symmetric] 
        pcompose_pCons)
  done

lemma Bernstein_Poly_span1: 
  assumes "c  d" and "degree P  p"
  shows "P  poly_vs.span {Bernstein_Poly x p c d | x. x  p}"
proof (subst Bernstein_sum[OF assms], rule poly_vs.span_sum)
  fix x :: nat
  assume "x  {0..p}"
  then have "n. Bernstein_Poly x p c d = Bernstein_Poly n p c d  n  p"
    by auto
  then have 
    "Bernstein_Poly x p c d  poly_vs.span {Bernstein_Poly n p c d |n. n  p}"
    by (simp add: poly_vs.span_base)
  thus "smult (inverse (real (p choose x)) *
        coeff (reciprocal_poly p (P p [:c, 1:] p [:0, d - c:]) p [:1, 1:])
        (p - x)) (Bernstein_Poly x p c d)
          poly_vs.span {Bernstein_Poly x p c d |x. x  p}"
    by (rule poly_vs.span_scale)
qed

lemma Bernstein_Poly_span: 
  assumes "c  d" 
  shows "poly_vs.span {Bernstein_Poly x p c d | x. x  p} = {x. degree x  p}"
proof (subst Bernstein_Poly_01_span[symmetric], subst poly_vs.span_eq, rule conjI)
  show "{Bernstein_Poly x p c d |x. x  p}
       poly_vs.span {Bernstein_Poly_01 x p |x. x  p}"
    apply (subst Setcompr_subset)
    apply (rule allI, rule impI)
    apply (rule Bernstein_Poly_01_span1)
    using assms by (auto simp: degree_Bernstein_le Bernstein_Poly_eq_rescale_01
                    degree_pcompose)

  show "{Bernstein_Poly_01 x p |x. x  p}
       poly_vs.span {Bernstein_Poly x p c d |x. x  p}"
    apply (subst Setcompr_subset)
    apply (rule allI, rule impI)
    apply (rule Bernstein_Poly_span1)
    using assms by (auto simp: degree_Bernstein_le)
qed

lemma Bernstein_Poly_independent: assumes "c  d" 
  shows "poly_vs.independent {Bernstein_Poly x p c d | x. x  {..p}}"
proof (rule poly_vs.card_le_dim_spanning)
  show "{Bernstein_Poly x p c d |x. x  {.. p}}  {x. degree x  p}"
    using assms 
    by (auto simp: degree_Bernstein Bernstein_Poly_eq_rescale_01 degree_pcompose)
  show "{x. degree x  p}  poly_vs.span {Bernstein_Poly x p c d |x. x  {..p}}"
    using assms by (auto simp: Bernstein_Poly_span1)
  show "finite {Bernstein_Poly x p c d |x. x  {..p}}" by fastforce
  show "card {Bernstein_Poly x p c d |x. x  {..p}}  poly_vs.dim {x. degree x  p}"
    apply (rule le_trans)
     apply (subst image_Collect[symmetric], rule card_image_le, force)
    by (force simp: dim_degree)
qed

subsection Bernstein coefficients and changes over any interval

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

lemma Bernstein_coeffs_eq_rescale: assumes "c  d"
  shows "Bernstein_coeffs p c d P = Bernstein_coeffs_01 p (P p [:c, 1:] p [:0, d-c:])"
  using assms by (auto simp: pcompose_pCons pcompose_assoc[symmetric]
                  Bernstein_coeffs_def Bernstein_coeffs_01_def)

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

lemma Bernstein_coeffs_sum: assumes "c  d" and hP: "degree P  p"
  shows "P = (j = 0..p. smult (nth_default 0 (Bernstein_coeffs p c d P) j)
         (Bernstein_Poly j p c d))"
  apply (subst nth_default_Bernstein_coeffs[OF hP])
  apply (subst Bernstein_sum[OF assms])
  by argo

definition Bernstein_changes :: "nat  real  real  real poly  int" where
  "Bernstein_changes p c d P = nat (changes (Bernstein_coeffs p c d P))"

lemma Bernstein_changes_eq_rescale: assumes "c  d" and "degree P  p"
  shows "Bernstein_changes p c d P =
         Bernstein_changes_01 p (P p [:c, 1:] p [:0, d-c:])"
  using assms by (auto simp: Bernstein_coeffs_eq_rescale Bernstein_changes_def
                  Bernstein_changes_01_def)

text This is related and mostly equivalent to previous Descartes test \cite{li2019counting}
lemma Bernstein_changes_test: 
  fixes P::"real poly"
  assumes "degree P  p" and "P  0" and "c < d"
  shows "proots_count P {x. c < x  x < d}  Bernstein_changes p c d P 
        even (Bernstein_changes p c d P - proots_count P {x. c < x  x < d})"
proof -
  define Q where "Q=P p [:c, 1:] p [:0, d - c:]"

  have "int (proots_count Q {x. 0 < x  x < 1}) 
           Bernstein_changes_01 p Q 
              even (Bernstein_changes_01 p Q - 
                  int (proots_count Q {x. 0 < x  x < 1}))"
    unfolding Q_def
    apply (rule Bernstein_changes_01_test)
    subgoal using assms by fastforce
    subgoal using assms by (auto simp: pcompose_eq_0)
    done
  moreover have "proots_count P {x. c < x  x < d} =
            proots_count Q {x. 0 < x  x < 1}"
    unfolding Q_def
  proof (subst proots_pcompose)
    have "poly [:c, 1:] ` poly [:0, d - c:] ` {x. 0 < x  x < 1} =
        {x. c < x  x < d}" (is "?L = ?R")
    proof
      have "c + x * (d - c) < d" if "x < 1" for x
      proof - 
        have "x * (d - c) < 1 * (d - c)"
          using c < d that by force
        then show ?thesis by fastforce
      qed
      then show "?L  ?R"
        using assms by auto
    next
      show "?R  ?L"
      proof
        fix x::real assume "x  ?R"
        hence "c < x" and "x < d" by auto
        thus "x  ?L"
        proof (subst image_eqI)
          show "x = poly [:c, 1:] (x - c)" by force
          assume "c < x" and "x < d"
          thus "x - c  poly [:0, d - c:] ` {x. 0 < x  x < 1}"
          proof (subst image_eqI)
            show "x - c = poly [:0, d - c:] ((x - c)/(d - c))"
              using assms by fastforce
            assume "c < x" and "x < d"
            thus "(x - c) / (d - c)  {x. 0 < x  x < 1}"
              by auto
          qed fast
        qed fast
      qed
    qed
    then show "proots_count P {x. c < x  x < d} =
        proots_count (P p [:c, 1:]) 
        (poly [:0, d - c:] ` {x. 0 < x  x < 1})"
      using assms by (auto simp:proots_pcompose)
    show "P p [:c, 1:]  0"
      by (simp add: pcompose_eq_0 assms(2))
    show "degree [:0, d - c:] = 1"
      using assms by auto
  qed
  moreover have " Bernstein_changes p c d P = Bernstein_changes_01 p Q"
    unfolding Q_def
    apply (rule Bernstein_changes_eq_rescale)
    using assms by auto
  ultimately show ?thesis by auto
qed

subsection The control polygon of a polynomial

definition control_points ::
  "nat  real  real  real poly  (real × real) list"
where
  "control_points p c d P = 
   [(((real i)*d + (real (p - i))*c)/p, 
      nth_default 0 (Bernstein_coeffs p c d P) i).
      i  [0..<(p+1)]]"

lemma line_above: 
  fixes a b c d :: real and p :: nat and P :: "real poly"
  assumes hline: "i. i  p  a * (((real i)*d + (real (p - i))*c)/p) + b 
                  nth_default 0 (Bernstein_coeffs p c d P) i"
  and hp: "p  0" and hcd: "c  d" and hP: "degree P  p"
  shows "x. c  x  x  d  a*x + b  poly P x"
proof -
  fix x
  assume hc: "c  x" and  hd: "x  d"

  have bern_eq:"Bernstein_coeffs p c d [:b, a:] =
           [a*(real i * d + real (p - i) * c)/p + b. i  [0..<(p+1)]]"
  proof -
    have "Bernstein_coeffs p c d [:b, a:] = map (nth_default 0
          (Bernstein_coeffs_01 p ([:b, a:] p [:c, 1:] p [:0, d - c:])))
         [0..<p+1]"
      apply (subst Bernstein_coeffs_eq_rescale["OF" hcd])
      apply (subst map_nth_default[symmetric])
      apply (subst length_Bernstein_coeffs_01)
      by blast
    also have 
      "... = map (λi. a * (real i * d + real (p - i) * c) / real p + b) [0..<p + 1]"
    proof (rule map_cong)
      fix x assume hx: "x  set [0..<p + 1]"
      have "nth_default 0 (Bernstein_coeffs_01 p
            ([:b, a:] p [:c, 1:] p [:0, d - c:])) x =
            nth_default 0 (Bernstein_coeffs_01 p
            (smult (b + a*c) 1 + smult (a*(d - c)) (monom 1 1))) x"
      proof-
        have "[:b, a:] p [:c, 1:] p [:0, d - c:] =
                  smult (b + a*c) 1 + smult (a*(d - c)) (monom 1 1)"
          by (simp add: monom_altdef pcompose_pCons)
        then show ?thesis by auto
      qed
      also have "... = 
          nth_default 0 (Bernstein_coeffs_01 p (smult (b + a * c) 1)) x +
          nth_default 0 (Bernstein_coeffs_01 p (smult (a * (d - c)) (monom 1 1))) x"
        apply (subst Bernstein_coeffs_01_add)
        using hp by (auto simp: degree_monom_eq)
      also have "...  =
            (b + a*c) * nth_default 0 (Bernstein_coeffs_01 p 1) x +
            (a*(d - c)) * nth_default 0 (Bernstein_coeffs_01 p (monom 1 1)) x"
        apply (subst Bernstein_coeffs_01_smult)
        using hp by (auto simp: Bernstein_coeffs_01_smult degree_monom_eq)
      also have "... =
          (b + a * c) * (if x < p + 1 then 1 else 0) +
           a * (d - c) * (real (nth_default 0 [0..<p + 1] x) / real p)" 
        apply (subst Bernstein_coeffs_01_1, subst Bernstein_coeffs_01_x[OF hp])
        apply (subst nth_default_replicate_eq, subst nth_default_map_eq[of _ 0])
        by auto
      also have "... =
              (b + a * c) * (if x < p + 1 then 1 else 0) +
              a * (d - c) * (real ([0..<p + 1] ! x) / real p)"
        apply (subst nth_default_nth)
        using hx by auto
      also have "... = (b + a * c) * (if x < p + 1 then 1 else 0) +
              a * (d - c) * (real (0 + x) / real p)"
        apply (subst nth_upt)
        using hx by auto
      also have "... = a * (real x * d + real (p - x) * c) / real p + b"
        apply (subst of_nat_diff)
        using hx hp by (auto simp: field_simps)
      finally show "nth_default 0 (Bernstein_coeffs_01 p
                    ([:b, a:] p [:c, 1:] p [:0, d - c:])) x =
                    a * (real x * d + real (p - x) * c) / real p + b" .
    qed blast
    finally show ?thesis .
  qed

  have nth_default_geq:"nth_default 0 (Bernstein_coeffs p c d [:b, a:]) i 
           nth_default 0 (Bernstein_coeffs p c d P) i" for i
  proof -
    show "nth_default 0 (Bernstein_coeffs p c d [:b, a:]) i 
          nth_default 0 (Bernstein_coeffs p c d P) i"
    proof cases
      define p1 where "p1  p+1"
      assume h: "i  p"
      hence "nth_default 0 (Bernstein_coeffs p c d P) i 
             a * (((real i)*d + (real (p - i))*c)/p) + b"
        by (rule hline)
      also have "... = nth_default 0 (map (λi. a * (real i * d 
          + real (p - i) * c) / real p + b) [0..<p + 1]) i"
        apply (subst p1_def[symmetric])
        using h apply (auto simp: nth_default_def)
        by (auto simp: p1_def)
      also have "... = nth_default 0 (Bernstein_coeffs p c d [:b, a:]) i"
        using bern_eq by simp
      finally show ?thesis .
    next
      assume h: "¬i  p"
      thus ?thesis
        using assms 
        by (auto simp: nth_default_def Bernstein_coeffs_eq_rescale
                        length_Bernstein_coeffs_01)
    qed
  qed
  
  have "poly P x = (k = 0..p.
        poly (smult (nth_default 0 (Bernstein_coeffs p c d P) k)
        (Bernstein_Poly k p c d)) x)"
    apply (subst Bernstein_coeffs_sum[OF hcd hP])
    by (rule poly_sum)
  also have "...  (k = 0..p.
      poly (smult (nth_default 0 (Bernstein_coeffs p c d [:b, a:]) k)
        (Bernstein_Poly k p c d)) x)"
    apply (rule sum_mono)
    using mult_right_mono[OF nth_default_geq] Bernstein_Poly_nonneg[OF hc hd]
    by auto
  also have "... = poly [:b, a:] x"
    apply (subst(2) Bernstein_coeffs_sum[of c d "[:b, a:]" p])
    using assms apply auto[2]
    by (rule poly_sum[symmetric])
  also have "... = a*x + b" by force
  finally show "poly P x  a*x + b" .
qed

end