Theory Bernstein

section Bernstein Polynomials over any finite interval

theory Bernstein
imports
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
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:
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))))"
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
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}"
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 ::
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"
show "degree [:0, d - c:] = 1"
using assms by auto
qed
moreover have
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)"
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"
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 _ ])
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
by (rule poly_sum[symmetric])
also have "... = a*x + b" by force
finally show "poly P x  a*x + b" .
qed

end