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]
  obtain r rs where closed_form:
      "n. (fps_of_poly p / fps_of_poly q) $ n =
        coeff r n + (c | poly (reflect_poly q) c = 0. poly (rs c) (of_nat n) * c ^ n)"
      "z. poly (reflect_poly q) z = 0  degree (rs z)  order z (reflect_poly q) - 1"
    by blast

  have "fps_nth (fps_of_poly p / fps_of_poly q) =
          (λn. coeff r n + (c | poly q' c = 0. poly (rs c) (of_nat n) * c ^ n))"
    by (intro ext, subst closed_form) (simp_all add: q'_def)
  also have "  O(λn. of_nat n ^ k * of_real R ^ n)"
  proof (intro sum_in_bigo big_sum_in_bigo)
    have "eventually (λn. coeff r n = 0) at_top"
      using MOST_nat coeff_eq_0 cofinite_eq_sequentially by force
    hence "coeff r  Θ(λ_. 0)" by (rule bigthetaI_cong)
    also have "(λ_. 0 :: complex)  O(λn. of_nat n ^ k * of_real R ^ n)"
      by simp
    finally show "coeff r  O(λn. of_nat n ^ k * of_real R ^ n)" .
  next
    fix c assume c: "c  {c. poly q' c = 0}"
    hence [simp]: "c  0" by (auto simp: q'_def)

    show "(λn. poly (rs c) n * c ^ n)  O(λn. of_nat n ^ k * of_real R ^ n)"
    proof (cases "norm c = R")
      case True ―‹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 > Suc 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 ^ l)"
    unfolding square_free_factorization_def prod.case by blast 
  with q  0 have [simp]: "b  0" by auto
  note sff = square_free_factorizationD[OF assms(1)]

  from sff(2)[of 0] have [simp]: "(0, x)  set cs" for x by auto

  from assms(1) have coprime: "c1 = c2" "m = n"
    if "¬coprime c1 c2" "(c1, m)  set cs" "(c2, n)  set cs" for c1 c2 m n
    using that by (auto simp: square_free_factorization_def case_prod_unfold)

  show ?thesis
  proof (rule ratfps_nth_bigo)
    fix z :: complex assume z: "z  ball 0 (1 / R)"
    show "poly q z  0"
    proof
      assume "poly q z = 0"
      then obtain c l where cl: "(c, l)  set cs" and "poly c z = 0"
        by (auto simp: q poly_prod image_iff)
      with roots1[of c l] and z show False by auto
    qed
  next
    fix z :: complex assume z: "z  sphere 0 (1 / R)"

    have order: "order z q = order z ((c, l)set cs. c ^ l)"
      by (simp add: order_smult q)
    also have " = (xset cs. order z (case x of (c, l)  c ^ l))"
      by (subst order_prod) (auto dest: coprime)
    also have " = ((c, l)set cs. l * order z c)"
      unfolding case_prod_unfold by (intro sum.cong refl, subst order_power) auto
    finally have "order z q = " .

    show "order z q  Suc k"
    proof (cases "c0 l0. (c0, l0)  set cs  poly c0 z = 0")
      case False
      have "order z q = ((c, l)set cs. l * order z c)" by fact
      also have "order z c = 0" if "(c, l)  set cs" for c l
        using False that by (auto simp: order_root)
      hence "((c, l)set cs. l * order z c) = 0"
        by (intro sum.neutral) auto
      finally show "order z q  Suc k" by simp
    next
      case True ―‹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. l * order z c)" by fact
      also have " = l0 * order z c0 + ((c, l)  set cs - {(c0, l0)}. l * order z c)"
        using cl0 by (subst sum.remove[of _ "(c0, l0)"]) auto
      also have "((c, l)  set cs - {(c0, l0)}. l * order z c) = 0"
      proof (intro sum.neutral ballI, goal_cases)
        case (1 cl)
        then obtain c l where [simp]: "cl = (c, l)" and cl: "(c, l)  set cs" "(c0, l0)  (c, l)"
          by (cases cl) auto
        from cl and cl0 and coprime[of c c0 l l0] have "coprime c c0"
          by auto
        with same_root_imp_not_coprime[of c z c0] and cl0 have "poly c z  0" by auto
        thus ?case by (auto simp: order_root)
      qed
      also have "square_free c0" using cl0 assms(1)
        by (auto simp: square_free_factorization_def)
      hence "rsquarefree c0" by (rule square_free_rsquarefree)
      with cl0 have "order z c0 = 1"
        by (auto simp: rsquarefree_def' order_root intro: antisym)
      finally have "order z q = l0" by simp

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

lemma proots_within_card_zero_iff:
  assumes "p  (0 :: 'a :: idom poly)"
  shows   "card (proots_within p A) = 0  (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 > Suc k) cs)"
  shows   "fps_nth (fps_of_poly p / fps_of_poly q)  O(λn. of_nat n ^ k * of_real R ^ n)"
proof (rule ratfps_nth_bigo_square_free_factorization[OF assms(1)])
  note sff = square_free_factorizationD[OF assms(1)]

  from sff(2)[of 0] have [simp]: "(0, x)  set cs" for x by auto
  from assms(1) have q: "q = smult b ((c, l)set cs. c ^ l)"
    unfolding square_free_factorization_def prod.case by blast 
  with q  0 have [simp]: "b  0" by auto

  show "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 > Suc k" for c l
  proof -
    from roots2 that have "card (proots_within c (sphere 0 (1 / R))) = 0"
      by (auto simp: proots_sphere_card_def list_all_def)
    with that show ?thesis by (subst (asm) proots_within_card_zero_iff) auto
  qed
qed fact+

definition ratfps_has_asymptotics where
  "ratfps_has_asymptotics q k R  q  0  R > 0 
     (let cs = snd (yun_factorization gcd q)
      in  list_all (λcl. proots_ball_card (fst cl) 0 (1 / R) = 0) cs 
          list_all (λcl. proots_sphere_card (fst cl) 0 (1 / R) = 0) (filter (λcl. snd cl > Suc k) cs))"

lemma ratfps_has_asymptotics_correct:
  assumes "ratfps_has_asymptotics q k R"
  shows   "fps_nth (fps_of_poly p / fps_of_poly q)  O(λn. of_nat n ^ k * of_real R ^ n)"
proof (rule ratfps_nth_bigo_square_free_factorization')
  show "square_free_factorization q (fst (yun_factorization gcd q), snd (yun_factorization gcd q))"
    by (rule yun_factorization) simp
qed (insert assms, auto simp: ratfps_has_asymptotics_def Let_def list_all_def)


value "map (fps_nth (fps_of_poly [:0, 1:] / fps_of_poly [:1, -1, -1 :: real:])) [0..<5]"




method ratfps_bigo = (rule ratfps_has_asymptotics_correct; eval)

lemma "fps_nth (fps_of_poly [:0, 1:] / fps_of_poly [:1, -1, -1 :: complex:]) 
         O(λn. of_nat n ^ 0 * complex_of_real 1.618034 ^ n)"
  by ratfps_bigo

lemma "fps_nth (fps_of_poly 1 / fps_of_poly [:1, -3, 3, -1 :: complex:]) 
         O(λn. of_nat n ^ 2 * complex_of_real 1 ^ n)"
  by ratfps_bigo

lemma "fps_nth (fps_of_poly f / fps_of_poly [:5, 4, 3, 2, 1 :: complex:]) 
         O(λn. of_nat n ^ 0 * complex_of_real 0.69202 ^ n)"
  by ratfps_bigo

end