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