# Theory RatFPS

(*
File:    RatFPS.thy
Author:  Manuel Eberl, TU München
*)
section ‹Rational formal power series›
theory RatFPS
imports
Complex_Main
"HOL-Computational_Algebra.Computational_Algebra"
"HOL-Computational_Algebra.Polynomial_Factorial"
begin

subsection ‹Some auxiliary›

abbreviation constant_term :: "'a poly ⇒ 'a::zero"
where "constant_term p ≡ coeff p 0"

lemma coeff_0_mult: "coeff (p * q) 0 = coeff p 0 * coeff q 0"
by (simp add: coeff_mult)

lemma coeff_0_div:
assumes "coeff p 0 ≠ 0"
assumes "(q :: 'a :: field poly) dvd p"
shows   "coeff (p div q) 0 = coeff p 0 div coeff q 0"
proof (cases "q = 0")
case False
from assms have "p = p div q * q" by simp
also have "coeff … 0 = coeff (p div q) 0 * coeff q 0" by (simp add: coeff_0_mult)
finally show ?thesis using assms by auto
qed simp_all

assumes "coeff (snd (quot_of_fract x)) 0 ≠ 0" "coeff (snd (quot_of_fract y)) 0 ≠ 0"
shows   "coeff (snd (quot_of_fract (x + y))) 0 ≠ 0"
proof -
define num where "num = fst (quot_of_fract x) * snd (quot_of_fract y) +
snd (quot_of_fract x) * fst (quot_of_fract y)"
define denom where "denom = snd (quot_of_fract x) * snd (quot_of_fract y)"
define z where "z = (num, denom)"
from assms have "snd z ≠ 0" by (auto simp: denom_def z_def)
from normalize_quotE'[OF this] guess d . note d = this
from assms have z: "coeff (snd z) 0 ≠ 0" by (simp add: z_def denom_def coeff_0_mult)

have "coeff (snd (quot_of_fract (x + y))) 0 = coeff (snd (normalize_quot z)) 0"
by (simp add: quot_of_fract_add Let_def case_prod_unfold z_def num_def denom_def)
also from z have "… ≠ 0" using d by (simp add: d coeff_0_mult)
finally show ?thesis .
qed

lemma coeff_0_normalize_quot_nonzero [simp]:
assumes "coeff (snd x) 0 ≠ 0"
shows   "coeff (snd (normalize_quot x)) 0 ≠ 0"
proof -
from assms have "snd x ≠ 0" by auto
from normalize_quotE'[OF this] guess d .
with assms show ?thesis by (auto simp: coeff_0_mult)
qed

abbreviation numerator :: "'a fract ⇒ 'a::{ring_gcd,idom_divide,semiring_gcd_mult_normalize}"
where "numerator x ≡ fst (quot_of_fract x)"

abbreviation denominator :: "'a fract ⇒ 'a::{ring_gcd,idom_divide,semiring_gcd_mult_normalize}"
where "denominator x ≡ snd (quot_of_fract x)"

declare unit_factor_snd_quot_of_fract [simp]
normalize_snd_quot_of_fract [simp]

lemma constant_term_denominator_nonzero_imp_constant_term_denominator_div_gcd_nonzero:
"constant_term (denominator x div gcd a (denominator x)) ≠ 0"
if "constant_term (denominator x) ≠ 0"
using that coeff_0_normalize_quot_nonzero [of "(a, denominator x)"]
normalize_quot_proj(2) [of "denominator x" a]
by simp

subsection ‹The type of rational formal power series›

typedef (overloaded) 'a :: field_gcd ratfps =
"{x :: 'a poly fract. constant_term (denominator x) ≠ 0}"
by (rule exI [of _ 0]) simp

setup_lifting type_definition_ratfps

instantiation ratfps :: (field_gcd) idom
begin

lift_definition zero_ratfps :: "'a ratfps" is "0" by simp

lift_definition one_ratfps :: "'a ratfps" is "1" by simp

lift_definition uminus_ratfps :: "'a ratfps ⇒ 'a ratfps" is "uminus"
by (simp add: quot_of_fract_uminus case_prod_unfold Let_def)

lift_definition plus_ratfps :: "'a ratfps ⇒ 'a ratfps ⇒ 'a ratfps" is "(+)"

lift_definition minus_ratfps :: "'a ratfps ⇒ 'a ratfps ⇒ 'a ratfps" is "(-)"
(simp_all add: quot_of_fract_uminus Let_def case_prod_unfold)

lift_definition times_ratfps :: "'a ratfps ⇒ 'a ratfps ⇒ 'a ratfps" is "(*)"
by (simp add: quot_of_fract_mult Let_def case_prod_unfold coeff_0_mult
constant_term_denominator_nonzero_imp_constant_term_denominator_div_gcd_nonzero)

instance
by (standard; transfer) (simp_all add: ring_distribs)

end

fun ratfps_nth_aux :: "('a::field) poly ⇒ nat ⇒ 'a"
where
"ratfps_nth_aux p 0 = inverse (coeff p 0)"
| "ratfps_nth_aux p n =
- inverse (coeff p 0) * sum (λi. coeff p i * ratfps_nth_aux p (n - i)) {1..n}"

lemma ratfps_nth_aux_correct: "ratfps_nth_aux p n = natfun_inverse (fps_of_poly p) n"
by (induction p n rule: ratfps_nth_aux.induct) simp_all

lift_definition ratfps_nth :: "'a :: field_gcd ratfps ⇒ nat ⇒ 'a" is
"λx n. let (a,b) = quot_of_fract x
in  (∑i = 0..n. coeff a i * ratfps_nth_aux b (n - i))" .

lift_definition ratfps_subdegree :: "'a :: field_gcd ratfps ⇒ nat" is
"λx. poly_subdegree (fst (quot_of_fract x))" .

context
includes lifting_syntax
begin

lemma RatFPS_parametric: "(rel_prod (=) (=) ===> (=))
(λ(p,q). if coeff q 0 = 0 then 0 else quot_to_fract (p, q))
(λ(p,q). if coeff q 0 = 0 then 0 else quot_to_fract (p, q))"
by transfer_prover

end

lemma normalize_quot_quot_of_fract [simp]:
"normalize_quot (quot_of_fract x) = quot_of_fract x"
by (rule normalize_quot_id, rule quot_of_fract_in_normalized_fracts)

context
assumes "SORT_CONSTRAINT('a::field_gcd)"
begin

lift_definition quot_of_ratfps :: "'a ratfps ⇒ ('a poly × 'a poly)" is
"quot_of_fract :: 'a poly fract ⇒ ('a poly × 'a poly)" .

lift_definition quot_to_ratfps :: "('a poly × 'a poly) ⇒ 'a ratfps" is
"λ(x,y). let (x',y') = normalize_quot (x,y)
in  if coeff y' 0 = 0 then 0 else quot_to_fract (x',y')"
by (simp add: case_prod_unfold Let_def quot_of_fract_quot_to_fract)

lemma quot_to_ratfps_quot_of_ratfps [code abstype]:
"quot_to_ratfps (quot_of_ratfps x) = x"
by transfer (simp add: case_prod_unfold Let_def)

lemma coeff_0_snd_quot_of_ratfps_nonzero [simp]:
"coeff (snd (quot_of_ratfps x)) 0 ≠ 0"
by transfer simp

lemma quot_of_ratfps_quot_to_ratfps:
"coeff (snd x) 0 ≠ 0 ⟹ x ∈ normalized_fracts ⟹ quot_of_ratfps (quot_to_ratfps x) = x"
by transfer (simp add: Let_def case_prod_unfold coeff_0_normalize_quot_nonzero
quot_of_fract_quot_to_fract normalize_quot_id)

lemma quot_of_ratfps_0 [simp, code abstract]: "quot_of_ratfps 0 = (0, 1)"
by transfer simp_all

lemma quot_of_ratfps_1 [simp, code abstract]: "quot_of_ratfps 1 = (1, 1)"
by transfer simp_all

lift_definition ratfps_of_poly :: "'a poly ⇒ 'a ratfps" is
"to_fract :: 'a poly ⇒ _"
by transfer simp

lemma ratfps_of_poly_code [code abstract]:
"quot_of_ratfps (ratfps_of_poly p) = (p, 1)"
by transfer' simp

lemmas zero_ratfps_code = quot_of_ratfps_0

lemmas one_ratfps_code = quot_of_ratfps_1

lemma uminus_ratfps_code [code abstract]:
"quot_of_ratfps (- x) = (let (a, b) = quot_of_ratfps x in (-a, b))"
by transfer (rule quot_of_fract_uminus)

lemma plus_ratfps_code [code abstract]:
"quot_of_ratfps (x + y) =
(let (a,b) = quot_of_ratfps x; (c,d) = quot_of_ratfps y
in  normalize_quot (a * d + b * c, b * d))"
by transfer' (rule quot_of_fract_add)

lemma minus_ratfps_code [code abstract]:
"quot_of_ratfps (x - y) =
(let (a,b) = quot_of_ratfps x; (c,d) = quot_of_ratfps y
in  normalize_quot (a * d - b * c, b * d))"
by transfer' (rule quot_of_fract_diff)

definition ratfps_cutoff :: "nat ⇒ 'a :: field_gcd ratfps ⇒ 'a poly" where
"ratfps_cutoff n x = poly_of_list (map (ratfps_nth x) [0..<n])"

definition ratfps_shift :: "nat ⇒ 'a :: field_gcd ratfps ⇒ 'a ratfps" where
"ratfps_shift n x = (let (a, b) = quot_of_ratfps (x - ratfps_of_poly (ratfps_cutoff n x))
in  quot_to_ratfps (poly_shift n a, b))"

lemma times_ratfps_code [code abstract]:
"quot_of_ratfps (x * y) =
(let (a,b) = quot_of_ratfps x; (c,d) = quot_of_ratfps y;
(e,f) = normalize_quot (a,d); (g,h) = normalize_quot (c,b)
in  (e*g, f*h))"
by transfer' (rule quot_of_fract_mult)

lemma ratfps_nth_code [code]:
"ratfps_nth x n =
(let (a,b) = quot_of_ratfps x
in  ∑i = 0..n. coeff a i * ratfps_nth_aux b (n - i))"
by transfer' simp

lemma ratfps_subdegree_code [code]:
"ratfps_subdegree x = poly_subdegree (fst (quot_of_ratfps x))"
by transfer simp

end

instantiation ratfps :: ("field_gcd") inverse
begin

lift_definition inverse_ratfps :: "'a ratfps ⇒ 'a ratfps" is
"λx. let (a,b) = quot_of_fract x
in  if coeff a 0 = 0 then 0 else inverse x"
by (auto simp: case_prod_unfold Let_def quot_of_fract_inverse)

lift_definition divide_ratfps :: "'a ratfps ⇒ 'a ratfps ⇒ 'a ratfps" is
"λf g. (if g = 0 then 0 else
let n = ratfps_subdegree g; h = ratfps_shift n g
in  ratfps_shift n (f * inverse h))" .

instance ..
end

lemma ratfps_inverse_code [code abstract]:
"quot_of_ratfps (inverse x) =
(let (a,b) = quot_of_ratfps x
in  if coeff a 0 = 0 then (0, 1)
else let u = unit_factor a in (b div u, a div u))"
by transfer' (simp_all add: Let_def case_prod_unfold quot_of_fract_inverse)

instantiation ratfps :: (equal) equal
begin

definition equal_ratfps :: "'a ratfps ⇒ 'a ratfps ⇒ bool" where
[simp]: "equal_ratfps x y ⟷ x = y"

instance by standard simp

end

lemma quot_of_fract_eq_iff [simp]: "quot_of_fract x = quot_of_fract y ⟷ x = y"
by transfer (auto simp: normalize_quot_eq_iff)

lemma equal_ratfps_code [code]: "HOL.equal x y ⟷ quot_of_ratfps x = quot_of_ratfps y"
unfolding equal_ratfps_def by transfer simp

