Session Linear_Recurrences

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

lemma coeff_0_add_fract_nonzero:
  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 "(+)"
  by (rule coeff_0_add_fract_nonzero)

lift_definition minus_ratfps :: "'a ratfps  'a ratfps  'a ratfps" is "(-)"
  by (simp only: diff_conv_add_uminus, rule coeff_0_add_fract_nonzero)
     (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 fA-{0}. ratfps_subdegree f))"
proof (rule sym, rule GcdI)
  fix f assume "f  A"
  thus "ratfps_of_poly (monom 1 (INF fA - {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 fA-{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 fA-{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 fA-{0}. ratfps_subdegree f)))"
  using ratfps_Gcd by auto

lemma ratfps_Lcm:
  assumes "A  {}" "0  A" "bdd_above (ratfps_subdegree`A)"
  shows   "Lcm A = ratfps_of_poly (monom 1 (SUP fA. 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 fA. 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 fA. 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 fA. 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_subdegree`A) then 0 else
      if A = {} then 1 else ratfps_of_poly (monom 1 (SUP fA. ratfps_subdegree f)))"
proof (cases "bdd_above (ratfps_subdegree`A)")
  assume unbounded: "¬bdd_above (ratfps_subdegree`A)"
  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 "yset 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. jsnd (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. (jsnd (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 = (jsnd (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) (xsnd 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) (xsnd 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) = (xA. 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) = (xA. 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 " = (xfctrs. 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 " = (xset 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) (xxs. [:- x, 1:] ^ Suc (order x p - 1))"
    by (simp add: interp_factorization_def o_def)
  also have "(xxs. [:- 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 (jsnd (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'. jsnd (?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. jsnd (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. jsnd (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. jsnd (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 (jsnd (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 + (psnd (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 = (pts. 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 " = (cmap 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 " = (pzip (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 [(imin 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) [(imin 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  (kN. c k * f (n + k)) = 0"
  assumes cN: "c N  0"
  defines "p  Poly [c (N - k). k  [0..<Suc N]]"
  defines "q  Poly [(imin 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 = (in. coeff p i * f (n - i))"
      by (simp add: fps_mult_nth atLeast0AtMost)
    also from True have " = (iN. coeff p i * f (n - i))"
      by (intro sum.mono_neutral_right) (auto simp: nth_default_def p_def)
    also have " = (iN. 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 " = (iN. 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 = (in. coeff p i * f (n - i))"
      by (simp add: fps_mult_nth atLeast0AtMost)
    also have " = (imin 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 " = (imin 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  (kN. 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. imin 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. (imin 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 "(imin N k. cs ! (N - i) * f (k - i)) = (imin N k. cs ! (N - i) * fs ! (k - i))"
      by simp
  }
  hence "map (λk. (imin N k. cs ! (N - i) * f (k - i))) [0..<length fs] =
           map (λk. (imin 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 (pds. [: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 "(klength 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 (psnd 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 [(imin 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) [(imin 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. (imin 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 "(imin N k. cs ! (N - i) * f (k - i)) = (imin N k. cs ! (N - i) * fs ! (k - i))"
      by simp
  }
  hence "map (λk. (imin N k. cs ! (N - i) * f (k - i)) - g k) [0..<length fs] =
           map (λk. (imin 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  (kN. 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 [(imin 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 = (in. coeff p i * f (n - i))"
      by (simp add: fps_mult_nth atLeast0AtMost)
    also from True have " = (iN. coeff p i * f (n - i))"
      by (intro sum.mono_neutral_right) (auto simp: nth_default_def p_def)
    also have " = (iN. 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 " = (iN. 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 = (in. coeff p i * f (n - i))"
      by (simp add: fps_mult_nth atLeast0AtMost)
    also have " = (imin 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 " = (imin 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  (kN. 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 [(imin 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)"
      proof (subst asymp_equiv_add_left)
        have "(λx. x * poly p x) ∼[F] (λx. x * (lead_coeff p * x ^ degree p))"
          by (intro asymp_equiv_intros *)
        also have " = (λx. lead_coeff p * x ^ n)" by (simp add: n_def mult_ac)
        finally show "(λx. x * poly p x) ∼[F] " .
      next
        have "filterlim (λx. x) at_infinity F"
          by (simp add: filterlim_def assms)
        hence "(λx. x ^ n)  ω[F](λ_. 1 :: 'a)" unfolding smallomega_1_conv_filterlim
          by (intro Limits.filterlim_power_at_infinity filterlim_ident) (auto simp: n_def)
        hence "(λx. a)  o[F](λx. x ^ n)" unfolding smallomega_iff_smallo[symmetric]
          by (cases "a = 0") auto
        thus "(λx. a)  o[F](λx. lead_coeff p * x ^ n)"
          by simp
      qed
      ultimately show ?thesis by (simp add: n_def)
    qed auto
  qed auto
qed

lemma poly_bigtheta:
  assumes "p  0" and "F  at_infinity"
  shows   "poly p  Θ[F](λx. x ^ degree p)"
proof -
  have "poly p ∼[F] (λx. lead_coeff p * x ^ degree p)"
    by (intro poly_asymp_equiv assms)
  thus ?thesis using assms by (auto dest!: asymp_equiv_imp_bigtheta)
qed

lemma poly_bigo:
  assumes "F  at_infinity" and "degree p  k"
  shows   "poly p  O[F](λx. x ^ k)"
proof (cases "p = 0")
  case True
  hence "poly p = (λ_. 0)" by (auto simp: fun_eq_iff)
  thus ?thesis by simp
next
  case False
  have *: "(λx. x ^ (k - degree p))  Ω[F](λx. 1)"
  proof (cases "k = degree p")
    case False
    hence "(λx. x ^ (k - degree p))  ω[F](λ_. 1)"
      unfolding smallomega_1_conv_filterlim using assms False
      by (intro Limits.filterlim_power_at_infinity filterlim_ident)
         (auto simp: filterlim_def)
    thus ?thesis by (rule landau_omega.small_imp_big)
  qed auto

  have "poly p  Θ[F](λx. x ^ degree p * 1)"
    using poly_bigtheta[OF False assms(1)] by simp
  also have "(λx. x ^ degree p * 1)  O[F](λx. x ^ degree p * x ^ (k - degree p))" using *
    by (intro landau_o.big.mult landau_o.big_refl) (auto simp: bigomega_iff_bigo)
  also have "(λx::'a. x ^ degree p * x ^ (k - degree p)) = (λx. x ^ k)"
    using assms by (simp add: power_add [symmetric])
  finally show ?thesis .
qed

lemma reflect_poly_dvdI:
  fixes p q :: "'a::{comm_semiring_1,semiring_no_zero_divisors} poly"
  assumes "p dvd q"
  shows   "reflect_poly p dvd reflect_poly q"
  using assms by (auto simp: reflect_poly_mult)

lemma smult_altdef: "smult c p = [:c:] * p"
  by (induction p) (auto simp: mult_ac)

lemma smult_power: "smult (c ^ n) (p ^ n) = (smult c p) ^ n"
proof -
  have "smult (c ^ n) (p ^ n) = [:c ^ n:] * p ^ n"
    by simp
  also have "[:c:] ^ n = [:c ^ n:]"
    by (induction n) (auto simp: mult_ac)
  hence "[:c ^ n:] = [:c:] ^ n" ..
  also have " * p ^ n = ([:c:] * p) ^ n"
    by (rule power_mult_distrib [symmetric])
  also have " = (smult c p) ^ n" by simp
  finally show ?thesis .
qed

lemma order_reflect_poly_ge:
  fixes c :: "'a :: field"
  assumes "c  0" and "p  0"
  shows   "order c (reflect_poly p)  order (1 / c) p"
proof -
  have "reflect_poly ([:-(1 / c), 1:] ^ order (1 / c) p) dvd reflect_poly p"
    by (intro reflect_poly_dvdI, subst order_divides) auto
  also have "reflect_poly ([:-(1 / c), 1:] ^ order (1 / c) p) =
               smult ((-1 / c) ^ order (1 / c) p) ([:-c, 1:] ^ order (1 / c) p)"
    using assms by (simp add: reflect_poly_power reflect_poly_pCons smult_power)
  finally have "([:-c, 1:] ^ order (1 / c) p) dvd reflect_poly p"
    by (rule smult_dvd_cancel)
  with p  0 show ?thesis by (subst (asm) order_divides) auto
qed

lemma order_reflect_poly:
  fixes c :: "'a :: field"
  assumes "c  0" and "coeff p 0  0"
  shows   "order c (reflect_poly p) = order (1 / c) p"
proof (rule antisym)
  from assms show "order c (reflect_poly p)  order (1 / c) p"
    by (intro order_reflect_poly_ge) auto
next
  from assms have "order (1 / (1 / c)) (reflect_poly p) 
                     order (1 / c) (reflect_poly (reflect_poly p))"
    by (intro order_reflect_poly_ge) auto
  with assms show "order c (reflect_poly p)  order (1 / c) p"
    by simp
qed

lemma poly_reflect_eq_0_iff:
  "poly (reflect_poly p) (x :: 'a :: field) = 0  p = 0  x  0  poly p (1 / x) = 0"
  by (cases "x = 0") (auto simp: poly_reflect_poly_nz inverse_eq_divide)
  

theorem ratfps_nth_bigo:
  fixes q :: "complex poly"
  assumes "R > 0" 
  assumes roots1: "z. z  ball 0 (1 / R)  poly q z  0"
  assumes roots2: "z. z  sphere 0 (1 / R)  poly q z = 0  order z q  Suc k"
  shows   "fps_nth (fps_of_poly p / fps_of_poly q)  O(λn. of_nat n ^ k * of_real R ^ n)"
proof -
  define q' where "q' = reflect_poly q"
  from roots1[of 0] and R > 0 have [simp]: "coeff q 0  0" "q  0"
    by (auto simp: poly_0_coeff_0)
  from ratfps_closed_form_exists[OF this(1), of p] 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  xball 0 (1 / R). poly c x  0"
  assumes roots2: "c l. (c, l)  set cs  l > k  xsphere 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 " = (xset 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  (xA. 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 "xball 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 "xsphere 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