Theory Rational_FPS_Asymptotics
theory Rational_FPS_Asymptotics
imports
"HOL-Library.Landau_Symbols"
"Polynomial_Factorization.Square_Free_Factorization"
"HOL-Real_Asymp.Real_Asymp"
"Count_Complex_Roots.Count_Complex_Roots"
Linear_Homogenous_Recurrences
Linear_Inhomogenous_Recurrences
RatFPS
Rational_FPS_Solver
"HOL-Library.Code_Target_Numeral"
begin
lemma poly_asymp_equiv:
assumes "p ≠ 0" and "F ≤ at_infinity"
shows "poly p ∼[F] (λx. lead_coeff p * x ^ degree p)"
proof -
have poly_pCons': "poly (pCons a q) = (λx. a + x * poly q x)" for a :: 'a and q
by (simp add: fun_eq_iff)
show ?thesis using assms(1)
proof (induction p)
case (pCons a p)
define n where "n = Suc (degree p)"
show ?case
proof (cases "p = 0")
case [simp]: False
hence *: "poly p ∼[F] (λx. lead_coeff p * x ^ degree p)"
by (intro pCons.IH)
have "poly (pCons a p) = (λx. a + x * poly p x)"
by (simp add: poly_pCons')
moreover have "… ∼[F] (λx. lead_coeff p * x ^ n)"
proof (subst asymp_equiv_add_left)
have "(λx. x * poly p x) ∼[F] (λx. x * (lead_coeff p * x ^ degree p))"
by (intro asymp_equiv_intros *)
also have "… = (λx. lead_coeff p * x ^ n)" by (simp add: n_def mult_ac)
finally show "(λx. x * poly p x) ∼[F] …" .
next
have "filterlim (λx. x) at_infinity F"
by (simp add: filterlim_def assms)
hence "(λx. x ^ n) ∈ ω[F](λ_. 1 :: 'a)" unfolding smallomega_1_conv_filterlim
by (intro Limits.filterlim_power_at_infinity filterlim_ident) (auto simp: n_def)
hence "(λx. a) ∈ o[F](λx. x ^ n)" unfolding smallomega_iff_smallo[symmetric]
by (cases "a = 0") auto
thus "(λx. a) ∈ o[F](λx. lead_coeff p * x ^ n)"
by simp
qed
ultimately show ?thesis by (simp add: n_def)
qed auto
qed auto
qed
lemma poly_bigtheta:
assumes "p ≠ 0" and "F ≤ at_infinity"
shows "poly p ∈ Θ[F](λx. x ^ degree p)"
proof -
have "poly p ∼[F] (λx. lead_coeff p * x ^ degree p)"
by (intro poly_asymp_equiv assms)
thus ?thesis using assms by (auto dest!: asymp_equiv_imp_bigtheta)
qed
lemma poly_bigo:
assumes "F ≤ at_infinity" and "degree p ≤ k"
shows "poly p ∈ O[F](λx. x ^ k)"
proof (cases "p = 0")
case True
hence "poly p = (λ_. 0)" by (auto simp: fun_eq_iff)
thus ?thesis by simp
next
case False
have *: "(λx. x ^ (k - degree p)) ∈ Ω[F](λx. 1)"
proof (cases "k = degree p")
case False
hence "(λx. x ^ (k - degree p)) ∈ ω[F](λ_. 1)"
unfolding smallomega_1_conv_filterlim using assms False
by (intro Limits.filterlim_power_at_infinity filterlim_ident)
(auto simp: filterlim_def)
thus ?thesis by (rule landau_omega.small_imp_big)
qed auto
have "poly p ∈ Θ[F](λx. x ^ degree p * 1)"
using poly_bigtheta[OF False assms(1)] by simp
also have "(λx. x ^ degree p * 1) ∈ O[F](λx. x ^ degree p * x ^ (k - degree p))" using *
by (intro landau_o.big.mult landau_o.big_refl) (auto simp: bigomega_iff_bigo)
also have "(λx::'a. x ^ degree p * x ^ (k - degree p)) = (λx. x ^ k)"
using assms by (simp add: power_add [symmetric])
finally show ?thesis .
qed
lemma reflect_poly_dvdI:
fixes p q :: "'a::{comm_semiring_1,semiring_no_zero_divisors} poly"
assumes "p dvd q"
shows "reflect_poly p dvd reflect_poly q"
using assms by (auto simp: reflect_poly_mult)
lemma smult_altdef: "smult c p = [:c:] * p"
by (induction p) (auto simp: mult_ac)
lemma smult_power: "smult (c ^ n) (p ^ n) = (smult c p) ^ n"
proof -
have "smult (c ^ n) (p ^ n) = [:c ^ n:] * p ^ n"
by simp
also have "[:c:] ^ n = [:c ^ n:]"
by (induction n) (auto simp: mult_ac)
hence "[:c ^ n:] = [:c:] ^ n" ..
also have "… * p ^ n = ([:c:] * p) ^ n"
by (rule power_mult_distrib [symmetric])
also have "… = (smult c p) ^ n" by simp
finally show ?thesis .
qed
lemma order_reflect_poly_ge:
fixes c :: "'a :: field"
assumes "c ≠ 0" and "p ≠ 0"
shows "order c (reflect_poly p) ≥ order (1 / c) p"
proof -
have "reflect_poly ([:-(1 / c), 1:] ^ order (1 / c) p) dvd reflect_poly p"
by (intro reflect_poly_dvdI, subst order_divides) auto
also have "reflect_poly ([:-(1 / c), 1:] ^ order (1 / c) p) =
smult ((-1 / c) ^ order (1 / c) p) ([:-c, 1:] ^ order (1 / c) p)"
using assms by (simp add: reflect_poly_power reflect_poly_pCons smult_power)
finally have "([:-c, 1:] ^ order (1 / c) p) dvd reflect_poly p"
by (rule smult_dvd_cancel)
with ‹p ≠ 0› show ?thesis by (subst (asm) order_divides) auto
qed
lemma order_reflect_poly:
fixes c :: "'a :: field"
assumes "c ≠ 0" and "coeff p 0 ≠ 0"
shows "order c (reflect_poly p) = order (1 / c) p"
proof (rule antisym)
from assms show "order c (reflect_poly p) ≥ order (1 / c) p"
by (intro order_reflect_poly_ge) auto
next
from assms have "order (1 / (1 / c)) (reflect_poly p) ≤
order (1 / c) (reflect_poly (reflect_poly p))"
by (intro order_reflect_poly_ge) auto
with assms show "order c (reflect_poly p) ≤ order (1 / c) p"
by simp
qed
lemma poly_reflect_eq_0_iff:
"poly (reflect_poly p) (x :: 'a :: field) = 0 ⟷ p = 0 ∨ x ≠ 0 ∧ poly p (1 / x) = 0"
by (cases "x = 0") (auto simp: poly_reflect_poly_nz inverse_eq_divide)
theorem ratfps_nth_bigo:
fixes q :: "complex poly"
assumes "R > 0"
assumes roots1: "⋀z. z ∈ ball 0 (1 / R) ⟹ poly q z ≠ 0"
assumes roots2: "⋀z. z ∈ sphere 0 (1 / R) ⟹ poly q z = 0 ⟹ order z q ≤ Suc k"
shows "fps_nth (fps_of_poly p / fps_of_poly q) ∈ O(λn. of_nat n ^ k * of_real R ^ n)"
proof -
define q' where "q' = reflect_poly q"
from roots1[of 0] and ‹R > 0› have [simp]: "coeff q 0 ≠ 0" "q ≠ 0"
by (auto simp: poly_0_coeff_0)
from ratfps_closed_form_exists[OF this(1), of p]
obtain r rs where closed_form:
"⋀n. (fps_of_poly p / fps_of_poly q) $ n =
coeff r n + (∑c | poly (reflect_poly q) c = 0. poly (rs c) (of_nat n) * c ^ n)"
"⋀z. poly (reflect_poly q) z = 0 ⟹ degree (rs z) ≤ order z (reflect_poly q) - 1"
by blast
have "fps_nth (fps_of_poly p / fps_of_poly q) =
(λn. coeff r n + (∑c | poly q' c = 0. poly (rs c) (of_nat n) * c ^ n))"
by (intro ext, subst closed_form) (simp_all add: q'_def)
also have "… ∈ O(λn. of_nat n ^ k * of_real R ^ n)"
proof (intro sum_in_bigo big_sum_in_bigo)
have "eventually (λn. coeff r n = 0) at_top"
using MOST_nat coeff_eq_0 cofinite_eq_sequentially by force
hence "coeff r ∈ Θ(λ_. 0)" by (rule bigthetaI_cong)
also have "(λ_. 0 :: complex) ∈ O(λn. of_nat n ^ k * of_real R ^ n)"
by simp
finally show "coeff r ∈ O(λn. of_nat n ^ k * of_real R ^ n)" .
next
fix c assume c: "c ∈ {c. poly q' c = 0}"
hence [simp]: "c ≠ 0" by (auto simp: q'_def)
show "(λn. poly (rs c) n * c ^ n) ∈ O(λn. of_nat n ^ k * of_real R ^ n)"
proof (cases "norm c = R")
case True
show ?thesis
proof (intro landau_o.big.mult landau_o.big.compose[OF poly_bigo tendsto_of_nat])
have "degree (rs c) ≤ order c (reflect_poly q) - 1"
using c by (intro closed_form(2)) (auto simp: q'_def)
also have "order c (reflect_poly q) = order (1 / c) q"
using c by (intro order_reflect_poly) (auto simp: q'_def)
also {
have "order (1 / c) q ≤ Suc k" using ‹R > 0› and True and c
by (intro roots2) (auto simp: q'_def norm_divide poly_reflect_eq_0_iff)
moreover have "order (1 / c) q ≠ 0"
using order_root[of q "1 / c"] c by (auto simp: q'_def poly_reflect_eq_0_iff)
ultimately have "order (1 / c) q - 1 ≤ k" by simp
}
finally show "degree (rs c) ≤ k" .
next
have "(λn. norm (c ^ n)) ∈ O(λn. norm (complex_of_real R ^ n))"
using True and ‹R > 0› by (simp add: norm_power)
thus "(λn. c ^ n) ∈ O(λn. complex_of_real R ^ n)"
by (subst (asm) landau_o.big.norm_iff)
qed auto
next
case False
hence "norm c < R" using c and roots1[of "1/c"] and ‹R > 0›
by (cases "norm c" R rule: linorder_cases)
(auto simp: q'_def poly_reflect_eq_0_iff norm_divide field_simps)
define l where "l = degree (rs c)"
have "(λn. poly (rs c) (of_nat n) * c ^ n) ∈ O(λn. of_nat n ^ l * c ^ n)"
by (intro landau_o.big.mult landau_o.big.compose[OF poly_bigo tendsto_of_nat])
(auto simp: l_def)
also have "(λn. of_nat n ^ l * c ^ n) ∈ O(λn. of_nat n ^ k * of_real R ^ n)"
proof (subst landau_o.big.norm_iff [symmetric])
have "(λn. real n ^ l) ∈ O(λn. real n ^ k * (R / norm c) ^ n)"
using ‹norm c < R› and ‹R > 0› by real_asymp
hence "(λn. real n ^ l * norm c ^ n) ∈ O(λn. real n ^ k * R ^ n)"
by (simp add: power_divide landau_o.big.divide_eq1)
thus "(λx. norm (of_nat x ^ l * c ^ x)) ∈
O(λx. norm (of_nat x ^ k * complex_of_real R ^ x))"
unfolding norm_power norm_mult using ‹R > 0› by simp
qed
finally show ?thesis .
qed
qed
finally show ?thesis .
qed
lemma order_power: "p ≠ 0 ⟹ order c (p ^ n) = n * order c p"
by (induction n) (auto simp: order_mult)
lemma same_root_imp_not_coprime:
assumes "poly p x = 0" and "poly q (x :: 'a :: {factorial_ring_gcd,semiring_gcd_mult_normalize}) = 0"
shows "¬coprime p q"
proof
assume "coprime p q"
from assms have "[:-x, 1:] dvd p" and "[:-x, 1:] dvd q"
by (simp_all add: poly_eq_0_iff_dvd)
hence "[:-x, 1:] dvd gcd p q" by (simp add: poly_eq_0_iff_dvd)
also from ‹coprime p q› have "gcd p q = 1"
by (rule coprime_imp_gcd_eq_1)
finally show False by (elim is_unit_polyE) auto
qed
lemma ratfps_nth_bigo_square_free_factorization:
fixes p :: "complex poly"
assumes "square_free_factorization q (b, cs)"
assumes "q ≠ 0" and "R > 0"
assumes roots1: "⋀c l. (c, l) ∈ set cs ⟹ ∀x∈ball 0 (1 / R). poly c x ≠ 0"
assumes roots2: "⋀c l. (c, l) ∈ set cs ⟹ l > Suc k ⟹ ∀x∈sphere 0 (1 / R). poly c x ≠ 0"
shows "fps_nth (fps_of_poly p / fps_of_poly q) ∈ O(λn. of_nat n ^ k * of_real R ^ n)"
proof -
from assms(1) have q: "q = smult b (∏(c, l)∈set cs. c ^ l)"
unfolding square_free_factorization_def prod.case by blast
with ‹q ≠ 0› have [simp]: "b ≠ 0" by auto
note sff = square_free_factorizationD[OF assms(1)]
from sff(2)[of 0] have [simp]: "(0, x) ∉ set cs" for x by auto
from assms(1) have coprime: "c1 = c2" "m = n"
if "¬coprime c1 c2" "(c1, m) ∈ set cs" "(c2, n) ∈ set cs" for c1 c2 m n
using that by (auto simp: square_free_factorization_def case_prod_unfold)
show ?thesis
proof (rule ratfps_nth_bigo)
fix z :: complex assume z: "z ∈ ball 0 (1 / R)"
show "poly q z ≠ 0"
proof
assume "poly q z = 0"
then obtain c l where cl: "(c, l) ∈ set cs" and "poly c z = 0"
by (auto simp: q poly_prod image_iff)
with roots1[of c l] and z show False by auto
qed
next
fix z :: complex assume z: "z ∈ sphere 0 (1 / R)"
have order: "order z q = order z (∏(c, l)∈set cs. c ^ l)"
by (simp add: order_smult q)
also have "… = (∑x∈set cs. order z (case x of (c, l) ⇒ c ^ l))"
by (subst order_prod) (auto dest: coprime)
also have "… = (∑(c, l)∈set cs. l * order z c)"
unfolding case_prod_unfold by (intro sum.cong refl, subst order_power) auto
finally have "order z q = …" .
show "order z q ≤ Suc k"
proof (cases "∃c0 l0. (c0, l0) ∈ set cs ∧ poly c0 z = 0")
case False
have "order z q = (∑(c, l)∈set cs. l * order z c)" by fact
also have "order z c = 0" if "(c, l) ∈ set cs" for c l
using False that by (auto simp: order_root)
hence "(∑(c, l)∈set cs. l * order z c) = 0"
by (intro sum.neutral) auto
finally show "order z q ≤ Suc k" by simp
next
case True
then obtain c0 l0 where cl0: "(c0, l0) ∈ set cs" "poly c0 z = 0"
by blast
have "order z q = (∑(c, l)∈set cs. l * order z c)" by fact
also have "… = l0 * order z c0 + (∑(c, l) ∈ set cs - {(c0, l0)}. l * order z c)"
using cl0 by (subst sum.remove[of _ "(c0, l0)"]) auto
also have "(∑(c, l) ∈ set cs - {(c0, l0)}. l * order z c) = 0"
proof (intro sum.neutral ballI, goal_cases)
case (1 cl)
then obtain c l where [simp]: "cl = (c, l)" and cl: "(c, l) ∈ set cs" "(c0, l0) ≠ (c, l)"
by (cases cl) auto
from cl and cl0 and coprime[of c c0 l l0] have "coprime c c0"
by auto
with same_root_imp_not_coprime[of c z c0] and cl0 have "poly c z ≠ 0" by auto
thus ?case by (auto simp: order_root)
qed
also have "square_free c0" using cl0 assms(1)
by (auto simp: square_free_factorization_def)
hence "rsquarefree c0" by (rule square_free_rsquarefree)
with cl0 have "order z c0 = 1"
by (auto simp: rsquarefree_def' order_root intro: antisym)
finally have "order z q = l0" by simp
also from roots2[OF cl0(1)] cl0(2) z have "l0 ≤ Suc k"
by (cases l0 "Suc k" rule: linorder_cases) auto
finally show "order z q ≤ Suc k" by simp
qed
qed fact+
qed
lemma proots_within_card_zero_iff:
assumes "p ≠ (0 :: 'a :: idom poly)"
shows "card (proots_within p A) = 0 ⟷ (∀x∈A. poly p x ≠ 0)"
using assms by (subst card_0_eq) (auto intro: finite_proots)
lemma ratfps_nth_bigo_square_free_factorization':
fixes p :: "complex poly"
assumes "square_free_factorization q (b, cs)"
assumes "q ≠ 0" and "R > 0"
assumes roots1: "list_all (λcl. proots_ball_card (fst cl) 0 (1 / R) = 0) cs"
assumes roots2: "list_all (λcl. proots_sphere_card (fst cl) 0 (1 / R) = 0)
(filter (λcl. snd cl > Suc k) cs)"
shows "fps_nth (fps_of_poly p / fps_of_poly q) ∈ O(λn. of_nat n ^ k * of_real R ^ n)"
proof (rule ratfps_nth_bigo_square_free_factorization[OF assms(1)])
note sff = square_free_factorizationD[OF assms(1)]
from sff(2)[of 0] have [simp]: "(0, x) ∉ set cs" for x by auto
from assms(1) have q: "q = smult b (∏(c, l)∈set cs. c ^ l)"
unfolding square_free_factorization_def prod.case by blast
with ‹q ≠ 0› have [simp]: "b ≠ 0" by auto
show "∀x∈ball 0 (1 / R). poly c x ≠ 0" if "(c, l) ∈ set cs" for c l
proof -
from roots1 that have "card (proots_within c (ball 0 (1 / R))) = 0"
by (auto simp: proots_ball_card_def list_all_def)
with that show ?thesis by (subst (asm) proots_within_card_zero_iff) auto
qed
show "∀x∈sphere 0 (1 / R). poly c x ≠ 0" if "(c, l) ∈ set cs" "l > Suc k" for c l
proof -
from roots2 that have "card (proots_within c (sphere 0 (1 / R))) = 0"
by (auto simp: proots_sphere_card_def list_all_def)
with that show ?thesis by (subst (asm) proots_within_card_zero_iff) auto
qed
qed fact+
definition ratfps_has_asymptotics where
"ratfps_has_asymptotics q k R ⟷ q ≠ 0 ∧ R > 0 ∧
(let cs = snd (yun_factorization gcd q)
in list_all (λcl. proots_ball_card (fst cl) 0 (1 / R) = 0) cs ∧
list_all (λcl. proots_sphere_card (fst cl) 0 (1 / R) = 0) (filter (λcl. snd cl > Suc k) cs))"
lemma ratfps_has_asymptotics_correct:
assumes "ratfps_has_asymptotics q k R"
shows "fps_nth (fps_of_poly p / fps_of_poly q) ∈ O(λn. of_nat n ^ k * of_real R ^ n)"
proof (rule ratfps_nth_bigo_square_free_factorization')
show "square_free_factorization q (fst (yun_factorization gcd q), snd (yun_factorization gcd q))"
by (rule yun_factorization) simp
qed (insert assms, auto simp: ratfps_has_asymptotics_def Let_def list_all_def)
value "map (fps_nth (fps_of_poly [:0, 1:] / fps_of_poly [:1, -1, -1 :: real:])) [0..<5]"
method ratfps_bigo = (rule ratfps_has_asymptotics_correct; eval)
lemma "fps_nth (fps_of_poly [:0, 1:] / fps_of_poly [:1, -1, -1 :: complex:]) ∈
O(λn. of_nat n ^ 0 * complex_of_real 1.618034 ^ n)"
by ratfps_bigo
lemma "fps_nth (fps_of_poly 1 / fps_of_poly [:1, -3, 3, -1 :: complex:]) ∈
O(λn. of_nat n ^ 2 * complex_of_real 1 ^ n)"
by ratfps_bigo
lemma "fps_nth (fps_of_poly f / fps_of_poly [:5, 4, 3, 2, 1 :: complex:]) ∈
O(λn. of_nat n ^ 0 * complex_of_real 0.69202 ^ n)"
by ratfps_bigo
end