lemma fps_of_poly_quot_normalize_quot [simp]:
"fps_of_poly (fst (normalize_quot x)) / fps_of_poly (snd (normalize_quot x)) =
fps_of_poly (fst x) / fps_of_poly (snd x)"
if "(snd x :: 'a :: field_gcd poly) ≠ 0"
proof -
from that obtain d where "fst x = fst (normalize_quot x) * d"
and "snd x = snd (normalize_quot x) * d" and "d ≠ 0"
by (rule normalize_quotE')
then show ?thesis
by (simp add: fps_of_poly_mult)
qed

lemma fps_of_poly_quot_normalize_quot' [simp]:
"fps_of_poly (fst (normalize_quot x)) / fps_of_poly (snd (normalize_quot x)) =
fps_of_poly (fst x) / fps_of_poly (snd x)"
if "coeff (snd x) 0 ≠ (0 :: 'a :: field_gcd)"
using that by (auto intro: fps_of_poly_quot_normalize_quot)

lift_definition fps_of_ratfps :: "'a :: field_gcd ratfps ⇒ 'a fps" is
"λx. fps_of_poly (numerator x) / fps_of_poly (denominator x)" .

lemma fps_of_ratfps_altdef:
"fps_of_ratfps x = (case quot_of_ratfps x of (a, b) ⇒ fps_of_poly a / fps_of_poly b)"
by transfer (simp add: case_prod_unfold)

lemma fps_of_ratfps_ratfps_of_poly [simp]: "fps_of_ratfps (ratfps_of_poly p) = fps_of_poly p"
by transfer simp

lemma fps_of_ratfps_0 [simp]: "fps_of_ratfps 0 = 0"
by transfer simp

lemma fps_of_ratfps_1 [simp]: "fps_of_ratfps 1 = 1"
by transfer simp

lemma fps_of_ratfps_uminus [simp]: "fps_of_ratfps (-x) = - fps_of_ratfps x"
by transfer (simp add: quot_of_fract_uminus case_prod_unfold Let_def fps_of_poly_simps dvd_neg_div)

lemma fps_of_ratfps_add [simp]: "fps_of_ratfps (x + y) = fps_of_ratfps x + fps_of_ratfps y"
by transfer (simp add: quot_of_fract_add Let_def case_prod_unfold fps_of_poly_simps)

lemma fps_of_ratfps_diff [simp]: "fps_of_ratfps (x - y) = fps_of_ratfps x - fps_of_ratfps y"
by transfer (simp add: quot_of_fract_diff Let_def case_prod_unfold fps_of_poly_simps)

lemma is_unit_div_div_commute: "is_unit b ⟹ is_unit c ⟹ a div b div c = a div c div b"
by (metis is_unit_div_mult2_eq mult.commute)

lemma fps_of_ratfps_mult [simp]: "fps_of_ratfps (x * y) = fps_of_ratfps x * fps_of_ratfps y"
proof (transfer, goal_cases)
case (1 x y)
moreover define x' y' where "x' = quot_of_fract x" and "y' = quot_of_fract y"
ultimately have assms: "coeff (snd x') 0 ≠ 0" "coeff (snd y') 0 ≠ 0"
by simp_all
moreover define w z where "w = normalize_quot (fst x', snd y')" and "z = normalize_quot (fst y', snd x')"
ultimately have unit: "coeff (snd x') 0 ≠ 0" "coeff (snd y') 0 ≠ 0"
"coeff (snd w) 0 ≠ 0" "coeff (snd z) 0 ≠ 0"
by (simp_all add: coeff_0_normalize_quot_nonzero)
have "fps_of_poly (fst w * fst z) / fps_of_poly (snd w * snd z) =
(fps_of_poly (fst w) / fps_of_poly (snd w)) *
(fps_of_poly (fst z) / fps_of_poly (snd z))" (is "_ = ?A * ?B")
by (simp add: is_unit_div_mult2_eq fps_of_poly_mult unit_div_mult_swap unit_div_commute unit)
also have "… = (fps_of_poly (fst x') / fps_of_poly (snd x')) *
(fps_of_poly (fst y') / fps_of_poly (snd y'))" using unit
by (simp add: w_def z_def unit_div_commute unit_div_mult_swap is_unit_div_div_commute)
finally show ?case
by (simp add: w_def z_def x'_def y'_def Let_def case_prod_unfold quot_of_fract_mult mult_ac)
qed

lemma div_const_unit_poly: "is_unit c ⟹ p div [:c:] = smult (1 div c) p"
by (simp add: is_unit_const_poly_iff unit_eq_div1)

lemma normalize_field:
"normalize (x :: 'a :: {normalization_semidom,field}) = (if x = 0 then 0 else 1)"
by (auto simp: normalize_1_iff dvd_field_iff)

lemma unit_factor_field [simp]:
"unit_factor (x :: 'a :: {normalization_semidom,field}) = x"
using unit_factor_mult_normalize[of x] normalize_field[of x]
by (simp split: if_splits)

lemma fps_of_poly_normalize_field:
"fps_of_poly (normalize (p :: 'a :: {field, normalization_semidom} poly)) =
fps_of_poly p * fps_const (inverse (lead_coeff p))"
by (cases "p = 0")
(simp_all add: normalize_poly_def div_const_unit_poly divide_simps dvd_field_iff)

lemma unit_factor_poly_altdef: "unit_factor p = monom (unit_factor (lead_coeff p)) 0"
by (simp add: unit_factor_poly_def monom_altdef)

lemma div_const_poly: "p div [:c::'a::field:] = smult (inverse c) p"
by (cases "c = 0") (simp_all add: unit_eq_div1 is_unit_triv)

lemma fps_of_ratfps_inverse [simp]: "fps_of_ratfps (inverse x) = inverse (fps_of_ratfps x)"
proof (transfer, goal_cases)
case (1 x)
hence "smult (lead_coeff (fst (quot_of_fract x))) (snd (quot_of_fract x)) div
unit_factor (fst (quot_of_fract x)) = snd (quot_of_fract x)"
if "fst (quot_of_fract x) ≠ 0" using that
by (simp add: unit_factor_poly_altdef monom_0 div_const_poly)
with 1 show ?case
by (auto simp: Let_def case_prod_unfold fps_divide_unit fps_inverse_mult
quot_of_fract_inverse mult_ac
fps_of_poly_simps fps_const_inverse
fps_of_poly_normalize_field div_smult_left [symmetric])
qed

context
includes fps_notation
begin

lemma ratfps_nth_altdef: "ratfps_nth x n = fps_of_ratfps x $n" by transfer (simp_all add: case_prod_unfold fps_divide_unit fps_times_def fps_inverse_def ratfps_nth_aux_correct Let_def) lemma fps_of_ratfps_is_unit: "fps_of_ratfps a$ 0 ≠ 0 ⟷ ratfps_nth a 0 ≠ 0"
by (simp add: ratfps_nth_altdef)

lemma ratfps_nth_0 [simp]: "ratfps_nth 0 n = 0"
by (simp add: ratfps_nth_altdef)

lemma fps_of_ratfps_cases:
obtains p q where "coeff q 0 ≠ 0" "fps_of_ratfps f = fps_of_poly p / fps_of_poly q"
by (rule that[of "snd (quot_of_ratfps f)" "fst (quot_of_ratfps f)"])
(simp_all add: fps_of_ratfps_altdef case_prod_unfold)

lemma fps_of_ratfps_cutoff [simp]:
"fps_of_poly (ratfps_cutoff n x) = fps_cutoff n (fps_of_ratfps x)"
by (simp add: fps_eq_iff ratfps_cutoff_def nth_default_def ratfps_nth_altdef)

lemma subdegree_fps_of_ratfps:
"subdegree (fps_of_ratfps x) = ratfps_subdegree x"
by transfer (simp_all add: case_prod_unfold subdegree_div_unit poly_subdegree_def)

lemma ratfps_subdegree_altdef:
"ratfps_subdegree x = subdegree (fps_of_ratfps x)"
using subdegree_fps_of_ratfps ..

end

code_datatype fps_of_ratfps

lemma fps_zero_code [code]: "0 = fps_of_ratfps 0" by simp

lemma fps_one_code [code]: "1 = fps_of_ratfps 1" by simp

lemma fps_const_code [code]: "fps_const c = fps_of_poly [:c:]" by simp

lemma fps_of_poly_code [code]: "fps_of_poly p = fps_of_ratfps (ratfps_of_poly p)" by simp

lemma fps_X_code [code]: "fps_X = fps_of_ratfps (ratfps_of_poly [:0,1:])" by simp

lemma fps_nth_code [code]: "fps_nth (fps_of_ratfps x) n = ratfps_nth x n"
by (simp add: ratfps_nth_altdef)

lemma fps_uminus_code [code]: "- fps_of_ratfps x = fps_of_ratfps (-x)" by simp

lemma fps_add_code [code]: "fps_of_ratfps x + fps_of_ratfps y = fps_of_ratfps (x + y)" by simp

lemma fps_diff_code [code]: "fps_of_ratfps x - fps_of_ratfps y = fps_of_ratfps (x - y)" by simp

lemma fps_mult_code [code]: "fps_of_ratfps x * fps_of_ratfps y = fps_of_ratfps (x * y)" by simp

lemma fps_inverse_code [code]: "inverse (fps_of_ratfps x) = fps_of_ratfps (inverse x)"
by simp

lemma fps_cutoff_code [code]: "fps_cutoff n (fps_of_ratfps x) = fps_of_poly (ratfps_cutoff n x)"
by simp

lemmas subdegree_code [code] = subdegree_fps_of_ratfps

lemma fractrel_normalize_quot:
"fractrel p p ⟹ fractrel q q ⟹
fractrel (normalize_quot p) (normalize_quot q) ⟷ fractrel p q"
by (subst fractrel_normalize_quot_left fractrel_normalize_quot_right, simp)+ (rule refl)

lemma fps_of_ratfps_eq_iff [simp]:
"fps_of_ratfps p = fps_of_ratfps q ⟷ p = q"
proof -
{
fix p q :: "'a poly fract"
assume "fractrel (quot_of_fract p) (quot_of_fract q)"
hence "p = q" by transfer (simp only: fractrel_normalize_quot)
} note A = this
show ?thesis
by transfer (auto simp: case_prod_unfold unit_eq_div1 unit_eq_div2 unit_div_commute intro: A)
qed

lemma fps_of_ratfps_eq_zero_iff [simp]:
"fps_of_ratfps p = 0 ⟷ p = 0"
by (simp del: fps_of_ratfps_0 add: fps_of_ratfps_0 [symmetric])

lemma unit_factor_snd_quot_of_ratfps [simp]:
"unit_factor (snd (quot_of_ratfps x)) = 1"
by transfer simp

lemma poly_shift_times_monom_le:
"n ≤ m ⟹ poly_shift n (monom c m * p) = monom c (m - n) * p"
by (intro poly_eqI) (auto simp: coeff_monom_mult coeff_poly_shift)

lemma poly_shift_times_monom_ge:
"n ≥ m ⟹ poly_shift n (monom c m * p) = smult c (poly_shift (n - m) p)"
by (intro poly_eqI) (auto simp: coeff_monom_mult coeff_poly_shift)

lemma poly_shift_times_monom:
"poly_shift n (monom c n * p) = smult c p"
by (intro poly_eqI) (auto simp: coeff_monom_mult coeff_poly_shift)

lemma monom_times_poly_shift:
assumes "poly_subdegree p ≥ n"
shows   "monom c n * poly_shift n p = smult c p" (is "?lhs = ?rhs")
proof (intro poly_eqI)
fix k
show "coeff ?lhs k = coeff ?rhs k"
proof (cases "k < n")
case True
with assms have "k < poly_subdegree p" by simp
hence "coeff p k = 0" by (simp add: coeff_less_poly_subdegree)
thus ?thesis by (auto simp: coeff_monom_mult coeff_poly_shift)
qed (auto simp: coeff_monom_mult coeff_poly_shift)
qed

lemma monom_times_poly_shift':
assumes "poly_subdegree p ≥ n"
shows   "monom (1 :: 'a :: comm_semiring_1) n * poly_shift n p = p"
by (simp add: monom_times_poly_shift[OF assms])

lemma subdegree_minus_cutoff_ge:
assumes "f - fps_cutoff n (f :: 'a :: ab_group_add fps) ≠ 0"
shows   "subdegree (f - fps_cutoff n f) ≥ n"
using assms by (rule subdegree_geI) simp_all

lemma fps_shift_times_X_power'': "fps_shift n (fps_X ^ n * f :: 'a :: comm_ring_1 fps) = f"
using fps_shift_times_fps_X_power'[of n f] by (simp add: mult.commute)

lemma
ratfps_shift_code [code abstract]:
"quot_of_ratfps (ratfps_shift n x) =
(let (a, b) = quot_of_ratfps (x - ratfps_of_poly (ratfps_cutoff n x))
in  (poly_shift n a, b))" (is "?lhs1 = ?rhs1") and
fps_of_ratfps_shift [simp]:
"fps_of_ratfps (ratfps_shift n x) = fps_shift n (fps_of_ratfps x)"
proof -
include fps_notation
define x' where "x' = ratfps_of_poly (ratfps_cutoff n x)"
define y where "y = quot_of_ratfps (x - x')"

have "coprime (fst y) (snd y)" unfolding y_def
by transfer (rule coprime_quot_of_fract)
also have fst_y: "fst y = monom 1 n * poly_shift n (fst y)"
proof (cases "x = x'")
case False
have "poly_subdegree (fst y) = subdegree (fps_of_poly (fst y))"
by (simp add: poly_subdegree_def)
also have "… = subdegree (fps_of_poly (fst y) / fps_of_poly (snd y))"
by (subst subdegree_div_unit) (simp_all add: y_def)
also have "fps_of_poly (fst y) / fps_of_poly (snd y) = fps_of_ratfps (x - x')"
unfolding y_def by transfer (simp add: case_prod_unfold)
also from False have "subdegree … ≥ n"
proof (intro subdegree_geI)
fix k assume "k < n"
thus "fps_of_ratfps (x - x') $k = 0" by (simp add: x'_def) qed simp_all finally show ?thesis by (rule monom_times_poly_shift' [symmetric]) qed (simp_all add: y_def) finally have coprime: "coprime (poly_shift n (fst y)) (snd y)" by simp have "quot_of_ratfps (ratfps_shift n x) = quot_of_ratfps (quot_to_ratfps (poly_shift n (fst y), snd y))" by (simp add: ratfps_shift_def Let_def case_prod_unfold x'_def y_def) also from coprime have "… = (poly_shift n (fst y), snd y)" by (intro quot_of_ratfps_quot_to_ratfps) (simp_all add: y_def normalized_fracts_def) also have "… = ?rhs1" by (simp add: case_prod_unfold Let_def y_def x'_def) finally show eq: "?lhs1 = ?rhs1" . have "fps_shift n (fps_of_ratfps x) = fps_shift n (fps_of_ratfps (x - x'))" by (intro fps_ext) (simp_all add: x'_def) also have "fps_of_ratfps (x - x') = fps_of_poly (fst y) / fps_of_poly (snd y)" by (simp add: fps_of_ratfps_altdef y_def case_prod_unfold) also have "fps_shift n … = fps_of_ratfps (ratfps_shift n x)" by (subst fst_y, subst fps_of_poly_mult, subst unit_div_mult_swap [symmetric]) (simp_all add: y_def fps_of_poly_monom fps_shift_times_X_power'' eq fps_of_ratfps_altdef case_prod_unfold Let_def x'_def) finally show "fps_of_ratfps (ratfps_shift n x) = fps_shift n (fps_of_ratfps x)" .. qed lemma fps_shift_code [code]: "fps_shift n (fps_of_ratfps x) = fps_of_ratfps (ratfps_shift n x)" by simp instantiation fps :: (equal) equal begin definition equal_fps :: "'a fps ⇒ 'a fps ⇒ bool" where [simp]: "equal_fps f g ⟷ f = g" instance by standard simp_all end lemma equal_fps_code [code]: "HOL.equal (fps_of_ratfps f) (fps_of_ratfps g) ⟷ f = g" by simp lemma fps_of_ratfps_divide [simp]: "fps_of_ratfps (f div g) = fps_of_ratfps f div fps_of_ratfps g" unfolding fps_divide_def Let_def by transfer' (simp add: Let_def ratfps_subdegree_altdef) lemma ratfps_eqI: "fps_of_ratfps x = fps_of_ratfps y ⟹ x = y" by simp instance ratfps :: ("field_gcd") algebraic_semidom by standard (auto intro: ratfps_eqI) lemma fps_of_ratfps_dvd [simp]: "fps_of_ratfps x dvd fps_of_ratfps y ⟷ x dvd y" proof assume "fps_of_ratfps x dvd fps_of_ratfps y" hence "fps_of_ratfps y = fps_of_ratfps y div fps_of_ratfps x * fps_of_ratfps x" by simp also have "… = fps_of_ratfps (y div x * x)" by simp finally have "y = y div x * x" by (subst (asm) fps_of_ratfps_eq_iff) thus "x dvd y" by (intro dvdI[of _ _ "y div x"]) (simp add: mult_ac) next assume "x dvd y" hence "y = y div x * x" by simp also have "fps_of_ratfps … = fps_of_ratfps (y div x) * fps_of_ratfps x" by simp finally show "fps_of_ratfps x dvd fps_of_ratfps y" by (simp del: fps_of_ratfps_divide) qed lemma is_unit_ratfps_iff [simp]: "is_unit x ⟷ ratfps_nth x 0 ≠ 0" proof assume "is_unit x" then obtain y where "1 = x * y" by (auto elim!: dvdE) hence "1 = fps_of_ratfps (x * y)" by (simp del: fps_of_ratfps_mult) also have "… = fps_of_ratfps x * fps_of_ratfps y" by simp finally have "is_unit (fps_of_ratfps x)" by (rule dvdI[of _ _ "fps_of_ratfps y"]) thus "ratfps_nth x 0 ≠ 0" by (simp add: ratfps_nth_altdef) next assume "ratfps_nth x 0 ≠ 0" hence "fps_of_ratfps (x * inverse x) = 1" by (simp add: ratfps_nth_altdef inverse_mult_eq_1') also have "… = fps_of_ratfps 1" by simp finally have "x * inverse x = 1" by (subst (asm) fps_of_ratfps_eq_iff) thus "is_unit x" by (intro dvdI[of _ _ "inverse x"]) simp_all qed instantiation ratfps :: ("field_gcd") normalization_semidom begin definition unit_factor_ratfps :: "'a ratfps ⇒ 'a ratfps" where "unit_factor x = ratfps_shift (ratfps_subdegree x) x" definition normalize_ratfps :: "'a ratfps ⇒ 'a ratfps" where "normalize x = (if x = 0 then 0 else ratfps_of_poly (monom 1 (ratfps_subdegree x)))" lemma fps_of_ratfps_unit_factor [simp]: "fps_of_ratfps (unit_factor x) = unit_factor (fps_of_ratfps x)" unfolding unit_factor_ratfps_def by (simp add: ratfps_subdegree_altdef) lemma fps_of_ratfps_normalize [simp]: "fps_of_ratfps (normalize x) = normalize (fps_of_ratfps x)" unfolding normalize_ratfps_def by (simp add: fps_of_poly_monom ratfps_subdegree_altdef) instance proof show "unit_factor x * normalize x = x" "normalize (0 :: 'a ratfps) = 0" "unit_factor (0 :: 'a ratfps) = 0" for x :: "'a ratfps" by (rule ratfps_eqI, simp add: ratfps_subdegree_code del: fps_of_ratfps_eq_iff fps_unit_factor_def fps_normalize_def)+ show "is_unit (unit_factor a)" if "a ≠ 0" for a :: "'a ratfps" using that by (auto simp: ratfps_nth_altdef) fix a b :: "'a ratfps" assume "is_unit a" thus "unit_factor (a * b) = a * unit_factor b" by (intro ratfps_eqI, unfold fps_of_ratfps_unit_factor fps_of_ratfps_mult, subst unit_factor_mult_unit_left) (auto simp: ratfps_nth_altdef) show "unit_factor a = a" if "is_unit a" for a :: "'a ratfps" by (rule ratfps_eqI) (insert that, auto simp: fps_of_ratfps_is_unit) qed end instance ratfps :: ("field_gcd") normalization_semidom_multiplicative proof show "unit_factor (a * b) = unit_factor a * unit_factor b" for a b :: "'a ratfps" by (rule ratfps_eqI, insert unit_factor_mult[of "fps_of_ratfps a" "fps_of_ratfps b"]) (simp del: fps_of_ratfps_eq_iff) qed instantiation ratfps :: ("field_gcd") semidom_modulo begin lift_definition modulo_ratfps :: "'a ratfps ⇒ 'a ratfps ⇒ 'a ratfps" is "λf g. if g = 0 then f else let n = ratfps_subdegree g; h = ratfps_shift n g in ratfps_of_poly (ratfps_cutoff n (f * inverse h)) * h" . lemma fps_of_ratfps_mod [simp]: "fps_of_ratfps (f mod g :: 'a ratfps) = fps_of_ratfps f mod fps_of_ratfps g" unfolding fps_mod_def by transfer' (simp add: Let_def ratfps_subdegree_altdef) instance by standard (auto intro: ratfps_eqI) end instantiation ratfps :: ("field_gcd") euclidean_ring begin definition euclidean_size_ratfps :: "'a ratfps ⇒ nat" where "euclidean_size_ratfps x = (if x = 0 then 0 else 2 ^ ratfps_subdegree x)" lemma fps_of_ratfps_euclidean_size [simp]: "euclidean_size x = euclidean_size (fps_of_ratfps x)" unfolding euclidean_size_ratfps_def fps_euclidean_size_def by (simp add: ratfps_subdegree_altdef) instance proof show "euclidean_size (0 :: 'a ratfps) = 0" by simp show "euclidean_size (a mod b) < euclidean_size b" "euclidean_size a ≤ euclidean_size (a * b)" if "b ≠ 0" for a b :: "'a ratfps" using that by (simp_all add: mod_size_less size_mult_mono) qed end instantiation ratfps :: ("field_gcd") euclidean_ring_cancel begin instance by standard (auto intro: ratfps_eqI) end lemma quot_of_ratfps_eq_iff [simp]: "quot_of_ratfps x = quot_of_ratfps y ⟷ x = y" by transfer simp lemma ratfps_eq_0_code: "x = 0 ⟷ fst (quot_of_ratfps x) = 0" proof assume "fst (quot_of_ratfps x) = 0" moreover have "coprime (fst (quot_of_ratfps x)) (snd (quot_of_ratfps x))" by transfer (simp add: coprime_quot_of_fract) moreover have "normalize (snd (quot_of_ratfps x)) = snd (quot_of_ratfps x)" by (simp add: div_unit_factor [symmetric] del: div_unit_factor) ultimately have "quot_of_ratfps x = (0,1)" by (simp add: prod_eq_iff normalize_idem_imp_is_unit_iff) also have "… = quot_of_ratfps 0" by simp finally show "x = 0" by (subst (asm) quot_of_ratfps_eq_iff) qed simp_all lemma fps_dvd_code [code_unfold]: "x dvd y ⟷ y = 0 ∨ ((x::'a::field_gcd fps) ≠ 0 ∧ subdegree x ≤ subdegree y)" using fps_dvd_iff[of x y] by (cases "x = 0") auto lemma ratfps_dvd_code [code_unfold]: "x dvd y ⟷ y = 0 ∨ (x ≠ 0 ∧ ratfps_subdegree x ≤ ratfps_subdegree y)" using fps_dvd_code [of "fps_of_ratfps x" "fps_of_ratfps y"] by (simp add: ratfps_subdegree_altdef) instance ratfps :: ("field_gcd") normalization_euclidean_semiring .. instantiation ratfps :: ("field_gcd") euclidean_ring_gcd begin definition "gcd_ratfps = (Euclidean_Algorithm.gcd :: 'a ratfps ⇒ _)" definition "lcm_ratfps = (Euclidean_Algorithm.lcm :: 'a ratfps ⇒ _)" definition "Gcd_ratfps = (Euclidean_Algorithm.Gcd :: 'a ratfps set ⇒ _)" definition "Lcm_ratfps = (Euclidean_Algorithm.Lcm:: 'a ratfps set ⇒ _)" instance by standard (simp_all add: gcd_ratfps_def lcm_ratfps_def Gcd_ratfps_def Lcm_ratfps_def) end lemma ratfps_eq_0_iff: "x = 0 ⟷ fps_of_ratfps x = 0" using fps_of_ratfps_eq_iff[of x 0] unfolding fps_of_ratfps_0 by simp lemma ratfps_of_poly_eq_0_iff: "ratfps_of_poly x = 0 ⟷ x = 0" by (auto simp: ratfps_eq_0_iff) lemma ratfps_gcd: assumes [simp]: "f ≠ 0" "g ≠ 0" shows "gcd f g = ratfps_of_poly (monom 1 (min (ratfps_subdegree f) (ratfps_subdegree g)))" by (rule sym, rule gcdI) (auto simp: ratfps_subdegree_altdef ratfps_dvd_code subdegree_fps_of_poly ratfps_of_poly_eq_0_iff normalize_ratfps_def) lemma ratfps_gcd_altdef: "gcd (f :: 'a :: field_gcd ratfps) g = (if f = 0 ∧ g = 0 then 0 else if f = 0 then ratfps_of_poly (monom 1 (ratfps_subdegree g)) else if g = 0 then ratfps_of_poly (monom 1 (ratfps_subdegree f)) else ratfps_of_poly (monom 1 (min (ratfps_subdegree f) (ratfps_subdegree g))))" by (simp add: ratfps_gcd normalize_ratfps_def) lemma ratfps_lcm: assumes [simp]: "f ≠ 0" "g ≠ 0" shows "lcm f g = ratfps_of_poly (monom 1 (max (ratfps_subdegree f) (ratfps_subdegree g)))" by (rule sym, rule lcmI) (auto simp: ratfps_subdegree_altdef ratfps_dvd_code subdegree_fps_of_poly ratfps_of_poly_eq_0_iff normalize_ratfps_def) lemma ratfps_lcm_altdef: "lcm (f :: 'a :: field_gcd ratfps) g = (if f = 0 ∨ g = 0 then 0 else ratfps_of_poly (monom 1 (max (ratfps_subdegree f) (ratfps_subdegree g))))" by (simp add: ratfps_lcm) lemma ratfps_Gcd: assumes "A - {0} ≠ {}" shows "Gcd A = ratfps_of_poly (monom 1 (INF f∈A-{0}. ratfps_subdegree f))" proof (rule sym, rule GcdI) fix f assume "f ∈ A" thus "ratfps_of_poly (monom 1 (INF f∈A - {0}. ratfps_subdegree f)) dvd f" by (cases "f = 0") (auto simp: ratfps_dvd_code ratfps_of_poly_eq_0_iff ratfps_subdegree_altdef subdegree_fps_of_poly intro!: cINF_lower) next fix d assume d: "⋀f. f ∈ A ⟹ d dvd f" from assms obtain f where "f ∈ A - {0}" by auto with d[of f] have [simp]: "d ≠ 0" by auto from d assms have "ratfps_subdegree d ≤ (INF f∈A-{0}. ratfps_subdegree f)" by (intro cINF_greatest) (auto simp: ratfps_dvd_code) with d assms show "d dvd ratfps_of_poly (monom 1 (INF f∈A-{0}. ratfps_subdegree f))" by (simp add: ratfps_dvd_code ratfps_subdegree_altdef subdegree_fps_of_poly) qed (simp_all add: ratfps_subdegree_altdef subdegree_fps_of_poly normalize_ratfps_def) lemma ratfps_Gcd_altdef: "Gcd (A :: 'a :: field_gcd ratfps set) = (if A ⊆ {0} then 0 else ratfps_of_poly (monom 1 (INF f∈A-{0}. ratfps_subdegree f)))" using ratfps_Gcd by auto lemma ratfps_Lcm: assumes "A ≠ {}" "0 ∉ A" "bdd_above (ratfps_subdegreeA)" shows "Lcm A = ratfps_of_poly (monom 1 (SUP f∈A. ratfps_subdegree f))" proof (rule sym, rule LcmI) fix f assume "f ∈ A" moreover from assms(3) have "bdd_above (ratfps_subdegree  A)" by auto ultimately show "f dvd ratfps_of_poly (monom 1 (SUP f∈A. ratfps_subdegree f))" using assms(2) by (cases "f = 0") (auto simp: ratfps_dvd_code ratfps_of_poly_eq_0_iff subdegree_fps_of_poly ratfps_subdegree_altdef [abs_def] intro!: cSUP_upper) next fix d assume d: "⋀f. f ∈ A ⟹ f dvd d" from assms obtain f where f: "f ∈ A" "f ≠ 0" by auto show "ratfps_of_poly (monom 1 (SUP f∈A. ratfps_subdegree f)) dvd d" proof (cases "d = 0") assume "d ≠ 0" moreover from d have "⋀f. f ∈ A ⟹ f ≠ 0 ⟹ f dvd d" by blast ultimately have "ratfps_subdegree d ≥ (SUP f∈A. ratfps_subdegree f)" using assms by (intro cSUP_least) (auto simp: ratfps_dvd_code) with ‹d ≠ 0› show ?thesis by (simp add: ratfps_dvd_code ratfps_of_poly_eq_0_iff ratfps_subdegree_altdef subdegree_fps_of_poly) qed simp_all qed (simp_all add: ratfps_subdegree_altdef subdegree_fps_of_poly normalize_ratfps_def) lemma ratfps_Lcm_altdef: "Lcm (A :: 'a :: field_gcd ratfps set) = (if 0 ∈ A ∨ ¬bdd_above (ratfps_subdegreeA) then 0 else if A = {} then 1 else ratfps_of_poly (monom 1 (SUP f∈A. ratfps_subdegree f)))" proof (cases "bdd_above (ratfps_subdegreeA)") assume unbounded: "¬bdd_above (ratfps_subdegreeA)" have "Lcm A = 0" proof (rule ccontr) assume "Lcm A ≠ 0" from unbounded obtain f where f: "f ∈ A" "ratfps_subdegree (Lcm A) < ratfps_subdegree f" unfolding bdd_above_def by (auto simp: not_le) moreover from this and ‹Lcm A ≠ 0› have "ratfps_subdegree f ≤ ratfps_subdegree (Lcm A)" using dvd_Lcm[of f A] by (auto simp: ratfps_dvd_code) ultimately show False by simp qed with unbounded show ?thesis by simp qed (simp_all add: ratfps_Lcm Lcm_eq_0_I) lemma fps_of_ratfps_quot_to_ratfps: "coeff y 0 ≠ 0 ⟹ fps_of_ratfps (quot_to_ratfps (x,y)) = fps_of_poly x / fps_of_poly y" proof (transfer, goal_cases) case (1 y x) define x' y' where "x' = fst (normalize_quot (x,y))" and "y' = snd (normalize_quot (x,y))" from 1 have nz: "y ≠ 0" by auto have eq: "normalize_quot (x', y') = (x', y')" by (simp add: x'_def y'_def) from normalize_quotE[OF nz, of x] guess d . note d = this [folded x'_def y'_def] have "(case quot_of_fract (if coeff y' 0 = 0 then 0 else quot_to_fract (x', y')) of (a, b) ⇒ fps_of_poly a / fps_of_poly b) = fps_of_poly x / fps_of_poly y" using d eq 1 by (auto simp: case_prod_unfold fps_of_poly_simps quot_of_fract_quot_to_fract Let_def coeff_0_mult) thus ?case by (auto simp add: Let_def case_prod_unfold x'_def y'_def) qed lemma fps_of_ratfps_quot_to_ratfps_code_post1: "fps_of_ratfps (quot_to_ratfps (x,pCons 1 y)) = fps_of_poly x / fps_of_poly (pCons 1 y)" "fps_of_ratfps (quot_to_ratfps (x,pCons (-1) y)) = fps_of_poly x / fps_of_poly (pCons (-1) y)" by (simp_all add: fps_of_ratfps_quot_to_ratfps) lemma fps_of_ratfps_quot_to_ratfps_code_post2: "fps_of_ratfps (quot_to_ratfps (x'::'a::{field_char_0,field_gcd} poly,pCons (numeral n) y')) = fps_of_poly x' / fps_of_poly (pCons (numeral n) y')" "fps_of_ratfps (quot_to_ratfps (x'::'a::{field_char_0,field_gcd} poly,pCons (-numeral n) y')) = fps_of_poly x' / fps_of_poly (pCons (-numeral n) y')" by (simp_all add: fps_of_ratfps_quot_to_ratfps) lemmas fps_of_ratfps_quot_to_ratfps_code_post [code_post] = fps_of_ratfps_quot_to_ratfps_code_post1 fps_of_ratfps_quot_to_ratfps_code_post2 lemma fps_dehorner: fixes a b c :: "'a :: semiring_1 fps" and d e f :: "'b :: ring_1 fps" shows "(b + c) * fps_X = b * fps_X + c * fps_X" "(a * fps_X) * fps_X = a * fps_X ^ 2" "a * fps_X ^ m * fps_X = a * fps_X ^ (Suc m)" "a * fps_X * fps_X ^ m = a * fps_X ^ (Suc m)" "a * fps_X^m * fps_X^n = a * fps_X^(m+n)" "a + (b + c) = a + b + c" "a * 1 = a" "1 * a = a" "d + - e = d - e" "(-d) * e = - (d * e)" "d + (e - f) = d + e - f" "(d - e) * fps_X = d * fps_X - e * fps_X" "fps_X * fps_X = fps_X^2" "fps_X * fps_X^m = fps_X^(Suc m)" "fps_X^m * fps_X = fps_X^Suc m" "fps_X^m * fps_X^n = fps_X^(m + n)" by (simp_all add: algebra_simps power2_eq_square power_add power_commutes) lemma fps_divide_1: "(a :: 'a :: field fps) / 1 = a" by simp lemmas fps_of_poly_code_post [code_post] = fps_of_poly_simps fps_const_0_eq_0 fps_const_1_eq_1 numeral_fps_const [symmetric] fps_const_neg [symmetric] fps_const_divide [symmetric] fps_dehorner Suc_numeral arith_simps fps_divide_1 context includes term_syntax begin definition valterm_ratfps :: "'a ::{field_gcd, typerep} poly × (unit ⇒ Code_Evaluation.term) ⇒ 'a poly × (unit ⇒ Code_Evaluation.term) ⇒ 'a ratfps × (unit ⇒ Code_Evaluation.term)" where [code_unfold]: "valterm_ratfps k l = Code_Evaluation.valtermify (/) {⋅} (Code_Evaluation.valtermify ratfps_of_poly {⋅} k) {⋅} (Code_Evaluation.valtermify ratfps_of_poly {⋅} l)" end instantiation ratfps :: ("{field_gcd,random}") random begin context includes state_combinator_syntax term_syntax begin definition "Quickcheck_Random.random i = Quickcheck_Random.random i ∘→ (λnum::'a poly × (unit ⇒ term). Quickcheck_Random.random i ∘→ (λdenom::'a poly × (unit ⇒ term). Pair (let denom = (if fst denom = 0 then Code_Evaluation.valtermify 1 else denom) in valterm_ratfps num denom)))" instance .. end end instantiation ratfps :: ("{field,factorial_ring_gcd,exhaustive}") exhaustive begin definition "exhaustive_ratfps f d = Quickcheck_Exhaustive.exhaustive (λnum. Quickcheck_Exhaustive.exhaustive (λdenom. f ( let denom = if denom = 0 then 1 else denom in ratfps_of_poly num / ratfps_of_poly denom)) d) d" instance .. end instantiation ratfps :: ("{field_gcd,full_exhaustive}") full_exhaustive begin definition "full_exhaustive_ratfps f d = Quickcheck_Exhaustive.full_exhaustive (λnum::'a poly × (unit ⇒ term). Quickcheck_Exhaustive.full_exhaustive (λdenom::'a poly × (unit ⇒ term). f (let denom = if fst denom = 0 then Code_Evaluation.valtermify 1 else denom in valterm_ratfps num denom)) d) d" instance .. end quickcheck_generator fps constructors: fps_of_ratfps end  # Theory Pochhammer_Polynomials (* File: Pochhammer_Polynomials.thy Author: Manuel Eberl, TU München *) section ‹Falling factorial as a polynomial› theory Pochhammer_Polynomials imports Complex_Main "HOL-Combinatorics.Stirling" "HOL-Computational_Algebra.Polynomial" begin definition pochhammer_poly :: "nat ⇒ 'a :: {comm_semiring_1} poly" where "pochhammer_poly n = Poly [of_nat (stirling n k). k ← [0..<Suc n]]" lemma pochhammer_poly_code [code abstract]: "coeffs (pochhammer_poly n) = map of_nat (stirling_row n)" by (simp add: pochhammer_poly_def stirling_row_def Let_def) lemma coeff_pochhammer_poly: "coeff (pochhammer_poly n) k = of_nat (stirling n k)" by (simp add: pochhammer_poly_def nth_default_def del: upt_Suc) lemma degree_pochhammer_poly [simp]: "degree (pochhammer_poly n) = n" by (simp add: degree_eq_length_coeffs pochhammer_poly_def) lemma pochhammer_poly_0 [simp]: "pochhammer_poly 0 = 1" by (simp add: pochhammer_poly_def) lemma pochhammer_poly_Suc: "pochhammer_poly (Suc n) = [:of_nat n,1:] * pochhammer_poly n" by (cases "n = 0") (simp_all add: poly_eq_iff coeff_pochhammer_poly coeff_pCons split: nat.split) lemma pochhammer_poly_altdef: "pochhammer_poly n = (∏i<n. [:of_nat i,1:])" by (induction n) (simp_all add: pochhammer_poly_Suc) lemma eval_pochhammer_poly: "poly (pochhammer_poly n) k = pochhammer k n" by (cases n) (auto simp add: pochhammer_poly_altdef poly_prod add_ac lessThan_Suc_atMost pochhammer_Suc_prod atLeast0AtMost) lemma pochhammer_poly_Suc': "pochhammer_poly (Suc n) = pCons 0 (pcompose (pochhammer_poly n) [:1,1:])" by (simp add: pochhammer_poly_altdef prod.lessThan_Suc_shift pcompose_prod pcompose_pCons add_ac del: prod.lessThan_Suc) end  # Theory Linear_Recurrences_Misc (* File: Linear_Recurrences_Misc.thy Author: Manuel Eberl, TU München *) section ‹Miscellaneous material required for linear recurrences› theory Linear_Recurrences_Misc imports Complex_Main "HOL-Computational_Algebra.Computational_Algebra" "HOL-Computational_Algebra.Polynomial_Factorial" begin fun zip_with where "zip_with f (x#xs) (y#ys) = f x y # zip_with f xs ys" | "zip_with f _ _ = []" lemma length_zip_with [simp]: "length (zip_with f xs ys) = min (length xs) (length ys)" by (induction f xs ys rule: zip_with.induct) simp_all lemma zip_with_altdef: "zip_with f xs ys = map (λ(x,y). f x y) (zip xs ys)" by (induction f xs ys rule: zip_with.induct) simp_all lemma zip_with_nth [simp]: "n < length xs ⟹ n < length ys ⟹ zip_with f xs ys ! n = f (xs!n) (ys!n)" by (simp add: zip_with_altdef) lemma take_zip_with: "take n (zip_with f xs ys) = zip_with f (take n xs) (take n ys)" proof (induction f xs ys arbitrary: n rule: zip_with.induct) case (1 f x xs y ys n) thus ?case by (cases n) simp_all qed simp_all lemma drop_zip_with: "drop n (zip_with f xs ys) = zip_with f (drop n xs) (drop n ys)" proof (induction f xs ys arbitrary: n rule: zip_with.induct) case (1 f x xs y ys n) thus ?case by (cases n) simp_all qed simp_all lemma map_zip_with: "map f (zip_with g xs ys) = zip_with (λx y. f (g x y)) xs ys" by (induction g xs ys rule: zip_with.induct) simp_all lemma zip_with_map: "zip_with f (map g xs) (map h ys) = zip_with (λx y. f (g x) (h y)) xs ys" by (induction "λx y. f (g x) (h y)" xs ys rule: zip_with.induct) simp_all lemma zip_with_map_left: "zip_with f (map g xs) ys = zip_with (λx y. f (g x) y) xs ys" using zip_with_map[of f g xs "λx. x" ys] by simp lemma zip_with_map_right: "zip_with f xs (map g ys) = zip_with (λx y. f x (g y)) xs ys" using zip_with_map[of f "λx. x" xs g ys] by simp lemma zip_with_swap: "zip_with (λx y. f y x) xs ys = zip_with f ys xs" by (induction f ys xs rule: zip_with.induct) simp_all lemma set_zip_with: "set (zip_with f xs ys) = (λ(x,y). f x y)  set (zip xs ys)" by (induction f xs ys rule: zip_with.induct) simp_all lemma zip_with_Pair: "zip_with Pair (xs :: 'a list) (ys :: 'b list) = zip xs ys" by (induction "Pair :: 'a ⇒ 'b ⇒ _" xs ys rule: zip_with.induct) simp_all lemma zip_with_altdef': "zip_with f xs ys = [f (xs!i) (ys!i). i ← [0..<min (length xs) (length ys)]]" by (induction f xs ys rule: zip_with.induct) (simp_all add: map_upt_Suc del: upt_Suc) lemma zip_altdef: "zip xs ys = [(xs!i, ys!i). i ← [0..<min (length xs) (length ys)]]" by (simp add: zip_with_Pair [symmetric] zip_with_altdef') lemma card_poly_roots_bound: fixes p :: "'a::{comm_ring_1,ring_no_zero_divisors} poly" assumes "p ≠ 0" shows "card {x. poly p x = 0} ≤ degree p" using assms proof (induction "degree p" arbitrary: p rule: less_induct) case (less p) show ?case proof (cases "∃x. poly p x = 0") case False hence "{x. poly p x = 0} = {}" by blast thus ?thesis by simp next case True then obtain x where x: "poly p x = 0" by blast hence "[:-x, 1:] dvd p" by (subst (asm) poly_eq_0_iff_dvd) then obtain q where q: "p = [:-x, 1:] * q" by (auto simp: dvd_def) with ‹p ≠ 0› have [simp]: "q ≠ 0" by auto have deg: "degree p = Suc (degree q)" by (subst q, subst degree_mult_eq) auto have "card {x. poly p x = 0} ≤ card (insert x {x. poly q x = 0})" by (intro card_mono) (auto intro: poly_roots_finite simp: q) also have "… ≤ Suc (card {x. poly q x = 0})" by (rule card_insert_le_m1) auto also from deg have "card {x. poly q x = 0} ≤ degree q" using ‹p ≠ 0› and q by (intro less) auto also have "Suc … = degree p" by (simp add: deg) finally show ?thesis by - simp_all qed qed lemma poly_eqI_degree: fixes p q :: "'a :: {comm_ring_1, ring_no_zero_divisors} poly" assumes "⋀x. x ∈ A ⟹ poly p x = poly q x" assumes "card A > degree p" "card A > degree q" shows "p = q" proof (rule ccontr) assume neq: "p ≠ q" have "degree (p - q) ≤ max (degree p) (degree q)" by (rule degree_diff_le_max) also from assms have "… < card A" by linarith also have "… ≤ card {x. poly (p - q) x = 0}" using neq and assms by (intro card_mono poly_roots_finite) auto finally have "degree (p - q) < card {x. poly (p - q) x = 0}" . moreover have "degree (p - q) ≥ card {x. poly (p - q) x = 0}" using neq by (intro card_poly_roots_bound) auto ultimately show False by linarith qed lemma poly_root_order_induct [case_names 0 no_roots root]: fixes p :: "'a :: idom poly" assumes "P 0" "⋀p. (⋀x. poly p x ≠ 0) ⟹ P p" "⋀p x n. n > 0 ⟹ poly p x ≠ 0 ⟹ P p ⟹ P ([:-x, 1:] ^ n * p)" shows "P p" proof (induction "degree p" arbitrary: p rule: less_induct) case (less p) consider "p = 0" | "p ≠ 0" "∃x. poly p x = 0" | "⋀x. poly p x ≠ 0" by blast thus ?case proof cases case 3 with assms(2)[of p] show ?thesis by simp next case 2 then obtain x where x: "poly p x = 0" by auto have "[:-x, 1:] ^ order x p dvd p" by (intro order_1) then obtain q where q: "p = [:-x, 1:] ^ order x p * q" by (auto simp: dvd_def) with 2 have [simp]: "q ≠ 0" by auto have order_pos: "order x p > 0" using ‹p ≠ 0› and x by (auto simp: order_root) have "order x p = order x p + order x q" by (subst q, subst order_mult) (auto simp: order_power_n_n) hence [simp]: "order x q = 0" by simp have deg: "degree p = order x p + degree q" by (subst q, subst degree_mult_eq) (auto simp: degree_power_eq) with order_pos have "degree q < degree p" by simp hence "P q" by (rule less) with order_pos have "P ([:-x, 1:] ^ order x p * q)" by (intro assms(3)) (auto simp: order_root) with q show ?thesis by simp qed (simp_all add: assms(1)) qed lemma complex_poly_decompose: "smult (lead_coeff p) (∏z|poly p z = 0. [:-z, 1:] ^ order z p) = (p :: complex poly)" proof (induction p rule: poly_root_order_induct) case (no_roots p) show ?case proof (cases "degree p = 0") case False hence "¬constant (poly p)" by (subst constant_degree) with fundamental_theorem_of_algebra and no_roots show ?thesis by blast qed (auto elim!: degree_eq_zeroE) next case (root p x n) from root have *: "{z. poly ([:- x, 1:] ^ n * p) z = 0} = insert x {z. poly p z = 0}" by auto have "smult (lead_coeff ([:-x, 1:] ^ n * p)) (∏z|poly ([:-x,1:] ^ n * p) z = 0. [:-z, 1:] ^ order z ([:- x, 1:] ^ n * p)) = [:- x, 1:] ^ order x ([:- x, 1:] ^ n * p) * smult (lead_coeff p) (∏z∈{z. poly p z = 0}. [:- z, 1:] ^ order z ([:- x, 1:] ^ n * p))" by (subst *, subst prod.insert) (insert root, auto intro: poly_roots_finite simp: mult_ac lead_coeff_mult lead_coeff_power) also have "order x ([:- x, 1:] ^ n * p) = n" using root by (subst order_mult) (auto simp: order_power_n_n order_0I) also have "(∏z∈{z. poly p z = 0}. [:- z, 1:] ^ order z ([:- x, 1:] ^ n * p)) = (∏z∈{z. poly p z = 0}. [:- z, 1:] ^ order z p)" proof (intro prod.cong refl, goal_cases) case (1 y) with root have "order y ([:-x,1:] ^ n) = 0" by (intro order_0I) auto thus ?case using root by (subst order_mult) auto qed also note root.IH finally show ?case . qed simp_all lemma normalize_field: "normalize (x :: 'a :: {normalization_semidom,field}) = (if x = 0 then 0 else 1)" by (auto simp: normalize_1_iff dvd_field_iff) lemma unit_factor_field [simp]: "unit_factor (x :: 'a :: {normalization_semidom,field}) = x" using unit_factor_mult_normalize[of x] normalize_field[of x] by (simp split: if_splits) lemma coprime_linear_poly: fixes c :: "'a :: field_gcd" assumes "c ≠ c'" shows "coprime [:c,1:] [:c',1:]" proof - have "gcd [:c,1:] [:c',1:] = gcd ([:c,1:] - [:c',1:]) [:c',1:]" by (rule gcd_diff1 [symmetric]) also have "[:c,1:] - [:c',1:] = [:c-c':]" by simp also from assms have "gcd … [:c',1:] = normalize [:c-c':]" by (intro gcd_proj1_if_dvd) (auto simp: const_poly_dvd_iff dvd_field_iff) also from assms have "… = 1" by (simp add: normalize_poly_def) finally show "coprime [:c,1:] [:c',1:]" by (simp add: gcd_eq_1_imp_coprime) qed lemma coprime_linear_poly': fixes c :: "'a :: field_gcd" assumes "c ≠ c'" "c ≠ 0" "c' ≠ 0" shows "coprime [:1,c:] [:1,c':]" proof - have "gcd [:1,c:] [:1,c':] = gcd ([:1,c:] mod [:1,c':]) [:1,c':]" by simp also have "eucl_rel_poly [:1, c:] [:1, c':] ([:c/c':], [:1-c/c':])" using assms by (auto simp add: eucl_rel_poly_iff one_pCons) hence "[:1,c:] mod [:1,c':] = [:1 - c / c':]" by (rule mod_poly_eq) also from assms have "gcd … [:1,c':] = normalize ([:1 - c / c':])" by (intro gcd_proj1_if_dvd) (auto simp: const_poly_dvd_iff dvd_field_iff) also from assms have "… = 1" by (auto simp: normalize_poly_def) finally show ?thesis by (rule gcd_eq_1_imp_coprime) qed end  # Theory Partial_Fraction_Decomposition (* File: Partial_Fraction_Decomposition.thy Author: Manuel Eberl <eberlm@in.tum.de> Partial fraction decomposition on Euclidean rings, i.e. decomposing a quotient into a sum of quotients where each denominator is a power of an irreducible element. (and possibly one summand that is an entire element) The most interesting setting is when the Euclidean ring is a polynomial ring. *) section ‹Partial Fraction Decomposition› theory Partial_Fraction_Decomposition imports Main "HOL-Computational_Algebra.Computational_Algebra" "HOL-Computational_Algebra.Polynomial_Factorial" "HOL-Library.Sublist" Linear_Recurrences_Misc begin subsection ‹Decomposition on general Euclidean rings› text ‹ Consider elements$x, y_1, \ldots, y_n$of a ring$R$, where the$y_i$are pairwise coprime. A \emph{Partial Fraction Decomposition} of these elements (or rather the formal quotient$x / (y_1 \ldots y_n)$that they represent) is a finite sum of summands of the form$a / y_i ^ k$. Obviously, the sum can be arranged such that there is at most one summand with denominator$y_i ^ n$for any combination of$i$and$n$; in particular, there is at most one summand with denominator 1. We can decompose the summands further by performing division with remainder until in all quotients, the numerator's Euclidean size is less than that of the denominator. › text ‹ The following function performs the first step of the above process: it takes the values$x$and$y_1,\ldots, y_n$and returns the numerators of the summands in the decomposition. (the denominators are simply the$y_i$from the input) › fun decompose :: "('a :: euclidean_ring_gcd) ⇒ 'a list ⇒ 'a list" where "decompose x [] = []" | "decompose x [y] = [x]" | "decompose x (y#ys) = (case bezout_coefficients y (prod_list ys) of (a, b) ⇒ (b*x) # decompose (a*x) ys)" lemma decompose_rec: "ys ≠ [] ⟹ decompose x (y#ys) = (case bezout_coefficients y (prod_list ys) of (a, b) ⇒ (b*x) # decompose (a*x) ys)" by (cases ys) simp_all lemma length_decompose [simp]: "length (decompose x ys) = length ys" proof (induction x ys rule: decompose.induct) case (3 x y z ys) obtain a b where ab: "(a,b) = bezout_coefficients y (prod_list (z#ys))" by (cases "bezout_coefficients y (z * prod_list ys)") simp_all from 3[OF ab] ab[symmetric] show ?case by simp qed simp_all fun decompose' :: "('a :: euclidean_ring_gcd) ⇒ 'a list ⇒ 'a list ⇒ 'a list" where "decompose' x [] _ = []" | "decompose' x [y] _ = [x]" | "decompose' _ _ [] = []" | "decompose' x (y#ys) (p#ps) = (case bezout_coefficients y p of (a, b) ⇒ (b*x) # decompose' (a*x) ys ps)" primrec decompose_aux :: "'a :: {ab_semigroup_mult,monoid_mult} ⇒ _" where "decompose_aux acc [] = [acc]" | "decompose_aux acc (x#xs) = acc # decompose_aux (x * acc) xs" lemma decompose_code [code]: "decompose x ys = decompose' x ys (tl (rev (decompose_aux 1 (rev ys))))" proof (induction x ys rule: decompose.induct) case (3 x y1 y2 ys) have [simp]: "decompose_aux acc xs = map (λx. prod_list x * acc) (prefixes xs)" for acc :: 'a and xs by (induction xs arbitrary: acc) (simp_all add: mult_ac) show ?case using 3[of "fst (bezout_coefficients y1 (y2 * prod_list ys))" "snd (bezout_coefficients y1 (y2 * prod_list ys))"] by (simp add: case_prod_unfold rev_map prefixes_conv_suffixes o_def mult_ac) qed simp_all text ‹ The next function performs the second step: Given a quotient of the form$x / y^n$, it returns a list of$x_0, \ldots, x_n$such that$x / y^n = x_0 / y^n + \ldots + x_{n-1} / y + x_n$and all$x_i$have a Euclidean size less than that of$y$. › fun normalise_decomp :: "('a :: semiring_modulo) ⇒ 'a ⇒ nat ⇒ 'a × ('a list)" where "normalise_decomp x y 0 = (x, [])" | "normalise_decomp x y (Suc n) = ( case normalise_decomp (x div y) y n of (z, rs) ⇒ (z, x mod y # rs))" lemma length_normalise_decomp [simp]: "length (snd (normalise_decomp x y n)) = n" by (induction x y n rule: normalise_decomp.induct) (auto split: prod.split) text ‹ The following constant implements the full process of partial fraction decomposition: The input is a quotient$x / (y_1 ^ {k_1} \ldots y_n ^ {k_n})$and the output is a sum of an entire element and terms of the form$a / y_i ^ k$where$a$has a Euclidean size less than$y_i$. › definition partial_fraction_decomposition :: "'a :: euclidean_ring_gcd ⇒ ('a × nat) list ⇒ 'a × 'a list list" where "partial_fraction_decomposition x ys = (if ys = [] then (x, []) else (let zs = [let (y, n) = ys ! i in normalise_decomp (decompose x (map (λ(y,n). y ^ Suc n) ys) ! i) y (Suc n). i ← [0..<length ys]] in (sum_list (map fst zs), map snd zs)))" lemma length_pfd1 [simp]: "length (snd (partial_fraction_decomposition x ys)) = length ys" by (simp add: partial_fraction_decomposition_def) lemma length_pfd2 [simp]: "i < length ys ⟹ length (snd (partial_fraction_decomposition x ys) ! i) = snd (ys ! i) + 1" by (auto simp: partial_fraction_decomposition_def case_prod_unfold Let_def) lemma size_normalise_decomp: "a ∈ set (snd (normalise_decomp x y n)) ⟹ y ≠ 0 ⟹ euclidean_size a < euclidean_size y" by (induction x y n rule: normalise_decomp.induct) (auto simp: case_prod_unfold Let_def mod_size_less) lemma size_partial_fraction_decomposition: "i < length xs ⟹ fst (xs ! i) ≠ 0 ⟹ x ∈ set (snd (partial_fraction_decomposition y xs) ! i) ⟹ euclidean_size x < euclidean_size (fst (xs ! i))" by (auto simp: partial_fraction_decomposition_def Let_def case_prod_unfold simp del: normalise_decomp.simps split: if_split_asm intro!: size_normalise_decomp) text ‹ A homomorphism$\varphi$from a Euclidean ring$R$into another ring$S$with a notion of division. We will show that for any$x,y\in R$such that$\phi(y)$is a unit, we can perform partial fraction decomposition on the quotient$\varphi(x) / \varphi(y)$. The obvious choice for$S$is the fraction field of$R$, but other choices may also make sense: If, for example,$R$is a ring of polynomials$K[X]$, then one could let$S = K$and$\varphi$the evaluation homomorphism. Or one could let$S = K[[X]]$(the ring of formal power series) and$\varphi$the canonical homomorphism from polynomials to formal power series. › locale pfd_homomorphism = fixes lift :: "('a :: euclidean_ring_gcd) ⇒ ('b :: euclidean_semiring_cancel)" assumes lift_add: "lift (a + b) = lift a + lift b" assumes lift_mult: "lift (a * b) = lift a * lift b" assumes lift_0 [simp]: "lift 0 = 0" assumes lift_1 [simp]: "lift 1 = 1" begin lemma lift_power: "lift (a ^ n) = lift a ^ n" by (induction n) (simp_all add: lift_mult) definition from_decomp :: "'a ⇒ 'a ⇒ nat ⇒ 'b" where "from_decomp x y n = lift x div lift y ^ n" lemma decompose: assumes "ys ≠ []" "pairwise coprime (set ys)" "distinct ys" "⋀y. y ∈ set ys ⟹ is_unit (lift y)" shows "(∑i<length ys. lift (decompose x ys ! i) div lift (ys ! i)) = lift x div lift (prod_list ys)" using assms proof (induction ys arbitrary: x rule: list_nonempty_induct) case (cons y ys x) from cons.prems have "coprime (prod_list ys) y" by (auto simp add: pairwise_insert intro: prod_list_coprime_left) from cons.prems have unit: "is_unit (lift y)" by simp moreover from cons.prems have "∀y∈set ys. is_unit (lift y)" by simp hence unit': "is_unit (lift (prod_list ys))" by (induction ys) (auto simp: lift_mult) ultimately have unit: "lift y dvd b" "lift (prod_list ys) dvd b" for b by auto obtain s t where st: "bezout_coefficients y (prod_list ys) = (s, t)" by (cases "bezout_coefficients y (prod_list ys)") simp_all from ‹pairwise coprime (set (y#ys))› have coprime:"pairwise coprime (set ys)" by (rule pairwise_subset) auto have "(∑i<length (y # ys). lift (decompose x (y # ys) ! i) div lift ((y # ys) ! i)) = lift (t * x) div lift y + lift (s * x) div lift (prod_list ys)" using cons.hyps cons.prems coprime unfolding length_Cons atLeast0LessThan [symmetric] by (subst sum.atLeast_Suc_lessThan, simp, subst sum.shift_bounds_Suc_ivl) (simp add: atLeast0LessThan decompose_rec st cons.IH lift_mult) also have "(lift (t * x) div lift y + lift (s * x) div lift (prod_list ys)) * lift (prod_list (y#ys)) = lift (prod_list ys) * (lift y * (lift (t * x) div lift y)) + lift y * (lift (prod_list ys) * (lift (s * x) div lift (prod_list ys)))" by (simp_all add: lift_mult algebra_simps) also have "… = lift (prod_list ys * t * x + y * s * x)" using assms unit by (simp add: lift_mult lift_add algebra_simps) finally have "(∑i<length (y # ys). lift (decompose x (y # ys) ! i) div lift ((y # ys) ! i)) = lift ((s * y + t * prod_list ys) * x) div lift (prod_list (y#ys))" using unit by (subst unit_eq_div2) (auto simp: lift_mult lift_add algebra_simps) also have "s * y + t * prod_list ys = gcd (prod_list ys) y" using bezout_coefficients_fst_snd[of y "prod_list ys"] by (simp add: st gcd.commute) also have "… = 1" using ‹coprime (prod_list ys) y› by simp finally show ?case by simp qed simp_all lemma normalise_decomp: fixes x y :: 'a and n :: nat assumes "is_unit (lift y)" defines "xs ≡ snd (normalise_decomp x y n)" shows "lift (fst (normalise_decomp x y n)) + (∑i<n. from_decomp (xs!i) y (n-i)) = lift x div lift y ^ n" using assms unfolding xs_def proof (induction x y n rule: normalise_decomp.induct, goal_cases) case (2 x y n) from 2(2) have unit: "is_unit (lift y ^ n)" by (simp add: is_unit_power_iff) obtain a b where ab: "normalise_decomp (x div y) y n = (a, b)" by (cases "normalise_decomp (x div y) y n") simp_all have "lift (fst (normalise_decomp x y (Suc n))) + (∑i<Suc n. from_decomp (snd (normalise_decomp x y (Suc n)) ! i) y (Suc n - i)) = lift a + (∑i<n. from_decomp (b ! i) y (n - i)) + from_decomp (x mod y) y (Suc n)" unfolding atLeast0LessThan[symmetric] apply (subst sum.atLeast_Suc_lessThan) apply simp apply (subst sum.shift_bounds_Suc_ivl) apply (simp add: ab atLeast0LessThan ac_simps) done also have "lift a + (∑i<n. from_decomp (b ! i) y (n - i)) = lift (x div y) div lift y ^ n" using 2 by (simp add: ab) also from 2(2) unit have "(… + from_decomp (x mod y) y (Suc n)) * lift y = (lift ((x div y) * y + x mod y) div lift y ^ n)" (is "?A * _ = ?B div _") unfolding lift_add lift_mult apply (subst div_add) apply (auto simp add: from_decomp_def algebra_simps dvd_div_mult2_eq unit_div_mult_swap dvd_div_mult2_eq[OF unit_imp_dvd] is_unit_mult_iff) done with 2(2) have "?A = … div lift y" by (subst eq_commute, subst dvd_div_eq_mult) auto also from 2(2) unit have "… = ?B div (lift y ^ Suc n)" by (subst is_unit_div_mult2_eq [symmetric]) (auto simp: mult_ac) also have "x div y * y + x mod y = x" by (rule div_mult_mod_eq) finally show ?case . qed simp_all lemma lift_prod_list: "lift (prod_list xs) = prod_list (map lift xs)" by (induction xs) (simp_all add: lift_mult) lemma lift_sum: "lift (sum f A) = sum (λx. lift (f x)) A" by (cases "finite A", induction A rule: finite_induct) (simp_all add: lift_add) lemma partial_fraction_decomposition: fixes ys :: "('a × nat) list" defines "ys' ≡ map (λ(x,n). x ^ Suc n) ys :: 'a list" assumes unit: "⋀y. y ∈ fst  set ys ⟹ is_unit (lift y)" assumes coprime: "pairwise coprime (set ys')" assumes distinct: "distinct ys'" assumes "partial_fraction_decomposition x ys = (a, zs)" shows "lift a + (∑i<length ys. ∑j≤snd (ys!i). from_decomp (zs!i!j) (fst (ys!i)) (snd (ys!i)+1 - j)) = lift x div lift (prod_list ys')" proof (cases "ys = []") assume [simp]: "ys ≠ []" define n where "n = length ys" have "lift x div lift (prod_list ys') = (∑i<n. lift (decompose x ys' ! i) div lift (ys' ! i))" using assms by (subst decompose [symmetric]) (force simp: lift_prod_list prod_list_zero_iff lift_power lift_mult o_def n_def is_unit_mult_iff is_unit_power_iff)+ also have "… = (∑i<n. lift (fst (normalise_decomp (decompose x ys' ! i) (fst (ys!i)) (snd (ys!i)+1)))) + (∑i<n. (∑j≤snd (ys!i). from_decomp (zs!i!j) (fst (ys!i)) (snd (ys!i)+1 - j)))" (is "_ = ?A + ?B") proof (subst sum.distrib [symmetric], intro sum.cong refl, goal_cases) case (1 i) from 1 have "lift (ys' ! i) = lift (fst (ys ! i)) ^ Suc (snd (ys ! i))" by (simp add: ys'_def n_def lift_power lift_mult split: prod.split) also from 1 have "lift (decompose x ys' ! i) div … = lift (fst (normalise_decomp (decompose x ys' ! i) (fst (ys!i)) (snd (ys!i)+1))) + (∑j<Suc (snd (ys ! i)). from_decomp (snd (normalise_decomp (decompose x ys' ! i) (fst (ys!i)) (snd (ys!i)+1)) ! j) (fst (ys ! i)) (snd (ys!i)+1 - j))" (is "_ = _ + ?C") by (subst normalise_decomp [symmetric]) (simp_all add: n_def unit) also have "?C = (∑j≤snd (ys!i). from_decomp (zs!i!j) (fst (ys!i)) (snd (ys!i)+1 - j))" using assms 1 by (intro sum.cong refl) (auto simp: partial_fraction_decomposition_def case_prod_unfold Let_def o_def n_def simp del: normalise_decomp.simps) finally show ?case . qed also from assms have "?A = lift a" by (auto simp: partial_fraction_decomposition_def o_def sum_list_sum_nth atLeast0LessThan case_prod_unfold Let_def lift_sum n_def intro!: sum.cong) finally show ?thesis by (simp add: n_def) qed (insert assms, simp add: partial_fraction_decomposition_def) end subsection ‹Specific results for polynomials› definition divmod_field_poly :: "'a :: field poly ⇒ 'a poly ⇒ 'a poly × 'a poly" where "divmod_field_poly p q = (p div q, p mod q)" lemma divmod_field_poly_code [code]: "divmod_field_poly p q = (let cg = coeffs q in if cg = [] then (0, p) else let cf = coeffs p; ilc = inverse (last cg); ch = map ((*) ilc) cg; (q, r) = divmod_poly_one_main_list [] (rev cf) (rev ch) (1 + length cf - length cg) in (poly_of_list (map ((*) ilc) q), poly_of_list (rev r)))" unfolding divmod_field_poly_def by (rule pdivmod_via_divmod_list) definition normalise_decomp_poly :: "'a::field_gcd poly ⇒ 'a poly ⇒ nat ⇒ 'a poly × 'a poly list" where [simp]: "normalise_decomp_poly (p :: _ poly) q n = normalise_decomp p q n" lemma normalise_decomp_poly_code [code]: "normalise_decomp_poly x y 0 = (x, [])" "normalise_decomp_poly x y (Suc n) = ( let (x', r) = divmod_field_poly x y; (z, rs) = normalise_decomp_poly x' y n in (z, r # rs))" by (simp_all add: divmod_field_poly_def) definition poly_pfd_simple where "poly_pfd_simple x cs = (if cs = [] then (x, []) else (let zs = [let (c, n) = cs ! i in normalise_decomp_poly (decompose x (map (λ(c,n). [:1,-c:] ^ Suc n) cs) ! i) [:1,-c:] (n+1). i ← [0..<length cs]] in (sum_list (map fst zs), map (map (λp. coeff p 0) ∘ snd) zs)))" lemma poly_pfd_simple_code [code]: "poly_pfd_simple x cs = (if cs = [] then (x, []) else let zs = zip_with (λ(c,n) decomp. normalise_decomp_poly decomp [:1,-c:] (n+1)) cs (decompose x (map (λ(c,n). [:1,-c:] ^ Suc n) cs)) in (sum_list (map fst zs), map (map (λp. coeff p 0) ∘ snd) zs))" unfolding poly_pfd_simple_def zip_with_altdef' by (simp add: Let_def case_prod_unfold) lemma fst_poly_pfd_simple: "fst (poly_pfd_simple x cs) = fst (partial_fraction_decomposition x (map (λ(c,n). ([:1,-c:],n)) cs))" by (auto simp: poly_pfd_simple_def partial_fraction_decomposition_def o_def case_prod_unfold Let_def sum_list_sum_nth intro!: sum.cong) lemma const_polyI: "degree p = 0 ⟹ [:coeff p 0:] = p" by (elim degree_eq_zeroE) simp_all lemma snd_poly_pfd_simple: "map (map (λc. [:c :: 'a :: field_gcd:])) (snd (poly_pfd_simple x cs)) = (snd (partial_fraction_decomposition x (map (λ(c,n). ([:1,-c:],n)) cs)))" proof - have "snd (poly_pfd_simple x cs) = map (map (λp. coeff p 0)) (snd (partial_fraction_decomposition x (map (λ(c,n). ([:1,-c:],n)) cs)))" (is "_ = map ?f ?B") by (auto simp: poly_pfd_simple_def partial_fraction_decomposition_def o_def case_prod_unfold Let_def sum_list_sum_nth intro!: sum.cong) also have "map (map (λc. [:c:])) (map ?f ?B) = map (map (λx. x)) ?B" unfolding map_map o_def proof (intro map_cong refl const_polyI, goal_cases) case (1 ys y) from 1 obtain i where i: "i < length cs" "ys = snd (partial_fraction_decomposition x (map (λ(c,n). ([:1,-c:],n)) cs)) ! i" by (auto simp: in_set_conv_nth) with 1 have "euclidean_size y < euclidean_size (fst (map (λ(c,n). ([:1,-c:],n)) cs ! i))" by (intro size_partial_fraction_decomposition[of i _ _ x]) (auto simp: case_prod_unfold Let_def) with i(1) have "euclidean_size y < 2" by (auto simp: case_prod_unfold Let_def euclidean_size_poly_def split: if_split_asm) thus ?case by (cases y rule: pCons_cases) (auto simp: euclidean_size_poly_def split: if_split_asm) qed finally show ?thesis by simp qed lemma poly_pfd_simple: "partial_fraction_decomposition x (map (λ(c,n). ([:1,-c:],n)) cs) = (fst (poly_pfd_simple x cs), map (map (λc. [:c:])) (snd (poly_pfd_simple x cs)))" by (simp add: fst_poly_pfd_simple snd_poly_pfd_simple) end  # Theory Factorizations (* File: Factorizations.thy Author: Manuel Eberl, TU München *) section ‹Factorizations of polynomials› theory Factorizations imports Complex_Main Linear_Recurrences_Misc "HOL-Computational_Algebra.Computational_Algebra" "HOL-Computational_Algebra.Polynomial_Factorial" begin text ‹ We view a factorisation of a polynomial as a pair consisting of the leading coefficient and a list of roots with multiplicities. This gives us a factorization into factors of the form$(X - c) ^ {n+1}$. › definition interp_factorization where "interp_factorization = (λ(a,cs). Polynomial.smult a (∏(c,n)←cs. [:-c,1:] ^ Suc n))" text ‹ An alternative way to factorise is as a pair of the leading coefficient and factors of the form$(1 - cX) ^ {n+1}$. › definition interp_alt_factorization where "interp_alt_factorization = (λ(a,cs). Polynomial.smult a (∏(c,n)←cs. [:1,-c:] ^ Suc n))" definition is_factorization_of where "is_factorization_of fctrs p = (interp_factorization fctrs = p ∧ distinct (map fst (snd fctrs)))" definition is_alt_factorization_of where "is_alt_factorization_of fctrs p = (interp_alt_factorization fctrs = p ∧ 0 ∉ set (map fst (snd fctrs)) ∧ distinct (map fst (snd fctrs)))" text ‹ Regular and alternative factorisations are related by reflecting the polynomial. › lemma interp_factorization_reflect: assumes "(0::'a::idom) ∉ fst  set (snd fctrs)" shows "reflect_poly (interp_factorization fctrs) = interp_alt_factorization fctrs" proof - have "reflect_poly (interp_factorization fctrs) = Polynomial.smult (fst fctrs) (∏x←snd fctrs. reflect_poly [:- fst x, 1:] ^ Suc (snd x))" by (simp add: interp_factorization_def interp_alt_factorization_def case_prod_unfold reflect_poly_smult reflect_poly_prod_list reflect_poly_power o_def del: power_Suc) also have "map (λx. reflect_poly [:- fst x, 1:] ^ Suc (snd x)) (snd fctrs) = map (λx. [:1, - fst x:] ^ Suc (snd x)) (snd fctrs)" using assms by (intro list.map_cong0, subst reflect_poly_pCons) auto also have "Polynomial.smult (fst fctrs) (prod_list …) = interp_alt_factorization fctrs" by (simp add: interp_alt_factorization_def case_prod_unfold) finally show ?thesis . qed lemma interp_alt_factorization_reflect: assumes "(0::'a::idom) ∉ fst  set (snd fctrs)" shows "reflect_poly (interp_alt_factorization fctrs) = interp_factorization fctrs" proof - have "reflect_poly (interp_alt_factorization fctrs) = Polynomial.smult (fst fctrs) (∏x←snd fctrs. reflect_poly [:1, - fst x:] ^ Suc (snd x))" by (simp add: interp_factorization_def interp_alt_factorization_def case_prod_unfold reflect_poly_smult reflect_poly_prod_list reflect_poly_power o_def del: power_Suc) also have "map (λx. reflect_poly [:1, - fst x:] ^ Suc (snd x)) (snd fctrs) = map (λx. [:- fst x, 1:] ^ Suc (snd x)) (snd fctrs)" proof (intro list.map_cong0, clarsimp simp del: power_Suc, goal_cases) fix c n assume "(c, n) ∈ set (snd fctrs)" with assms have "c ≠ 0" by force thus "reflect_poly [:1, -c:] ^ Suc n = [:-c, 1:] ^ Suc n" by (simp add: reflect_poly_pCons del: power_Suc) qed also have "Polynomial.smult (fst fctrs) (prod_list …) = interp_factorization fctrs" by (simp add: interp_factorization_def case_prod_unfold) finally show ?thesis . qed lemma coeff_0_interp_factorization: "coeff (interp_factorization fctrs) 0 = (0 :: 'a :: idom) ⟷ fst fctrs = 0 ∨ 0 ∈ fst  set (snd fctrs)" by (force simp: interp_factorization_def case_prod_unfold coeff_0_prod_list o_def coeff_0_power prod_list_zero_iff simp del: power_Suc) lemma reflect_factorization: assumes "coeff p 0 ≠ (0::'a::idom)" assumes "is_factorization_of fctrs p" shows "is_alt_factorization_of fctrs (reflect_poly p)" using assms by (force simp: interp_factorization_reflect is_factorization_of_def is_alt_factorization_of_def coeff_0_interp_factorization) lemma reflect_factorization': assumes "coeff p 0 ≠ (0::'a::idom)" assumes "is_alt_factorization_of fctrs p" shows "is_factorization_of fctrs (reflect_poly p)" using assms by (force simp: interp_alt_factorization_reflect is_factorization_of_def is_alt_factorization_of_def coeff_0_interp_factorization) lemma zero_in_factorization_iff: assumes "is_factorization_of fctrs p" shows "coeff p 0 = 0 ⟷ p = 0 ∨ (0::'a::idom) ∈ fst  set (snd fctrs)" proof (cases "p = 0") assume "p ≠ 0" with assms have [simp]: "fst fctrs ≠ 0" by (auto simp: is_factorization_of_def interp_factorization_def case_prod_unfold) from assms have "p = interp_factorization fctrs" by (simp add: is_factorization_of_def) also have "coeff … 0 = 0 ⟷ 0 ∈ fst  set (snd fctrs)" by (force simp add: interp_factorization_def case_prod_unfold coeff_0_prod_list prod_list_zero_iff o_def coeff_0_power) finally show ?thesis using ‹p ≠ 0› by blast next assume p: "p = 0" with assms have "interp_factorization fctrs = 0" by (simp add: is_factorization_of_def) also have "interp_factorization fctrs = 0 ⟷ fst fctrs = 0 ∨ (∏(c,n)←snd fctrs. [:-c,1:]^Suc n) = 0" by (simp add: interp_factorization_def case_prod_unfold) also have "(∏(c,n)←snd fctrs. [:-c,1:]^Suc n) = 0 ⟷ False" by (auto simp: prod_list_zero_iff simp del: power_Suc) finally show ?thesis by (simp add: ‹p = 0›) qed lemma poly_prod_list [simp]: "poly (prod_list ps) x = prod_list (map (λp. poly p x) ps)" by (induction ps) auto lemma is_factorization_of_roots: fixes a :: "'a :: idom" assumes "is_factorization_of (a, fctrs) p" "p ≠ 0" shows "set (map fst fctrs) = {x. poly p x = 0}" using assms by (force simp: is_factorization_of_def interp_factorization_def o_def case_prod_unfold prod_list_zero_iff simp del: power_Suc) lemma (in monoid_mult) prod_list_prod_nth: "prod_list xs = (∏i<length xs. xs ! i)" by (induction xs) (auto simp: prod.lessThan_Suc_shift simp del: prod.lessThan_Suc) lemma order_prod: assumes "⋀x. x ∈ A ⟹ f x ≠ 0" assumes "⋀x y. x ∈ A ⟹ y ∈ A ⟹ x ≠ y ⟹ coprime (f x) (f y)" shows "order c (prod f A) = (∑x∈A. order c (f x))" using assms proof (induction A rule: infinite_finite_induct) case (insert x A) from insert.hyps have "order c (prod f (insert x A)) = order c (f x * prod f A)" by simp also have "… = order c (f x) + order c (prod f A)" using insert.prems and insert.hyps by (intro order_mult) auto also have "order c (prod f A) = (∑x∈A. order c (f x))" using insert.prems and insert.hyps by (intro insert.IH) auto finally show ?case using insert.hyps by simp qed auto lemma is_factorization_of_order: fixes p :: "'a :: field_gcd poly" assumes "p ≠ 0" assumes "is_factorization_of (a, fctrs) p" assumes "(c, n) ∈ set fctrs" shows "order c p = Suc n" proof - from assms have distinct: "distinct (map fst (fctrs))" by (simp add: is_factorization_of_def) from assms have [simp]: "a ≠ 0" by (auto simp: is_factorization_of_def interp_factorization_def) from assms(2) have "p = interp_factorization (a, fctrs)" unfolding is_factorization_of_def by simp also have "order c … = order c (∏(c,n)←fctrs. [:-c, 1:] ^ Suc n)" unfolding interp_factorization_def by (simp add: order_smult) also have "(∏(c,n)←fctrs. [:-c, 1:] ^ Suc n) = (∏i∈{..<length fctrs}. [:-fst (fctrs ! i), 1:] ^ Suc (snd (fctrs ! i)))" by (simp add: prod_list_prod_nth case_prod_unfold) also have "order c … = (∑x<length fctrs. order c ([:- fst (fctrs ! x), 1:] ^ Suc (snd (fctrs ! x))))" proof (rule order_prod) fix i assume "i ∈ {..<length fctrs}" then show "[:- fst (fctrs ! i), 1:] ^ Suc (snd (fctrs ! i)) ≠ 0" by (simp only: power_eq_0_iff) simp next fix i j :: nat assume "i ≠ j" "i ∈ {..<length fctrs}" "j ∈ {..<length fctrs}" then have "fst (fctrs ! i) ≠ fst (fctrs ! j)" using nth_eq_iff_index_eq [OF distinct, of i j] by simp then show "coprime ([:- fst (fctrs ! i), 1:] ^ Suc (snd (fctrs ! i))) ([:- fst (fctrs ! j), 1:] ^ Suc (snd (fctrs ! j)))" by (simp only: coprime_power_left_iff coprime_power_right_iff) (auto simp add: coprime_linear_poly) qed also have "… = (∑(c',n')←fctrs. order c ([:-c', 1:] ^ Suc n'))" by (simp add: sum_list_sum_nth case_prod_unfold atLeast0LessThan) also have "… = (∑(c',n')←fctrs. if c = c' then Suc n' else 0)" by (intro arg_cong[OF map_cong]) (auto simp add: order_power_n_n order_0I simp del: power_Suc) also have "… = (∑x←fctrs. if x = (c, n) then Suc (snd x) else 0)" using distinct assms by (intro arg_cong[OF map_cong]) (force simp: distinct_map inj_on_def)+ also from distinct have "… = (∑x∈set fctrs. if x = (c, n) then Suc (snd x) else 0)" by (intro sum_list_distinct_conv_sum_set) (simp_all add: distinct_map) also from assms have "… = Suc n" by simp finally show ?thesis . qed text ‹ For complex polynomials, a factorisation in the above sense always exists. › lemma complex_factorization_exists: "∃fctrs. is_factorization_of fctrs (p :: complex poly)" proof (cases "p = 0") case True thus ?thesis by (intro exI[of _ "(0, [])"]) (auto simp: is_factorization_of_def interp_factorization_def) next case False hence "∃xs. set xs = {x. poly p x = 0} ∧ distinct xs" by (intro finite_distinct_list poly_roots_finite) then obtain xs where [simp]: "set xs = {x. poly p x = 0}" "distinct xs" by blast have "interp_factorization (lead_coeff p, map (λx. (x, order x p - 1)) xs) = smult (lead_coeff p) (∏x←xs. [:- x, 1:] ^ Suc (order x p - 1))" by (simp add: interp_factorization_def o_def) also have "(∏x←xs. [:- x, 1:] ^ Suc (order x p - 1)) = (∏x|poly p x = 0. [:- x, 1:] ^ Suc (order x p - 1))" by (subst prod.distinct_set_conv_list [symmetric]) simp_all also have "… = (∏x|poly p x = 0. [:- x, 1:] ^ order x p)" proof (intro prod.cong refl, goal_cases) case (1 x) with False have "order x p ≠ 0" by (subst (asm) order_root) auto hence *: "Suc (order x p - 1) = order x p" by simp show ?case by (simp only: *) qed also have "smult (lead_coeff p) … = p" by (rule complex_poly_decompose) finally have "is_factorization_of (lead_coeff p, map (λx. (x, order x p - 1)) xs) p" by (auto simp: is_factorization_of_def o_def) thus ?thesis .. qed text ‹ By reflecting the polynomial, this means that for complex polynomials with non-zero constant coefficient, the alternative factorisation also exists. › corollary complex_alt_factorization_exists: assumes "coeff p 0 ≠ 0" shows "∃fctrs. is_alt_factorization_of fctrs (p :: complex poly)" proof - from assms have "coeff (reflect_poly p) 0 ≠ 0" by auto moreover from complex_factorization_exists [of "reflect_poly p"] obtain fctrs where "is_factorization_of fctrs (reflect_poly p)" .. ultimately have "is_alt_factorization_of fctrs (reflect_poly (reflect_poly p))" by (rule reflect_factorization) also from assms have "reflect_poly (reflect_poly p) = p" by simp finally show ?thesis .. qed end # Theory Rational_FPS_Solver (* File: Rational_FPS_Solver.thy Author: Manuel Eberl, TU München *) section ‹Solver for rational formal power series› theory Rational_FPS_Solver imports Complex_Main Pochhammer_Polynomials Partial_Fraction_Decomposition Factorizations "HOL-Computational_Algebra.Field_as_Ring" begin text ‹ We can determine the$k$-th coefficient of an FPS of the form$d/(1-cX)^n$, which is an important step in solving linear recurrences. The$k$-th coefficient of such an FPS is always of the form$p(k) c^k$where$p$is the following polynomial: › definition inverse_irred_power_poly :: "'a :: field_char_0 ⇒ nat ⇒ 'a poly" where "inverse_irred_power_poly d n = Poly [(d * of_nat (stirling n (k+1))) / (fact (n - 1)). k ← [0..<n]]" lemma one_minus_const_fps_X_neg_power'': fixes c :: "'a :: field_char_0" assumes n: "n > 0" shows "fps_const d / ((1 - fps_const (c :: 'a :: field_char_0) * fps_X) ^ n) = Abs_fps (λk. poly (inverse_irred_power_poly d n) (of_nat k) * c^k)" (is "?lhs = ?rhs") proof (rule fps_ext) include fps_notation fix k :: nat let ?p = "smult (d / (fact (n - 1))) (pcompose (pochhammer_poly (n - 1)) [:1,1:])" from n have "?lhs = fps_const d * inverse ((1 - fps_const c * fps_X) ^ n)" by (subst fps_divide_unit) auto also have "inverse ((1 - fps_const c * fps_X) ^ n) = Abs_fps (λk. of_nat ((n + k - 1) choose k) * c^k)" by (intro one_minus_const_fps_X_neg_power' n) also have "(fps_const d * …)$ k  = d * of_nat ((n + k - 1) choose k) * c^k" by simp
also from n have "(n + k - 1 choose k) = (n + k - 1 choose (n - 1))"
by (subst binomial_symmetric) simp_all
also from n have "of_nat … = (pochhammer (of_nat k + 1) (n - 1) / fact (n - 1) :: 'a)"
by (simp_all add: binomial_gbinomial gbinomial_pochhammer' of_nat_diff)
also have "d * … = poly ?p (of_nat k)"
by (simp add: divide_inverse eval_pochhammer_poly poly_pcompose add_ac)
also {
from assms have "pCons 0 (pcompose (pochhammer_poly (n-1)) [:1,1::'a:]) = pochhammer_poly n"
by (subst pochhammer_poly_Suc' [symmetric]) simp
also from assms have "… = pCons 0 (Poly [of_nat (stirling n (k+1)). k ← [0..<Suc n]])"
unfolding pochhammer_poly_def
by (auto simp add: poly_eq_iff nth_default_def coeff_pCons
split: nat.split simp del: upt_Suc )
finally have "pcompose (pochhammer_poly (n-1)) [:1,1::'a:] =
Poly [of_nat (stirling n (k+1)). k ← [0..<Suc n]]" by simp
}
also have "smult (d / fact (n - 1)) (Poly [of_nat (stirling n (k+1)). k ← [0..<Suc n]]) =
inverse_irred_power_poly d n"
by (auto simp: poly_eq_iff inverse_irred_power_poly_def nth_default_def)
also have "poly … (of_nat k) * c ^ k = ?rhs $k" by simp finally show "?lhs$ k = ?rhs $k" . qed lemma inverse_irred_power_poly_code [code abstract]: "coeffs (inverse_irred_power_poly d n) = (if n = 0 ∨ d = 0 then [] else let e = d / (fact (n - 1)) in [e * of_nat x. x ← tl (stirling_row n)])" proof (cases "n = 0 ∨ d = 0") case False define e where "e = d / (fact (n - 1))" from False have "coeffs (inverse_irred_power_poly d n) = [e * of_nat (stirling n (k+1)). k ← [0..<n]]" by (auto simp: inverse_irred_power_poly_def Let_def divide_inverse mult_ac last_map stirling_row_def map_tl [symmetric] tl_upt e_def no_trailing_unfold) also have "… = [e * of_nat x. x ← tl (stirling_row n)]" by (simp add: stirling_row_def map_tl [symmetric] o_def tl_upt map_Suc_upt [symmetric] del: upt_Suc) finally show ?thesis using False by (simp add: Let_def e_def) qed (auto simp: inverse_irred_power_poly_def) lemma solve_rat_fps_aux: fixes p :: "'a :: {field_char_0,field_gcd} poly" and cs :: "('a × nat) list" assumes distinct: "distinct (map fst cs)" assumes azs: "(a, zs) = poly_pfd_simple p cs" assumes nz: "0 ∉ fst  set cs" shows "fps_of_poly p / fps_of_poly (∏(c,n)←cs. [:1,-c:]^Suc n) = Abs_fps (λk. coeff a k + (∑i<length cs. poly (∑j≤snd (cs ! i). (inverse_irred_power_poly (zs ! i ! j) (snd (cs ! i)+1 - j))) (of_nat k) * (fst (cs ! i)) ^ k))" (is "_ = ?rhs") proof - interpret pfd_homomorphism "fps_of_poly :: 'a poly ⇒ 'a fps" by standard (auto simp: fps_of_poly_add fps_of_poly_mult) from distinct have distinct': "(a, b1) ∈ set cs ⟹ (a, b2) ∈ set cs ⟹ b1 = b2" for a b1 b2 by (metis (no_types, hide_lams) Some_eq_map_of_iff image_set in_set_zipE insert_iff list.simps(15) map_of_Cons_code(2) map_of_SomeD nz snd_conv) from nz have nz': "(0, b) ∉ set cs" for b by (auto simp add: image_iff) define n where "n = length cs" let ?g = "λ(c, n). [:1, - c:] ^ Suc n" have "inj_on ?g (set cs)" proof fix x y assume "x ∈ set cs" "y ∈ set cs" "?g x = ?g y" moreover obtain c1 n1 c2 n2 where [simp]: "x = (c1, n1)" "y = (c2, n2)" by (cases x, cases y) ultimately have in_cs: "(c1, n1) ∈ set cs" "(c2, n2) ∈ set cs" and eq: "[:1, - c1:] ^ Suc n1 = [:1, - c2:] ^ Suc n2" by simp_all with nz have [simp]: "c1 ≠ 0" "c2 ≠ 0" by (auto simp add: image_iff) have "Suc n1 = degree ([:1, - c1:] ^ Suc n1)" by (simp add: degree_power_eq del: power_Suc) also have "… = degree ([:1, - c2:] ^ Suc n2)" using eq by simp also have "… = Suc n2" by (simp add: degree_power_eq del: power_Suc) finally have "n1 = n2" by simp then have "0 = poly ([:1, - c1:] ^ Suc n1) (1 / c1)" by simp also have "… = poly ([:1, - c2:] ^ Suc n2) (1 / c1)" using eq by simp finally show "x = y" using ‹n1 = n2› by (auto simp: field_simps) qed with distinct have distinct': "distinct (map ?g cs)" by (simp add: distinct_map del: power_Suc) from nz' distinct have coprime: "pairwise coprime (?g  set cs)" by (auto intro!: pairwise_imageI coprime_linear_poly' simp add: eq_key_imp_eq_value simp del: power_Suc) have [simp]: "length zs = n" using assms by (simp add: poly_pfd_simple_def n_def split: if_split_asm) have [simp]: "i < length cs ⟹ length (zs!i) = snd (cs!i)+1" for i using assms by (simp add: poly_pfd_simple_def Let_def case_prod_unfold split: if_split_asm) let ?f = "λ(c, n). ([:1,-c:], n)" let ?cs' = "map ?f cs" have "fps_of_poly (fst (poly_pfd_simple p cs)) + (∑i<length ?cs'. ∑j≤snd (?cs' ! i). from_decomp (map (map (λc. [:c:])) (snd (poly_pfd_simple p cs)) ! i ! j) (fst (?cs' ! i)) (snd (?cs' ! i)+1 - j)) = fps_of_poly p / fps_of_poly (∏(x, n)←?cs'. x ^ Suc n)" (is "?A = ?B") using nz distinct' coprime by (intro partial_fraction_decomposition poly_pfd_simple) (force simp: o_def case_prod_unfold simp del: power_Suc)+ note this [symmetric] also from azs [symmetric] have "?A = fps_of_poly a + (∑i<n. ∑j≤snd (cs ! i). from_decomp (map (map (λc. [:c:])) zs ! i ! j) [:1,-fst (cs ! i):] (snd (cs ! i)+1 - j))" (is "_ = _ + ?S") by (simp add: case_prod_unfold Let_def n_def) also have "?S = (∑i<length cs. ∑j≤snd (cs ! i). fps_const (zs ! i ! j) / ((1 - fps_const (fst (cs!i))*fps_X) ^ (snd (cs!i)+1 - j)))" by (intro sum.cong refl) (auto simp: from_decomp_def map_nth n_def fps_of_poly_linear' fps_of_poly_simps fps_const_neg [symmetric] mult_ac simp del: fps_const_neg) also have "… = (∑i<length cs. ∑j≤snd (cs ! i) . Abs_fps (λk. poly (inverse_irred_power_poly (zs ! i ! j) (snd (cs ! i)+1 - j)) (of_nat k) * (fst (cs ! i)) ^ k))" using nz by (intro sum.cong refl one_minus_const_fps_X_neg_power'') auto also have "fps_of_poly a + … = ?rhs" by (intro fps_ext) (simp_all add: sum_distrib_right fps_sum_nth poly_sum) finally show ?thesis by (simp add: o_def case_prod_unfold) qed definition solve_factored_ratfps :: "('a :: {field_char_0,field_gcd}) poly ⇒ ('a × nat) list ⇒ 'a poly × ('a poly × 'a) list" where "solve_factored_ratfps p cs = (let n = length cs in case poly_pfd_simple p cs of (a, zs) ⇒ (a, zip_with (λzs (c,n). ((∑(z,j) ← zip zs [0..<Suc n]. inverse_irred_power_poly z (n + 1 - j)), c)) zs cs))" lemma length_snd_poly_pfd_simple [simp]: "length (snd (poly_pfd_simple p cs)) = length cs" by (simp add: poly_pfd_simple_def) lemma length_nth_snd_poly_pfd_simple [simp]: "i < length cs ⟹ length (snd (poly_pfd_simple p cs) ! i) = snd (cs!i) + 1" by (auto simp: poly_pfd_simple_def case_prod_unfold Let_def) lemma solve_factored_ratfps_roots: "map snd (snd (solve_factored_ratfps p cs)) = map fst cs" by (rule nth_equalityI) (simp_all add: solve_factored_ratfps_def poly_pfd_simple case_prod_unfold Let_def zip_with_altdef o_def) definition interp_ratfps_solution where "interp_ratfps_solution = (λ(p,cs) n. coeff p n + (∑(q,c)←cs. poly q (of_nat n) * c ^ n))" lemma solve_factored_ratfps: fixes p :: "'a :: {field_char_0,field_gcd} poly" and cs :: "('a × nat) list" assumes distinct: "distinct (map fst cs)" assumes nz: "0 ∉ fst  set cs" shows "fps_of_poly p / fps_of_poly (∏(c,n)←cs. [:1,-c:]^Suc n) = Abs_fps (interp_ratfps_solution (solve_factored_ratfps p cs))" (is "?lhs = ?rhs") proof - obtain a zs where azs: "(a, zs) = solve_factored_ratfps p cs" using prod.exhaust by metis from azs have a: "a = fst (poly_pfd_simple p cs)" by (simp add: solve_factored_ratfps_def Let_def case_prod_unfold) define zs' where "zs' = snd (poly_pfd_simple p cs)" with a have azs': "(a, zs') = poly_pfd_simple p cs" by simp from azs have zs: "zs = snd (solve_factored_ratfps p cs)" by (auto simp add: snd_def split: prod.split) have "?lhs = Abs_fps (λk. coeff a k + (∑i<length cs. poly (∑j≤snd (cs ! i). inverse_irred_power_poly (zs' ! i ! j) (snd (cs ! i)+1 - j)) (of_nat k) * (fst (cs ! i)) ^ k))" by (rule solve_rat_fps_aux[OF distinct azs' nz]) also from azs have "… = ?rhs" unfolding interp_ratfps_solution_def by (auto simp: a zs solve_factored_ratfps_def Let_def case_prod_unfold zip_altdef zip_with_altdef' sum_list_sum_nth atLeast0LessThan zs'_def lessThan_Suc_atMost intro!: fps_ext sum.cong simp del: upt_Suc) finally show ?thesis . qed definition solve_factored_ratfps' where "solve_factored_ratfps' = (λp (a,cs). solve_factored_ratfps (smult (inverse a) p) cs)" lemma solve_factored_ratfps': assumes "is_alt_factorization_of fctrs q" "q ≠ 0" shows "Abs_fps (interp_ratfps_solution (solve_factored_ratfps' p fctrs)) = fps_of_poly p / fps_of_poly q" proof - from assms have q: "q = interp_alt_factorization fctrs" by (simp add: is_alt_factorization_of_def) from assms(2) have nz: "fst fctrs ≠ 0" by (subst (asm) q) (auto simp: interp_alt_factorization_def case_prod_unfold) note q also from nz have "coeff (interp_alt_factorization fctrs) 0 ≠ 0" by (auto simp: interp_alt_factorization_def case_prod_unfold coeff_0_prod_list o_def coeff_0_power prod_list_zero_iff) finally have "coeff q 0 ≠ 0" . obtain a cs where fctrs: "fctrs = (a, cs)" by (cases fctrs) simp_all obtain b zs where sol: "solve_factored_ratfps' p fctrs = (b, zs)" using prod.exhaust by metis from assms have [simp]: "a ≠ 0" by (auto simp: is_alt_factorization_of_def interp_alt_factorization_def fctrs) have "fps_of_poly p / fps_of_poly (smult a (∏(c, n)←cs. [:1, - c:] ^ Suc n)) = fps_of_poly p / (fps_const a * fps_of_poly (∏(c, n)←cs. [:1, - c:] ^ Suc n))" by (simp_all add: fps_of_poly_smult case_prod_unfold del: power_Suc) also have "… = fps_of_poly p / fps_const a / fps_of_poly (∏(c, n)←cs. [:1, - c:] ^ Suc n)" by (subst is_unit_div_mult2_eq) (auto simp: coeff_0_power coeff_0_prod_list prod_list_zero_iff) also have "fps_of_poly p / fps_const a = fps_of_poly (smult (inverse a) p)" by (simp add: fps_const_inverse fps_divide_unit) also from assms have "smult a (∏(c, n)←cs. [:1, - c:] ^ Suc n) = q" by (simp add: is_alt_factorization_of_def interp_alt_factorization_def fctrs del: power_Suc) also have "fps_of_poly (smult (inverse a) p) / fps_of_poly (∏(c, n)←cs. [:1, - c:] ^ Suc n) = Abs_fps (interp_ratfps_solution (solve_factored_ratfps (smult (inverse a) p) cs))" (is "?lhs = _") using assms by (intro solve_factored_ratfps) (simp_all add: is_alt_factorization_of_def fctrs solve_factored_ratfps'_def) also have "… = Abs_fps (interp_ratfps_solution (solve_factored_ratfps' p fctrs))" by (simp add: solve_factored_ratfps'_def fctrs) finally show ?thesis .. qed lemma degree_Poly_eq: assumes "xs = [] ∨ last xs ≠ 0" shows "degree (Poly xs) = length xs - 1" proof - from assms consider "xs = []" | "xs ≠ []" "last xs ≠ 0" by blast thus ?thesis proof cases assume "last xs ≠ 0" "xs ≠ []" hence "no_trailing ((=) 0) xs" by (auto simp: no_trailing_unfold) thus ?thesis by (simp add: degree_eq_length_coeffs) qed auto qed lemma degree_Poly': "degree (Poly xs) ≤ length xs - 1" using length_strip_while_le[of "(=) 0" xs] by (simp add: degree_eq_length_coeffs) lemma degree_inverse_irred_power_poly_le: "degree (inverse_irred_power_poly c n) ≤ n - 1" by (auto simp: inverse_irred_power_poly_def intro: order.trans[OF degree_Poly']) lemma degree_inverse_irred_power_poly: assumes "c ≠ 0" shows "degree (inverse_irred_power_poly c n) = n - 1" unfolding inverse_irred_power_poly_def using assms by (subst degree_Poly_eq) (auto simp: last_conv_nth) lemma reflect_poly_0_iff [simp]: "reflect_poly p = 0 ⟷ p = 0" using coeff_0_reflect_poly_0_iff[of p] by fastforce lemma degree_sum_list_le: "(⋀p. p ∈ set ps ⟹ degree p ≤ T) ⟹ degree (sum_list ps) ≤ T" by (induction ps) (auto intro: degree_add_le) theorem ratfps_closed_form_exists: fixes q :: "complex poly" assumes nz: "coeff q 0 ≠ 0" defines "q' ≡ reflect_poly q" obtains r rs where "⋀n. 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)" and "⋀z. poly q' z = 0 ⟹ degree (rs z) ≤ order z q' - 1" proof - from assms have nz': "q ≠ 0" by auto from complex_alt_factorization_exists [OF nz] obtain fctrs where fctrs: "is_alt_factorization_of fctrs q" .. with nz have fctrs': "is_factorization_of fctrs q'" unfolding q'_def by (rule reflect_factorization') define r where "r = fst (solve_factored_ratfps' p fctrs)" define ts where "ts = snd (solve_factored_ratfps' p fctrs)" define rs where "rs = the ∘ map_of (map (λ(x,y). (y,x)) ts)" from nz' have "q' ≠ 0" by (simp add: q'_def) hence roots: "{z. poly q' z = 0} = set (map fst (snd fctrs))" using is_factorization_of_roots [of "fst fctrs" "snd fctrs" q'] fctrs' by simp have rs: "rs c = r" if "(r, c) ∈ set ts" for c r proof - have "map_of (map (λ(x,y). (y, x)) (snd (solve_factored_ratfps' p fctrs))) c = Some r" using that fctrs by (intro map_of_is_SomeI) (force simp: o_def case_prod_unfold solve_factored_ratfps'_def ts_def solve_factored_ratfps_roots is_alt_factorization_of_def)+ thus ?thesis by (simp add: rs_def ts_def) qed have [simp]: "length ts = length (snd fctrs)" by (auto simp: ts_def solve_factored_ratfps'_def case_prod_unfold solve_factored_ratfps_def) { fix n :: nat have "fps_of_poly p / fps_of_poly q = Abs_fps (interp_ratfps_solution (solve_factored_ratfps' p fctrs))" using solve_factored_ratfps' [OF fctrs nz'] .. also have "fps_nth … n = interp_ratfps_solution (solve_factored_ratfps' p fctrs) n" by simp also have "… = coeff r n + (∑p←snd (solve_factored_ratfps' p fctrs). poly (fst p) (of_nat n) * snd p ^ n)" (is "_ = _ + ?A") unfolding interp_ratfps_solution_def case_prod_unfold r_def by simp also have "?A = (∑p←ts. poly (rs (snd p)) (of_nat n) * snd p ^ n)" by (intro arg_cong[OF map_cong] refl) (auto simp: rs ts_def) also have "… = (∑c←map snd ts. poly (rs c) (of_nat n) * c ^ n)" by (simp add: o_def) also have "map snd ts = map fst (snd fctrs)" unfolding solve_factored_ratfps'_def case_prod_unfold ts_def by (rule solve_factored_ratfps_roots) also have "(∑c←…. poly (rs c) (of_nat n) * c ^ n) = (∑c | poly q' c = 0. poly (rs c) (of_nat n) * c ^ n)" unfolding roots using fctrs by (intro sum_list_distinct_conv_sum_set) (auto simp: is_alt_factorization_of_def) finally have "fps_nth (fps_of_poly p / fps_of_poly q) n = coeff r n + (∑c∈{z. poly q' z = 0}. poly (rs c) (of_nat n) * c ^ n)" . } moreover { fix z assume "poly q' z = 0" hence "z ∈ set (map fst (snd fctrs))" using roots by blast then obtain i where i: "i < length (snd fctrs)" and [simp]: "z = fst (snd fctrs ! i)" by (auto simp: set_conv_nth) from i have "(fst (ts ! i), snd (ts ! i)) ∈ set ts" by (auto simp: set_conv_nth) also from i have "snd (ts ! i) = z" by (simp add: ts_def solve_factored_ratfps'_def case_prod_unfold solve_factored_ratfps_def) finally have "rs z = fst (ts ! i)" by (intro rs) auto also have "… = (∑p←zip (snd (poly_pfd_simple (smult (inverse (fst fctrs)) p) (snd fctrs)) ! i) [0..<Suc (snd (snd fctrs ! i))]. inverse_irred_power_poly (fst p) (Suc (snd (snd fctrs ! i)) - snd p))" using i by (auto simp: ts_def solve_factored_ratfps'_def solve_factored_ratfps_def o_def case_prod_unfold Let_def simp del: upt_Suc power_Suc) also have "degree … ≤ snd (snd fctrs ! i)" by (intro degree_sum_list_le) (auto intro!: order.trans [OF degree_inverse_irred_power_poly_le]) also have "order z q' = Suc …" using nz' fctrs' i by (intro is_factorization_of_order[of q' "fst fctrs" "snd fctrs"]) (auto simp: q'_def) hence "snd (snd fctrs ! i) = order z q' - 1" by simp finally have "degree (rs z) ≤ …" . } ultimately show ?thesis using that[of r rs] by blast qed end  # Theory Linear_Recurrences_Common (* File: Linear_Recurrence_Common.thy Author: Manuel Eberl, TU München *) section ‹Material common to homogenous and inhomogenous linear recurrences› theory Linear_Recurrences_Common imports Complex_Main "HOL-Computational_Algebra.Computational_Algebra" begin definition lr_fps_denominator where "lr_fps_denominator cs = Poly (rev cs)" lemma lr_fps_denominator_code [code abstract]: "coeffs (lr_fps_denominator cs) = rev (dropWhile ((=) 0) cs)" by (simp add: lr_fps_denominator_def) definition lr_fps_denominator' where "lr_fps_denominator' cs = Poly cs" lemma lr_fps_denominator'_code [code abstract]: "coeffs (lr_fps_denominator' cs) = strip_while ((=) 0) cs" by (simp add: lr_fps_denominator'_def) lemma lr_fps_denominator_nz: "last cs ≠ 0 ⟹ cs ≠ [] ⟹ lr_fps_denominator cs ≠ 0" unfolding lr_fps_denominator_def by (subst coeffs_eq_iff) (auto simp: poly_eq_iff intro!: bexI[of _ "last cs"]) lemma lr_fps_denominator'_nz: "last cs ≠ 0 ⟹ cs ≠ [] ⟹ lr_fps_denominator' cs ≠ 0" unfolding lr_fps_denominator'_def by (subst coeffs_eq_iff) (auto simp: poly_eq_iff intro!: bexI[of _ "last cs"]) end  # Theory Linear_Homogenous_Recurrences (* File: Linear_Homogenous_Recurrences.thy Author: Manuel Eberl, TU München *) section ‹Homogenous linear recurrences› theory Linear_Homogenous_Recurrences imports Complex_Main RatFPS Rational_FPS_Solver Linear_Recurrences_Common begin text ‹ The following is the numerator of the rational generating function of a linear homogenous recurrence. › definition lhr_fps_numerator where "lhr_fps_numerator m cs f = (let N = length cs - 1 in Poly [(∑i≤min N k. cs ! (N - i) * f (k - i)). k ← [0..<N+m]])" lemma lhr_fps_numerator_code [code abstract]: "coeffs (lhr_fps_numerator m cs f) = (let N = length cs - 1 in strip_while ((=) 0) [(∑i≤min N k. cs ! (N - i) * f (k - i)). k ← [0..<N+m]])" by (simp add: lhr_fps_numerator_def Let_def) lemma lhr_fps_aux: fixes f :: "nat ⇒ 'a :: field" assumes "⋀n. n ≥ m ⟹ (∑k≤N. c k * f (n + k)) = 0" assumes cN: "c N ≠ 0" defines "p ≡ Poly [c (N - k). k ← [0..<Suc N]]" defines "q ≡ Poly [(∑i≤min N k. c (N - i) * f (k - i)). k ← [0..<N+m]]" shows "Abs_fps f = fps_of_poly q / fps_of_poly p" proof - include fps_notation define F where "F = Abs_fps f" have [simp]: "F$ n = f n" for n by (simp add: F_def)
have [simp]: "coeff p 0 = c N"
by (simp add: p_def nth_default_def del: upt_Suc)

have "(fps_of_poly p * F) $n = coeff q n" for n proof (cases "n ≥ N + m") case True let ?f = "λi. N - i" have "(fps_of_poly p * F)$ n = (∑i≤n. coeff p i * f (n - i))"
by (simp add: fps_mult_nth atLeast0AtMost)
also from True have "… = (∑i≤N. coeff p i * f (n - i))"
by (intro sum.mono_neutral_right) (auto simp: nth_default_def p_def)
also have "… = (∑i≤N. c (N - i) * f (n - i))"
by (intro sum.cong) (auto simp: nth_default_def p_def simp del: upt_Suc)
also from True have "… = (∑i≤N. c i * f (n - N + i))"
by (intro sum.reindex_bij_witness[of _ ?f ?f]) auto
also from True have "… = 0" by (intro assms) simp_all
also from True have "… = coeff q n"
by (simp add: q_def nth_default_def del: upt_Suc)
finally show ?thesis .
next
case False
hence "(fps_of_poly p * F) $n = (∑i≤n. coeff p i * f (n - i))" by (simp add: fps_mult_nth atLeast0AtMost) also have "… = (∑i≤min N n. coeff p i * f (n - i))" by (intro sum.mono_neutral_right) (auto simp: p_def nth_default_def simp del: upt_Suc) also have "… = (∑i≤min N n. c (N - i) * f (n - i))" by (intro sum.cong) (simp_all add: p_def nth_default_def del: upt_Suc) also from False have "… = coeff q n" by (simp add: q_def nth_default_def) finally show ?thesis . qed hence "fps_of_poly p * F = fps_of_poly q" by (intro fps_ext) simp with cN show "F = fps_of_poly q / fps_of_poly p" by (subst unit_eq_div2) (simp_all add: mult_ac) qed lemma lhr_fps: fixes f :: "nat ⇒ 'a :: field" and cs :: "'a list" defines "N ≡ length cs - 1" assumes cs: "cs ≠ []" assumes "⋀n. n ≥ m ⟹ (∑k≤N. cs ! k * f (n + k)) = 0" assumes cN: "last cs ≠ 0" shows "Abs_fps f = fps_of_poly (lhr_fps_numerator m cs f) / fps_of_poly (lr_fps_denominator cs)" proof - define p and q where "p = Poly (map (λk. ∑i≤min N k. cs ! (N - i) * f (k - i)) [0..<N + m])" and "q = Poly (map (λk. cs ! (N - k)) [0..<Suc N])" from assms have "Abs_fps f = fps_of_poly p / fps_of_poly q" unfolding p_def q_def by (intro lhr_fps_aux) (simp_all add: last_conv_nth) also have "p = lhr_fps_numerator m cs f" unfolding p_def lhr_fps_numerator_def by (auto simp: Let_def N_def) also from cN have "q = lr_fps_denominator cs" unfolding q_def lr_fps_denominator_def by (intro poly_eqI) (auto simp add: nth_default_def rev_nth N_def not_less cs simp del: upt_Suc) finally show ?thesis . qed (* TODO: Do I even need this? *) fun lhr where "lhr cs fs n = (if (cs :: 'a :: field list) = [] ∨ last cs = 0 ∨ length fs < length cs - 1 then undefined else (if n < length fs then fs ! n else (∑k<length cs - 1. cs ! k * lhr cs fs (n + 1 - length cs + k)) / -last cs))" declare lhr.simps [simp del] lemma lhr_rec: assumes "cs ≠ []" "last cs ≠ 0" "length fs ≥ length cs - 1" "n ≥ length fs" shows "(∑k<length cs. cs ! k * lhr cs fs (n + 1 - length cs + k)) = 0" proof - from assms have "{..<length cs} = insert (length cs - 1) {..<length cs - 1}" by auto also have "(∑k∈… . cs ! k * lhr cs fs (n + 1 - length cs + k)) = (∑k<length cs - 1. cs ! k * lhr cs fs (n + 1 - length cs + k)) + last cs * lhr cs fs n" using assms by (cases cs) (simp_all add: algebra_simps last_conv_nth) also from assms have "… = 0" by (subst (2) lhr.simps) (simp_all add: field_simps) finally show ?thesis . qed lemma lhrI: assumes "cs ≠ []" "last cs ≠ 0" "length fs ≥ length cs - 1" assumes "⋀n. n < length fs ⟹ f n = fs ! n" assumes "⋀n. n ≥ length fs ⟹ (∑k<length cs. cs ! k * f (n + 1 - length cs + k)) = 0" shows "f n = lhr cs fs n" using assms proof (induction cs fs n rule: lhr.induct) case (1 cs fs n) show ?case proof (cases "n < length fs") case False with 1 have "0 = (∑k<length cs. cs ! k * f (n + 1 - length cs + k))" by simp also from 1 have "{..<length cs} = insert (length cs - 1) {..<length cs - 1}" by auto also have "(∑k∈… . cs ! k * f (n + 1 - length cs + k)) = (∑k<length cs - 1. cs ! k * f (n + 1 - length cs + k)) + last cs * f n" using 1 False by (cases cs) (simp_all add: algebra_simps last_conv_nth) also have "(∑k<length cs - 1. cs ! k * f (n + 1 - length cs + k)) = (∑k<length cs - 1. cs ! k * lhr cs fs (n + 1 - length cs + k))" using False 1 by (intro sum.cong refl) simp finally have "f n = (∑k<length cs - 1. cs ! k * lhr cs fs (n + 1 - length cs + k)) / -last cs" using ‹last cs ≠ 0› by (simp add: field_simps eq_neg_iff_add_eq_0) also from 1(2-4) False have "… = lhr cs fs n" by (subst lhr.simps) simp finally show ?thesis . qed (insert 1(2-5), simp add: lhr.simps) qed (* END TODO *) locale linear_homogenous_recurrence = fixes f :: "nat ⇒ 'a :: comm_semiring_0" and cs fs :: "'a list" assumes base: "n < length fs ⟹ f n = fs ! n" assumes cs_not_null [simp]: "cs ≠ []" and last_cs [simp]: "last cs ≠ 0" and hd_cs [simp]: "hd cs ≠ 0" and enough_base: "length fs + 1 ≥ length cs" assumes rec: "n ≥ length fs - length cs ⟹ (∑k<length cs. cs ! k * f (n + k)) = 0" begin lemma lhr_fps_numerator_altdef: "lhr_fps_numerator (length fs + 1 - length cs) cs f = lhr_fps_numerator (length fs + 1 - length cs) cs ((!) fs)" proof - define N where "N = length cs - 1" define m where "m = length fs + 1 - length cs" have "lhr_fps_numerator m cs f = Poly (map (λk. (∑i≤min N k. cs ! (N - i) * f (k - i))) [0..<N + m])" by (simp add: lhr_fps_numerator_def Let_def N_def) also from enough_base have "N + m = length fs" by (cases cs) (simp_all add: N_def m_def algebra_simps) also { fix k assume k: "k ∈ {0..<length fs}" hence "f (k - i) = fs ! (k - i)" if "i ≤ min N k" for i using enough_base that by (intro base) (auto simp: Suc_le_eq N_def m_def algebra_simps) hence "(∑i≤min N k. cs ! (N - i) * f (k - i)) = (∑i≤min N k. cs ! (N - i) * fs ! (k - i))" by simp } hence "map (λk. (∑i≤min N k. cs ! (N - i) * f (k - i))) [0..<length fs] = map (λk. (∑i≤min N k. cs ! (N - i) * fs ! (k - i))) [0..<length fs]" by (intro map_cong) simp_all also have "Poly … = lhr_fps_numerator m cs ((!) fs)" using enough_base by (cases cs) (simp_all add: lhr_fps_numerator_def Let_def m_def N_def) finally show ?thesis unfolding m_def . qed end (* TODO Duplication *) lemma solve_lhr_aux: assumes "linear_homogenous_recurrence f cs fs" assumes "is_factorization_of fctrs (lr_fps_denominator' cs)" shows "f = interp_ratfps_solution (solve_factored_ratfps' (lhr_fps_numerator (length fs + 1 - length cs) cs ((!) fs)) fctrs)" proof - interpret linear_homogenous_recurrence f cs fs by fact note assms(2) hence "is_alt_factorization_of fctrs (reflect_poly (lr_fps_denominator' cs))" by (intro reflect_factorization) (simp_all add: lr_fps_denominator'_def nth_default_def hd_conv_nth [symmetric]) also have "reflect_poly (lr_fps_denominator' cs) = lr_fps_denominator cs" unfolding lr_fps_denominator_def lr_fps_denominator'_def by (subst coeffs_eq_iff) (simp add: coeffs_reflect_poly strip_while_rev [symmetric] no_trailing_unfold last_rev del: strip_while_rev) finally have factorization: "is_alt_factorization_of fctrs (lr_fps_denominator cs)" . define m where "m = length fs + 1 - length cs" obtain a ds where fctrs: "fctrs = (a, ds)" by (cases fctrs) simp_all define p and p' where "p = lhr_fps_numerator m cs ((!) fs)" and "p' = smult (inverse a) p" obtain b es where sol: "solve_factored_ratfps' p fctrs = (b, es)" by (cases "solve_factored_ratfps' p fctrs") simp_all have sol': "(b, es) = solve_factored_ratfps p' ds" by (subst sol [symmetric]) (simp add: fctrs p'_def solve_factored_ratfps_def solve_factored_ratfps'_def case_prod_unfold) have factorization': "lr_fps_denominator cs = interp_alt_factorization fctrs" using factorization by (simp add: is_alt_factorization_of_def) from assms(2) have distinct: "distinct (map fst ds)" by (simp add: fctrs is_factorization_of_def) have coeff_0_denom: "coeff (lr_fps_denominator cs) 0 ≠ 0" by (simp add: lr_fps_denominator_def nth_default_def hd_conv_nth [symmetric] hd_rev) have "coeff (lr_fps_denominator' cs) 0 ≠ 0" by (simp add: lr_fps_denominator'_def nth_default_def hd_conv_nth [symmetric]) with assms(2) have no_zero: "0 ∉ fst  set ds" by (simp add: zero_in_factorization_iff fctrs) from assms(2) have a_nz [simp]: "a ≠ 0" by (auto simp: fctrs interp_factorization_def is_factorization_of_def lr_fps_denominator'_nz) hence unit1: "is_unit (fps_const a)" by simp moreover have "is_unit (fps_of_poly (interp_alt_factorization fctrs))" by (simp add: coeff_0_denom factorization' [symmetric]) ultimately have unit2: "is_unit (fps_of_poly (∏p←ds. [:1, - fst p:] ^ Suc (snd p)))" by (simp add: fctrs case_prod_unfold interp_alt_factorization_def del: power_Suc) have "Abs_fps f = fps_of_poly (lhr_fps_numerator m cs f) / fps_of_poly (lr_fps_denominator cs)" proof (intro lhr_fps) fix n assume n: "n ≥ m" have "{..length cs - 1} = {..<length cs}" by (cases cs) auto also from n have "(∑k∈… . cs ! k * f (n + k)) = 0" by (intro rec) (simp_all add: m_def algebra_simps) finally show "(∑k≤length cs - 1. cs ! k * f (n + k)) = 0" . qed (simp_all add: m_def) also have "lhr_fps_numerator m cs f = lhr_fps_numerator m cs ((!) fs)" unfolding lhr_fps_numerator_def using enough_base by (auto simp: Let_def poly_eq_iff nth_default_def base m_def Suc_le_eq intro!: sum.cong) also have "fps_of_poly … / fps_of_poly (lr_fps_denominator cs) = fps_of_poly (lhr_fps_numerator m cs ((!) fs)) / (fps_const (fst fctrs) * fps_of_poly (∏p←snd fctrs. [:1, - fst p:] ^ Suc (snd p)))" unfolding assms factorization' interp_alt_factorization_def by (simp add: case_prod_unfold Let_def fps_of_poly_smult) also from unit1 unit2 have "… = fps_of_poly p / fps_const a / fps_of_poly (∏(c,n)←ds. [:1, -c:]^Suc n)" by (subst is_unit_div_mult2_eq) (simp_all add: fctrs case_prod_unfold p_def) also from unit1 have "fps_of_poly p / fps_const a = fps_of_poly p'" by (simp add: fps_divide_unit fps_of_poly_smult fps_const_inverse p'_def) also from distinct no_zero have "… / fps_of_poly (∏(c,n)←ds. [:1, -c:]^Suc n) = Abs_fps (interp_ratfps_solution (solve_factored_ratfps' p fctrs))" by (subst solve_factored_ratfps) (simp_all add: case_prod_unfold sol' sol) finally show ?thesis unfolding p_def m_def by (intro ext) (simp add: fps_eq_iff) qed definition "lhr_fps as fs = ( let m = length fs + 1 - length as; p = lhr_fps_numerator m as (λn. fs ! n); q = lr_fps_denominator as in ratfps_of_poly p / ratfps_of_poly q)" lemma lhr_fps_correct: fixes f :: "nat ⇒ 'a :: {field_char_0,field_gcd}" assumes "linear_homogenous_recurrence f cs fs" shows "fps_of_ratfps (lhr_fps cs fs) = Abs_fps f" proof - interpret linear_homogenous_recurrence f cs fs by fact define m where "m = length fs + 1 - length cs" let ?num = "lhr_fps_numerator m cs f" let ?num' = "lhr_fps_numerator m cs ((!) fs)" let ?denom = "lr_fps_denominator cs" have "{..length cs - 1} = {..<length cs}" by (cases cs) auto moreover have "length cs ≥ 1" by (cases cs) auto ultimately have "Abs_fps f = fps_of_poly ?num / fps_of_poly ?denom" by (intro lhr_fps) (insert rec, simp_all add: m_def) also have "?num = ?num'" by (rule lhr_fps_numerator_altdef [folded m_def]) also have "fps_of_poly ?num' / fps_of_poly ?denom = fps_of_ratfps (ratfps_of_poly ?num' / ratfps_of_poly ?denom)" by simp also from enough_base have "… = fps_of_ratfps (lhr_fps cs fs)" by (cases cs) (simp_all add: base fps_of_ratfps_def case_prod_unfold lhr_fps_def m_def) finally show ?thesis .. qed end  # Theory Eulerian_Polynomials (* File: Eulerian_Polynomials.thy Author: Manuel Eberl, TU München *) section ‹Eulerian polynomials› theory Eulerian_Polynomials imports Complex_Main "HOL-Combinatorics.Stirling" "HOL-Computational_Algebra.Computational_Algebra" begin text ‹ The Eulerian polynomials are a sequence of polynomials that is related to the closed forms of the power series $\sum_{n=0}^\infty n^k X^n$ for a fixed$k$. › primrec eulerian_poly :: "nat ⇒ 'a :: idom poly" where "eulerian_poly 0 = 1" | "eulerian_poly (Suc n) = (let p = eulerian_poly n in [:0,1,-1:] * pderiv p + p * [:1, of_nat n:])" lemmas eulerian_poly_Suc [simp del] = eulerian_poly.simps(2) lemma eulerian_poly: "fps_of_poly (eulerian_poly k :: 'a :: field poly) = Abs_fps (λn. of_nat (n+1) ^ k) * (1 - fps_X) ^ (k + 1)" proof (induction k) case 0 have "Abs_fps (λ_. 1 :: 'a) = inverse (1 - fps_X)" by (rule fps_inverse_unique [symmetric]) (simp add: inverse_mult_eq_1 fps_inverse_gp' [symmetric]) thus ?case by (simp add: inverse_mult_eq_1) next case (Suc k) define p :: "'a fps" where "p = fps_of_poly (eulerian_poly k)" define F :: "'a fps" where "F = Abs_fps (λn. of_nat (n+1) ^ k)" have p: "p = F * (1 - fps_X) ^ (k+1)" by (simp add: p_def Suc F_def) have p': "fps_deriv p = fps_deriv F * (1 - fps_X) ^ (k + 1) - F * (1 - fps_X) ^ k * of_nat (k + 1)" by (simp add: p fps_deriv_power algebra_simps fps_const_neg [symmetric] fps_of_nat del: power_Suc of_nat_Suc fps_const_neg) have "fps_of_poly (eulerian_poly (Suc k)) = (fps_X * fps_deriv F + F) * (1 - fps_X) ^ (Suc k + 1)" apply (simp add: Let_def p_def [symmetric] fps_of_poly_simps eulerian_poly_Suc del: power_Suc) apply (simp add: p p' fps_deriv_power fps_const_neg [symmetric] fps_of_nat del: power_Suc of_nat_Suc fps_const_neg) apply (simp add: algebra_simps) done also have "fps_X * fps_deriv F + F = Abs_fps (λn. of_nat (n + 1) ^ Suc k)" unfolding F_def by (intro fps_ext) (auto simp: algebra_simps) finally show ?case . qed lemma eulerian_poly': "Abs_fps (λn. of_nat (n+1) ^ k) = fps_of_poly (eulerian_poly k :: 'a :: field poly) / (1 - fps_X) ^ (k + 1)" by (subst eulerian_poly) simp lemma eulerian_poly'': assumes k: "k > 0" shows "Abs_fps (λn. of_nat n ^ k) = fps_of_poly (pCons 0 (eulerian_poly k :: 'a :: field poly)) / (1 - fps_X) ^ (k + 1)" proof - from assms have "Abs_fps (λn. of_nat n ^ k :: 'a) = fps_X * Abs_fps (λn. of_nat (n + 1) ^ k)" by (intro fps_ext) (auto simp: of_nat_diff) also have "Abs_fps (λn. of_nat (n + 1) ^ k :: 'a) = fps_of_poly (eulerian_poly k) / (1 - fps_X) ^ (k + 1)" by (rule eulerian_poly') also have "fps_X * … = fps_of_poly (pCons 0 (eulerian_poly k)) / (1 - fps_X) ^ (k + 1)" by (simp add: fps_of_poly_pCons fps_divide_unit) finally show ?thesis . qed definition fps_monom_poly :: "'a :: field ⇒ nat ⇒ 'a poly" where "fps_monom_poly c k = (if k = 0 then 1 else pcompose (pCons 0 (eulerian_poly k)) [:0,c:])" primrec fps_monom_poly_aux :: "'a :: field ⇒ nat ⇒ 'a poly" where "fps_monom_poly_aux c 0 = [:c:]" | "fps_monom_poly_aux c (Suc k) = (let p = fps_monom_poly_aux c k in [:0,1,-c:] * pderiv p + [:1, of_nat k * c:] * p)" lemma fps_monom_poly_aux: "fps_monom_poly_aux c k = smult c (pcompose (eulerian_poly k) [:0,c:])" by (induction k) (simp_all add: eulerian_poly_Suc Let_def pderiv_pcompose pcompose_pCons pcompose_add pcompose_smult pcompose_uminus smult_add_right pderiv_pCons pderiv_smult algebra_simps one_pCons) lemma fps_monom_poly_code [code]: "fps_monom_poly c k = (if k = 0 then 1 else pCons 0 (fps_monom_poly_aux c k))" by (simp add: fps_monom_poly_def fps_monom_poly_aux pcompose_pCons) lemma fps_monom_aux: "Abs_fps (λn. of_nat n ^ k) = fps_of_poly (fps_monom_poly 1 k) / (1 - fps_X) ^ (k+1)" proof (cases "k = 0") assume [simp]: "k = 0" hence "Abs_fps (λn. of_nat n ^ k :: 'a) = Abs_fps (λ_. 1)" by simp also have "… = 1 / (1 - fps_X)" by (subst gp [symmetric]) simp_all finally show ?thesis by (simp add: fps_monom_poly_def) qed (insert eulerian_poly''[of k, where ?'a = 'a], simp add: fps_monom_poly_def) lemma fps_monom: "Abs_fps (λn. of_nat n ^ k * c ^ n) = fps_of_poly (fps_monom_poly c k) / (1 - fps_const c * fps_X) ^ (k+1)" proof - have "Abs_fps (λn. of_nat n ^ k * c ^ n) = fps_compose (Abs_fps (λn. of_nat n ^ k)) (fps_const c * fps_X)" by (subst fps_compose_linear) (simp add: mult_ac) also have "Abs_fps (λn. of_nat n ^ k) = fps_of_poly (fps_monom_poly 1 k) / (1 - fps_X) ^ (k+1)" by (rule fps_monom_aux) also have "fps_compose … (fps_const c * fps_X) = (fps_of_poly (fps_monom_poly 1 k) oo fps_const c * fps_X) / ((1 - fps_X) ^ (k + 1) oo fps_const c * fps_X)" by (intro fps_compose_divide_distrib) (simp_all add: fps_compose_power [symmetric] fps_compose_sub_distrib del: power_Suc) also have "fps_of_poly (fps_monom_poly 1 k) oo (fps_const c * fps_X) = fps_of_poly (fps_monom_poly c k)" by (simp add: fps_monom_poly_def fps_of_poly_pcompose fps_of_poly_simps fps_of_poly_pCons mult_ac) also have "((1 - fps_X) ^ (k + 1) oo fps_const c * fps_X) = (1 - fps_const c * fps_X) ^ (k + 1)" by (simp add: fps_compose_power [symmetric] fps_compose_sub_distrib del: power_Suc) finally show ?thesis . qed end  # Theory Linear_Inhomogenous_Recurrences (* File: Linear_Inhomogenous_Recurrences.thy Author: Manuel Eberl, TU München *) section ‹Inhomogenous linear recurrences› theory Linear_Inhomogenous_Recurrences imports Complex_Main Linear_Homogenous_Recurrences Eulerian_Polynomials RatFPS begin definition lir_fps_numerator where "lir_fps_numerator m cs f g = (let N = length cs - 1 in Poly [(∑i≤min N k. cs ! (N - i) * f (k - i)) - g k. k ← [0..<N+m]])" lemma lir_fps_numerator_code [code abstract]: "coeffs (lir_fps_numerator m cs f g) = (let N = length cs - 1 in strip_while ((=) 0) [(∑i≤min N k. cs ! (N - i) * f (k - i)) - g k. k ← [0..<N+m]])" by (simp add: lir_fps_numerator_def Let_def) locale linear_inhomogenous_recurrence = fixes f g :: "nat ⇒ 'a :: comm_ring" and cs fs :: "'a list" assumes base: "n < length fs ⟹ f n = fs ! n" assumes cs_not_null [simp]: "cs ≠ []" and last_cs [simp]: "last cs ≠ 0" and hd_cs [simp]: "hd cs ≠ 0" and enough_base: "length fs + 1 ≥ length cs" assumes rec: "n ≥ length fs + 1 - length cs ⟹ (∑k<length cs. cs ! k * f (n + k)) = g (n + length cs - 1)" begin lemma coeff_0_lr_fps_denominator [simp]: "coeff (lr_fps_denominator cs) 0 = last cs" by (auto simp: lr_fps_denominator_def nth_default_def nth_Cons hd_conv_nth [symmetric] hd_rev) lemma lir_fps_numerator_altdef: "lir_fps_numerator (length fs + 1 - length cs) cs f g = lir_fps_numerator (length fs + 1 - length cs) cs ((!) fs) g" proof - define N where "N = length cs - 1" define m where "m = length fs + 1 - length cs" have "lir_fps_numerator m cs f g = Poly (map (λk. (∑i≤min N k. cs ! (N - i) * f (k - i)) - g k) [0..<N + m])" by (simp add: lir_fps_numerator_def Let_def N_def) also from enough_base have "N + m = length fs" by (cases cs) (simp_all add: N_def m_def algebra_simps) also { fix k assume k: "k ∈ {0..<length fs}" hence "f (k - i) = fs ! (k - i)" if "i ≤ min N k" for i using enough_base that by (intro base) (auto simp: Suc_le_eq N_def m_def algebra_simps) hence "(∑i≤min N k. cs ! (N - i) * f (k - i)) = (∑i≤min N k. cs ! (N - i) * fs ! (k - i))" by simp } hence "map (λk. (∑i≤min N k. cs ! (N - i) * f (k - i)) - g k) [0..<length fs] = map (λk. (∑i≤min N k. cs ! (N - i) * fs ! (k - i)) - g k) [0..<length fs]" by (intro map_cong) simp_all also have "Poly … = lir_fps_numerator m cs ((!) fs) g" using enough_base by (cases cs) (simp_all add: lir_fps_numerator_def Let_def m_def N_def) finally show ?thesis unfolding m_def . qed end context begin private lemma lir_fps_aux: fixes f :: "nat ⇒ 'a :: field" assumes rec: "⋀n. n ≥ m ⟹ (∑k≤N. c k * f (n + k)) = g (n + N)" assumes cN: "c N ≠ 0" defines "p ≡ Poly [c (N - k). k ← [0..<Suc N]]" defines "q ≡ Poly [(∑i≤min N k. c (N - i) * f (k - i)) - g k. k ← [0..<N+m]]" shows "Abs_fps f = (fps_of_poly q + Abs_fps g) / fps_of_poly p" proof - include fps_notation define F where "F = Abs_fps f" have [simp]: "F$ n = f n" for n by (simp add: F_def)
have [simp]: "coeff p 0 = c N"
by (simp add: p_def nth_default_def del: upt_Suc)

have "(fps_of_poly p * F) $n = coeff q n + g n" for n proof (cases "n ≥ N + m") case True let ?f = "λi. N - i" have "(fps_of_poly p * F)$ n = (∑i≤n. coeff p i * f (n - i))"
by (simp add: fps_mult_nth atLeast0AtMost)
also from True have "… = (∑i≤N. coeff p i * f (n - i))"
by (intro sum.mono_neutral_right) (auto simp: nth_default_def p_def)
also have "… = (∑i≤N. c (N - i) * f (n - i))"
by (intro sum.cong) (auto simp: nth_default_def p_def simp del: upt_Suc)
also from True have "… = (∑i≤N. c i * f (n - N + i))"
by (intro sum.reindex_bij_witness[of _ ?f ?f]) auto
also from True have "… = g (n - N + N)" by (intro rec) simp_all
also from True have "… = coeff q n + g n"
by (simp add: q_def nth_default_def del: upt_Suc)
finally show ?thesis .
next
case False
hence "(fps_of_poly p * F) \$ n = (∑i≤n. coeff p i * f (n - i))"
by (simp add: fps_mult_nth atLeast0AtMost)
also have "… = (∑i≤min N n. coeff p i * f (n - i))"
by (intro sum.mono_neutral_right)
(auto simp: p_def nth_default_def simp del: upt_Suc)
also have "… = (∑i≤min N n. c (N - i) * f (n - i))"
by (intro sum.cong) (simp_all add: p_def nth_default_def del: upt_Suc)
also from False have "… = coeff q n + g n" by (simp add: q_def nth_default_def)
finally show ?thesis .
qed
hence "fps_of_poly p * F = fps_of_poly q + Abs_fps g"
by (intro fps_ext) (simp add:)
with cN show "F = (fps_of_poly q + Abs_fps g) / fps_of_poly p"
by (subst unit_eq_div2) (simp_all add: mult_ac)
qed

lemma lir_fps:
fixes f g :: "nat ⇒ 'a :: field" and cs :: "'a list"
defines "N ≡ length cs - 1"
assumes cs: "cs ≠ []"
assumes "⋀n. n ≥ m ⟹ (∑k≤N. cs ! k * f (n + k)) = g (n + N)"
assumes cN: "last cs ≠ 0"
shows   "Abs_fps f = (fps_of_poly (lir_fps_numerator m cs f g) + Abs_fps g) /
fps_of_poly (lr_fps_denominator cs)"
proof -
define p and q
where "p = Poly [(∑i≤min N k. cs ! (N - i) * f (k - i)) - g k. k ← [0..<N+m]]"
and "q = Poly (map (λk. cs ! (N - k)) [0..<Suc N])"
from assms have "Abs_fps f = (fps_of_poly p + Abs_fps g) / fps_of_poly q"
unfolding p_def q_def by (intro lir_fps_aux) (simp_all add: last_conv_nth)
also have "p = lir_fps_numerator m cs f g"
unfolding p_def lir_fps_numerator_def by (auto simp: Let_def N_def)
also from cN have "q = lr_fps_denominator cs"
unfolding q_def lr_fps_denominator_def
by (intro poly_eqI)
(auto simp add: nth_default_def rev_nth N_def not_less cs simp del: upt_Suc)
finally show ?thesis .
qed

end

type_synonym 'a polyexp = "('a × nat × 'a) list"

definition eval_polyexp :: "('a::semiring_1) polyexp ⇒ nat ⇒ 'a" where
"eval_polyexp xs = (λn. ∑(a,k,b)←xs. a * of_nat n ^ k * b ^ n)"

lemma eval_polyexp_Nil [simp]: "eval_polyexp [] = (λ_. 0)"
by (simp add: eval_polyexp_def)

lemma eval_polyexp_Cons:
"eval_polyexp (x#xs) = (λn. (case x of (a,k,b) ⇒ a * of_nat n ^ k * b ^ n) + eval_polyexp xs n)"
by (simp add: eval_polyexp_def)

definition polyexp_fps :: "('a :: field) polyexp ⇒ 'a fps" where
"polyexp_fps xs =
(∑(a,k,b)←xs. fps_of_poly (Polynomial.smult a (fps_monom_poly b k)) /
(1 - fps_const b * fps_X) ^ (k + 1))"

lemma polyexp_fps_Nil [simp]: "polyexp_fps [] = 0"
by (simp add: polyexp_fps_def)

lemma polyexp_fps_Cons:
"polyexp_fps (x#xs) = (case x of (a,k,b) ⇒
fps_of_poly (Polynomial.smult a (fps_monom_poly b k)) / (1 - fps_const b * fps_X) ^ (k + 1)) +
polyexp_fps xs"
by (simp add: polyexp_fps_def)

definition polyexp_ratfps :: "('a :: field_gcd) polyexp ⇒ 'a ratfps" where
"polyexp_ratfps xs =
(∑(a,k,b)←xs. ratfps_of_poly (Polynomial.smult a (fps_monom_poly b k)) /
ratfps_of_poly ([:1, -b:] ^ (k + 1)))"

lemma polyexp_ratfps_Nil [simp]: "polyexp_ratfps [] = 0"
by (simp add: polyexp_ratfps_def)

lemma polyexp_ratfps_Cons: "polyexp_ratfps (x#xs) = (case x of (a,k,b) ⇒
ratfps_of_poly (Polynomial.smult a (fps_monom_poly b k)) /
ratfps_of_poly ([:1, -b:] ^ (k + 1))) + polyexp_ratfps xs"
by (simp add: polyexp_ratfps_def)

lemma polyexp_fps: "Abs_fps (eval_polyexp xs) = polyexp_fps xs"
proof (induction xs)
case (Cons x xs)
obtain a k b where [simp]: "x = (a, k, b)" by (metis prod.exhaust)
have "Abs_fps (eval_polyexp (x#xs)) =
fps_const a * Abs_fps (λn. of_nat n ^ k * b ^ n) + Abs_fps (eval_polyexp xs)"
by (simp add: eval_polyexp_Cons fps_plus_def mult_ac)
also have "Abs_fps (λn. of_nat n ^ k * b ^ n) =
fps_of_poly (fps_monom_poly b k) / (1 - fps_const b * fps_X) ^ (k + 1)"
(is "_ = ?A / ?B")
by (rule fps_monom)
also have "fps_const a * (?A / ?B) = (fps_const a * ?A) / ?B"
by (intro unit_div_mult_swap) simp_all
also have "fps_const a * ?A = fps_of_poly (Polynomial.smult a (fps_monom_poly b k))"
by simp
also note Cons.IH
finally show ?case by (simp add: polyexp_fps_Cons)
qed (simp_all add: fps_zero_def)

lemma polyexp_ratfps [simp]: "fps_of_ratfps (polyexp_ratfps xs) = polyexp_fps xs"
by (induction xs)
(auto simp del: power_Suc fps_const_neg
simp: coeff_0_power fps_of_poly_power fps_of_poly_smult fps_of_poly_pCons
fps_const_neg [symmetric] mult_ac polyexp_ratfps_Cons polyexp_fps_Cons)

definition lir_fps ::
"'a :: field_gcd list ⇒ 'a list ⇒ 'a polyexp ⇒ ('a ratfps) option" where
"lir_fps cs fs g = (if cs = [] ∨ length fs < length cs - 1 then None else
let m = length fs + 1 - length cs;
p = lir_fps_numerator m cs (λn. fs ! n) (eval_polyexp g);
q = lr_fps_denominator cs
in  Some ((ratfps_of_poly p + polyexp_ratfps g) * inverse (ratfps_of_poly q)))"

lemma lir_fps_correct:
fixes   f :: "nat ⇒ 'a :: field_gcd"
assumes "linear_inhomogenous_recurrence f (eval_polyexp g) cs fs"
shows   "map_option fps_of_ratfps (lir_fps cs fs g) = Some (Abs_fps f)"
proof -
interpret linear_inhomogenous_recurrence f "eval_polyexp g" cs fs by fact
define m where "m = length fs + 1 - length cs"
let ?num = "lir_fps_numerator m cs f (eval_polyexp g)"
let ?num' = "lir_fps_numerator m cs ((!) fs) (eval_polyexp g)"
let ?denom = "lr_fps_denominator cs"

have "{..length cs - 1} = {..<length cs}" by (cases cs) auto
moreover have "length cs ≥ 1" by (cases cs) auto
ultimately have "Abs_fps f = (fps_of_poly ?num + Abs_fps (eval_polyexp g)) / fps_of_poly ?denom"
by (intro lir_fps) (insert rec, simp_all add: m_def)
also have "?num = ?num'" by (rule lir_fps_numerator_altdef [folded m_def])
also have "(fps_of_poly ?num' + Abs_fps (eval_polyexp g)) / fps_of_poly ?denom =
fps_of_ratfps ((ratfps_of_poly ?num' + polyexp_ratfps g) *
inverse (ratfps_of_poly ?denom))"
by (simp add: polyexp_fps fps_divide_unit)
also from enough_base have "Some … = map_option fps_of_ratfps (lir_fps cs fs g)"
by (cases cs) (simp_all add: base fps_of_ratfps_def case_prod_unfold lir_fps_def m_def)
finally show ?thesis ..
qed

end


# Theory Rational_FPS_Asymptotics

(*
File:    Rational_FPS_Asymptotics.thy
Author:  Manuel Eberl, TU München
*)
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)"
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] guess r rs . note closed_form = this

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 ―‹The case of a root at the border of the disc›
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  ―‹The case of a root in the interior of the disc›
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 > 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 ^ Suc l)"
unfolding square_free_factorization_def prod.case by blast
with ‹q ≠ 0› have [simp]: "b ≠ 0" by auto

from assms(1) have [simp]: "(0, x) ∉ set cs" for x
by (auto simp: square_free_factorization_def)
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 ^ Suc l)"
by (simp add: order_smult q)
also have "… = (∑x∈set cs. order z (case x of (c, l) ⇒ c ^ Suc l))"
by (subst order_prod) (auto dest: coprime)
also have "… = (∑(c, l)∈set cs. Suc 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. Suc 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. Suc l * order z c) = 0"
by (intro sum.neutral) auto
finally show "order z q ≤ Suc k" by simp
next
case True ―‹The order of a root is determined by the unique polynomial in the
square-free factorisation that contains it.›
then obtain c0 l0 where cl0: "(c0, l0) ∈ set cs" "poly c0 z = 0"
by blast
have "order z q = (∑(c, l)∈set cs. Suc l * order z c)" by fact
also have "… = Suc l0 * order z c0 + (∑(c, l) ∈ set cs - {(c0, l0)}. Suc l * order z c)"
using cl0 by (subst sum.remove[of _ "(c0, l0)"]) auto
also have "(∑(c, l) ∈ set cs - {(c0, l0)}. Suc 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 = Suc l0" by simp

also from roots2[of c0 l0] cl0 z have "l0 ≤ k"
by (cases l0 k rule: linorder_cases) auto
finally show "order z q ≤ Suc k" by simp
qed
qed fact+
qed

find_consts name:"Count_Complex"
term proots_ball_card
term proots_sphere_card

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 > 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)])
from assms(1) have q: "q = smult b (∏(c, l)∈set cs. c ^ Suc l)"
unfolding square_free_factorization_def prod.case by blast
with ‹q ≠ 0› have [simp]: "b ≠ 0" by auto
from assms(1) have [simp]: "(0, x) ∉ set cs" for x
by (auto simp: square_free_factorization_def)

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 > 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 > 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