Theory HOL-Real_Asymp.Multiseries_Expansion

(*
  File:    Multiseries_Expansion.thy
  Author:  Manuel Eberl, TU München

  Asymptotic bases and Multiseries expansions of real-valued functions.
  Contains automation to prove limits and asymptotic relationships between such functions.
*)
section ‹Multiseries expansions›
theory Multiseries_Expansion
imports
  "HOL-Library.BNF_Corec"
  "HOL-Library.Landau_Symbols"
  Lazy_Eval
  Inst_Existentials
  Eventuallize
begin

(* TODO Move *)
lemma real_powr_at_top: 
  assumes "(p::real) > 0"
  shows   "filterlim (λx. x powr p) at_top at_top"
proof (subst filterlim_cong[OF refl refl])
  show "LIM x at_top. exp (p * ln x) :> at_top"
    by (rule filterlim_compose[OF exp_at_top filterlim_tendsto_pos_mult_at_top[OF tendsto_const]])
       (simp_all add: ln_at_top assms)
  show "eventually (λx. x powr p = exp (p * ln x)) at_top"
    using eventually_gt_at_top[of 0] by eventually_elim (simp add: powr_def)
qed


subsection ‹Streams and lazy lists›

codatatype 'a msstream = MSSCons 'a "'a msstream"

primrec mssnth :: "'a msstream  nat  'a" where
  "mssnth xs 0 = (case xs of MSSCons x xs  x)"
| "mssnth xs (Suc n) = (case xs of MSSCons x xs  mssnth xs n)"

codatatype 'a msllist = MSLNil | MSLCons 'a "'a msllist"
  for map: mslmap

fun lcoeff where
  "lcoeff MSLNil n = 0"
| "lcoeff (MSLCons x xs) 0 = x"
| "lcoeff (MSLCons x xs) (Suc n) = lcoeff xs n"

primcorec msllist_of_msstream :: "'a msstream  'a msllist" where
  "msllist_of_msstream xs = (case xs of MSSCons x xs  MSLCons x (msllist_of_msstream xs))"
  
lemma msllist_of_msstream_MSSCons [simp]: 
  "msllist_of_msstream (MSSCons x xs) = MSLCons x (msllist_of_msstream xs)"
  by (subst msllist_of_msstream.code) simp

lemma lcoeff_llist_of_stream [simp]: "lcoeff (msllist_of_msstream xs) n = mssnth xs n"
  by (induction "msllist_of_msstream xs" n arbitrary: xs rule: lcoeff.induct;
      subst msllist_of_msstream.code) (auto split: msstream.splits)


subsection ‹Convergent power series›

subsubsection ‹Definition›

definition convergent_powser :: "real msllist  bool" where
  "convergent_powser cs  (R>0. x. abs x < R  summable (λn. lcoeff cs n * x ^ n))"
 
lemma convergent_powser_stl:
  assumes "convergent_powser (MSLCons c cs)"
  shows   "convergent_powser cs"
proof -
  from assms obtain R where "R > 0" "x. abs x < R  summable (λn. lcoeff (MSLCons c cs) n * x ^ n)"
    unfolding convergent_powser_def by blast
  hence "summable (λn. lcoeff cs n * x ^ n)" if "abs x < R" for x
    using that by (subst (asm) summable_powser_split_head [symmetric]) simp
  thus ?thesis using R > 0 by (auto simp: convergent_powser_def)
qed
  
  
definition powser :: "real msllist  real  real" where
  "powser cs x = suminf (λn. lcoeff cs n * x ^ n)"

lemma isCont_powser:
  assumes "convergent_powser cs"
  shows   "isCont (powser cs) 0"
proof -
  from assms obtain R where R: "R > 0" "x. abs x < R  summable (λn. lcoeff cs n * x ^ n)"
    unfolding convergent_powser_def by blast
  hence *: "summable (λn. lcoeff cs n * (R/2) ^ n)" by (intro R) simp_all
  from R > 0 show ?thesis unfolding powser_def
    by (intro isCont_powser[OF *]) simp_all
qed

lemma powser_MSLNil [simp]: "powser MSLNil = (λ_. 0)"
  by (simp add: fun_eq_iff powser_def)

lemma powser_MSLCons:
  assumes "convergent_powser (MSLCons c cs)"
  shows   "eventually (λx. powser (MSLCons c cs) x = x * powser cs x + c) (nhds 0)"
proof -
  from assms obtain R where R: "R > 0" "x. abs x < R  summable (λn. lcoeff (MSLCons c cs) n * x ^ n)"
    unfolding convergent_powser_def by blast
  from R have "powser (MSLCons c cs) x = x * powser cs x + c" if "abs x < R" for x
    using that unfolding powser_def by (subst powser_split_head) simp_all
  moreover have "eventually (λx. abs x < R) (nhds 0)"
    using R > 0 filterlim_ident[of "nhds (0::real)"]
    by (subst (asm) tendsto_iff) (simp add: dist_real_def)
  ultimately show ?thesis by (auto elim: eventually_mono)
qed

definition convergent_powser' :: "real msllist  (real  real)  bool" where
  "convergent_powser' cs f  (R>0. x. abs x < R  (λn. lcoeff cs n * x ^ n) sums f x)"

lemma convergent_powser'_imp_convergent_powser:
  "convergent_powser' cs f  convergent_powser cs" 
  unfolding convergent_powser_def convergent_powser'_def by (auto simp: sums_iff)

lemma convergent_powser'_eventually:
  assumes "convergent_powser' cs f"
  shows   "eventually (λx. powser cs x = f x) (nhds 0)"
proof -
  from assms obtain R where "R > 0" "x. abs x < R  (λn. lcoeff cs n * x ^ n) sums f x"
    unfolding convergent_powser'_def by blast
  hence "powser cs x = f x" if "abs x < R" for x
    using that by (simp add: powser_def sums_iff)
  moreover have "eventually (λx. abs x < R) (nhds 0)"
    using R > 0 filterlim_ident[of "nhds (0::real)"]
    by (subst (asm) tendsto_iff) (simp add: dist_real_def)
  ultimately show ?thesis by (auto elim: eventually_mono)
qed  


subsubsection ‹Geometric series›

primcorec mssalternate :: "'a  'a  'a msstream" where
  "mssalternate a b = MSSCons a (mssalternate b a)"

lemma case_mssalternate [simp]: 
  "(case mssalternate a b of MSSCons c d  f c d) = f a (mssalternate b a)"
  by (subst mssalternate.code) simp

lemma mssnth_mssalternate: "mssnth (mssalternate a b) n = (if even n then a else b)"
  by (induction n arbitrary: a b; subst mssalternate.code) simp_all
  
lemma convergent_powser'_geometric:
  "convergent_powser' (msllist_of_msstream (mssalternate 1 (-1))) (λx. inverse (1 + x))"
proof -
  have "(λn. lcoeff (msllist_of_msstream (mssalternate 1 (-1))) n * x ^ n) sums (inverse (1 + x))"
    if "abs x < 1" for x :: real
  proof -
    have "(λn. (-1) ^ n * x ^ n) sums (inverse (1 + x))"
      using that geometric_sums[of "-x"] by (simp add: field_simps power_mult_distrib [symmetric])
    also have "(λn. (-1) ^ n * x ^ n) =
                 (λn. lcoeff (msllist_of_msstream (mssalternate 1 (-1))) n * x ^ n)"
      by (auto simp add: mssnth_mssalternate fun_eq_iff)
    finally show ?thesis .
  qed
  thus ?thesis unfolding convergent_powser'_def by (force intro!: exI[of _ 1])
qed


subsubsection ‹Exponential series›

primcorec exp_series_stream_aux :: "real  real  real msstream" where
  "exp_series_stream_aux acc n = MSSCons acc (exp_series_stream_aux (acc / n) (n + 1))"
  
lemma mssnth_exp_series_stream_aux:
  "mssnth (exp_series_stream_aux (1 / fact n) (n + 1)) m = 1 / fact (n + m)"
proof (induction m arbitrary: n)
  case (0 n)
  thus ?case by (subst exp_series_stream_aux.code) simp
next
  case (Suc m n)
  from Suc.IH[of "n + 1"] show ?case
    by (subst exp_series_stream_aux.code) (simp add: algebra_simps)
qed

definition exp_series_stream :: "real msstream" where
  "exp_series_stream = exp_series_stream_aux 1 1"
  
lemma mssnth_exp_series_stream:
  "mssnth exp_series_stream n = 1 / fact n"
  unfolding exp_series_stream_def using mssnth_exp_series_stream_aux[of 0 n] by simp

lemma convergent_powser'_exp:
  "convergent_powser' (msllist_of_msstream exp_series_stream) exp"
proof -
  have "(λn. lcoeff (msllist_of_msstream exp_series_stream) n * x ^ n) sums exp x" for x :: real
    using exp_converges[of x] by (simp add: mssnth_exp_series_stream field_split_simps)
  thus ?thesis by (auto intro: exI[of _ "1::real"] simp: convergent_powser'_def)
qed
  

subsubsection ‹Logarithm series›

primcorec ln_series_stream_aux :: "bool  real  real msstream" where
  "ln_series_stream_aux b n = 
     MSSCons (if b then -inverse n else inverse n) (ln_series_stream_aux (¬b) (n+1))"

lemma mssnth_ln_series_stream_aux:
  "mssnth (ln_series_stream_aux b x) n = 
     (if b then - ((-1)^n) * inverse (x + real n) else ((-1)^n) * inverse (x + real n))"
  by (induction n arbitrary: b x; subst ln_series_stream_aux.code) simp_all

definition ln_series_stream :: "real msstream" where
  "ln_series_stream = MSSCons 0 (ln_series_stream_aux False 1)"
  
lemma mssnth_ln_series_stream:
  "mssnth ln_series_stream n = (-1) ^ Suc n / real n"
  unfolding ln_series_stream_def 
  by (cases n) (simp_all add: mssnth_ln_series_stream_aux field_simps)

lemma ln_series': 
  assumes "abs (x::real) < 1"
  shows   "(λn. - ((-x)^n) / of_nat n) sums ln (1 + x)"
proof -
  have "n1. norm (-((-x)^n)) / of_nat n  norm x ^ n / 1"
    by (intro exI[of _ "1::nat"] allI impI frac_le) (simp_all add: power_abs)
  hence "N. nN. norm (-((-x)^n) / of_nat n)  norm x ^ n" 
    by (intro exI[of _ 1]) simp_all
  moreover from assms have "summable (λn. norm x ^ n)" 
    by (intro summable_geometric) simp_all
  ultimately have *: "summable (λn. - ((-x)^n) / of_nat n)"
    by (rule summable_comparison_test)
  moreover from assms have "0 < 1 + x" "1 + x < 2" by simp_all
  from ln_series[OF this] 
    have "ln (1 + x) = (n. - ((-x) ^ Suc n) / real (Suc n))"
    by (simp_all add: power_minus' mult_ac)
  hence "ln (1 + x) = (n. - ((-x) ^ n / real n))"
    by (subst (asm) suminf_split_head[OF *]) simp_all
  ultimately show ?thesis by (simp add: sums_iff)
qed

lemma convergent_powser'_ln:
  "convergent_powser' (msllist_of_msstream ln_series_stream) (λx. ln (1 + x))"
proof -
  have "(λn. lcoeff (msllist_of_msstream ln_series_stream)n * x ^ n) sums ln (1 + x)"
    if "abs x < 1" for x
  proof -
    from that have "(λn. - ((- x) ^ n) / real n) sums ln (1 + x)" by (rule ln_series')
    also have "(λn. - ((- x) ^ n) / real n) =
                 (λn. lcoeff (msllist_of_msstream ln_series_stream) n * x ^ n)"
      by (auto simp: fun_eq_iff mssnth_ln_series_stream power_mult_distrib [symmetric])
    finally show ?thesis .
  qed
  thus ?thesis unfolding convergent_powser'_def by (force intro!: exI[of _ 1])
qed


subsubsection ‹Generalized binomial series›

primcorec gbinomial_series_aux :: "bool  real  real  real  real msllist" where
  "gbinomial_series_aux abort x n acc =
     (if abort  acc = 0 then MSLNil else 
        MSLCons acc (gbinomial_series_aux abort x (n + 1) ((x - n) / (n + 1) * acc)))"

lemma gbinomial_series_abort [simp]: "gbinomial_series_aux True x n 0 = MSLNil"
  by (subst gbinomial_series_aux.code) simp

definition gbinomial_series :: "bool  real  real msllist" where
  "gbinomial_series abort x = gbinomial_series_aux abort x 0 1"


(* TODO Move *)
lemma gbinomial_Suc_rec:
  "(x gchoose (Suc k)) = (x gchoose k) * ((x - of_nat k) / (of_nat k + 1))"
proof -
  have "((x - 1) + 1) gchoose (Suc k) = x * (x - 1 gchoose k) / (1 + of_nat k)"
    by (subst gbinomial_factors) simp
  also have "x * (x - 1 gchoose k) = (x - of_nat k) * (x gchoose k)"
    by (rule gbinomial_absorb_comp [symmetric])
  finally show ?thesis by (simp add: algebra_simps)
qed

lemma lcoeff_gbinomial_series_aux:
  "lcoeff (gbinomial_series_aux abort x m (x gchoose m)) n = (x gchoose (n + m))"
proof (induction n arbitrary: m)
  case 0
  show ?case by (subst gbinomial_series_aux.code) simp
next
  case (Suc n m)
  show ?case
  proof (cases "abort  (x gchoose m) = 0")
    case True
    with gbinomial_mult_fact[of m x] obtain k where "x = real k" "k < m"
      by auto
    hence "(x gchoose Suc (n + m)) = 0"
      using gbinomial_mult_fact[of "Suc (n + m)" x]
      by (simp add: gbinomial_altdef_of_nat)
    with True show ?thesis by (subst gbinomial_series_aux.code) simp
  next
    case False
    hence "lcoeff (gbinomial_series_aux abort x m (x gchoose m)) (Suc n) = 
             lcoeff (gbinomial_series_aux abort x (Suc m) ((x gchoose m) *
             ((x - real m) / (real m + 1)))) n"
      by (subst gbinomial_series_aux.code) (auto simp: algebra_simps)
    also have "((x gchoose m) * ((x - real m) / (real m + 1))) = x gchoose (Suc m)" 
      by (rule gbinomial_Suc_rec [symmetric])
    also have "lcoeff (gbinomial_series_aux abort x (Suc m) ) n = x gchoose (n + Suc m)"
      by (rule Suc.IH)
    finally show ?thesis by simp
  qed
qed

lemma lcoeff_gbinomial_series [simp]:
  "lcoeff (gbinomial_series abort x) n = (x gchoose n)"
  using lcoeff_gbinomial_series_aux[of abort x 0 n] by (simp add: gbinomial_series_def)

lemma gbinomial_ratio_limit:
  fixes a :: "'a :: real_normed_field"
  assumes "a  "
  shows "(λn. (a gchoose n) / (a gchoose Suc n))  -1"
proof (rule Lim_transform_eventually)
  let ?f = "λn. inverse (a / of_nat (Suc n) - of_nat n / of_nat (Suc n))"
  from eventually_gt_at_top[of "0::nat"]
    show "eventually (λn. ?f n = (a gchoose n) /(a gchoose Suc n)) sequentially"
  proof eventually_elim
    fix n :: nat assume n: "n > 0"
    then obtain q where q: "n = Suc q" by (cases n) blast
    let ?P = "i=0..<n. a - of_nat i"
    from n have "(a gchoose n) / (a gchoose Suc n) = (of_nat (Suc n) :: 'a) *
                   (?P / (i=0..n. a - of_nat i))"
      by (simp add: gbinomial_prod_rev atLeastLessThanSuc_atLeastAtMost)
    also from q have "(i=0..n. a - of_nat i) = ?P * (a - of_nat n)"
      by (simp add: prod.atLeast0_atMost_Suc atLeastLessThanSuc_atLeastAtMost)
    also have "?P /  = (?P / ?P) / (a - of_nat n)" by (rule divide_divide_eq_left[symmetric])
    also from assms have "?P / ?P = 1" by auto
    also have "of_nat (Suc n) * (1 / (a - of_nat n)) =
                   inverse (inverse (of_nat (Suc n)) * (a - of_nat n))" by (simp add: field_simps)
    also have "inverse (of_nat (Suc n)) * (a - of_nat n) =
                 a / of_nat (Suc n) - of_nat n / of_nat (Suc n)"
      by (simp add: field_simps del: of_nat_Suc)
    finally show "?f n = (a gchoose n) / (a gchoose Suc n)" by simp
  qed

  have "(λn. norm a / (of_nat (Suc n)))  0"
    unfolding divide_inverse
    by (intro tendsto_mult_right_zero LIMSEQ_inverse_real_of_nat)
  hence "(λn. a / of_nat (Suc n))  0"
    by (subst tendsto_norm_zero_iff[symmetric]) (simp add: norm_divide del: of_nat_Suc)
  hence "?f  inverse (0 - 1)"
    by (intro tendsto_inverse tendsto_diff LIMSEQ_n_over_Suc_n) simp_all
  thus "?f  -1" by simp
qed
  
lemma summable_gbinomial:
  fixes a z :: real
  assumes "norm z < 1"
  shows   "summable (λn. (a gchoose n) * z ^ n)"
proof (cases "z = 0  a  ")
  case False
  define r where "r = (norm z + 1) / 2"
  from assms have r: "r > norm z" "r < 1" by (simp_all add: r_def field_simps)
  from False have "abs z > 0" by auto
  from False have "a  " by (auto elim!: Nats_cases)
  hence *: "(λn. norm (z / ((a gchoose n) / (a gchoose (Suc n)))))  norm (z / (-1))"
    by (intro tendsto_norm tendsto_divide tendsto_const gbinomial_ratio_limit) simp_all
  hence "F x in at_top. norm (z / ((a gchoose x) / (a gchoose Suc x))) > 0"
    using assms False by (intro order_tendstoD) auto
  hence nz: "F x in at_top. (a gchoose x)  0" by eventually_elim auto
      
  have "F x in at_top. norm (z / ((a gchoose x) / (a gchoose Suc x))) < r"
    using assms r by (intro order_tendstoD(2)[OF *]) simp_all
  with nz have "F n in at_top. norm ((a gchoose (Suc n)) * z)  r * norm (a gchoose n)"
    by eventually_elim (simp add: field_simps abs_mult split: if_splits)
  hence "F n in at_top. norm (z ^ n) * norm ((a gchoose (Suc n)) * z) 
                           norm (z ^ n) * (r * norm (a gchoose n))"
    by eventually_elim (insert False, intro mult_left_mono, auto)
  hence "F n in at_top. norm ((a gchoose (Suc n)) * z ^ (Suc n))  
                           r * norm ((a gchoose n) * z ^ n)"
    by (simp add: abs_mult mult_ac)
  then obtain N where N: "n. n  N  norm ((a gchoose (Suc n)) * z ^ (Suc n)) 
                                          r * norm ((a gchoose n) * z ^ n)"
    unfolding eventually_at_top_linorder by blast
  show "summable (λn. (a gchoose n) * z ^ n)"
    by (intro summable_ratio_test[OF r(2), of N] N)
next
  case True
  thus ?thesis
  proof
    assume "a  "
    then obtain n where [simp]: "a = of_nat n" by (auto elim: Nats_cases)
    from eventually_gt_at_top[of n] 
      have "eventually (λn. (a gchoose n) * z ^ n = 0) at_top"
      by eventually_elim (simp add: binomial_gbinomial [symmetric])
    from summable_cong[OF this] show ?thesis by simp
  qed auto
qed

lemma gen_binomial_real:
  fixes z :: real
  assumes "norm z < 1"
  shows   "(λn. (a gchoose n) * z^n) sums (1 + z) powr a"
proof -
  define K where "K = 1 - (1 - norm z) / 2"
  from assms have K: "K > 0" "K < 1" "norm z < K"
     unfolding K_def by (auto simp: field_simps intro!: add_pos_nonneg)
  let ?f = "λn. a gchoose n" and ?f' = "diffs (λn. a gchoose n)"
  have summable_strong: "summable (λn. ?f n * z ^ n)" if "norm z < 1" for z using that
    by (intro summable_gbinomial)
  with K have summable: "summable (λn. ?f n * z ^ n)" if "norm z < K" for z using that by auto
  hence summable': "summable (λn. ?f' n * z ^ n)" if "norm z < K" for z using that
    by (intro termdiff_converges[of _ K]) simp_all

  define f f' where [abs_def]: "f z = (n. ?f n * z ^ n)" "f' z = (n. ?f' n * z ^ n)" for z
  {
    fix z :: real assume z: "norm z < K"
    from summable_mult2[OF summable'[OF z], of z]
      have summable1: "summable (λn. ?f' n * z ^ Suc n)" by (simp add: mult_ac)
    hence summable2: "summable (λn. of_nat n * ?f n * z^n)"
      unfolding diffs_def by (subst (asm) summable_Suc_iff)

    have "(1 + z) * f' z = (n. ?f' n * z^n) + (n. ?f' n * z^Suc n)"
      unfolding f_f'_def using summable' z by (simp add: algebra_simps suminf_mult)
    also have "(n. ?f' n * z^n) = (n. of_nat (Suc n) * ?f (Suc n) * z^n)"
      by (intro suminf_cong) (simp add: diffs_def)
    also have "(n. ?f' n * z^Suc n) = (n. of_nat n * ?f n * z ^ n)"
      using summable1 suminf_split_initial_segment[OF summable1] unfolding diffs_def
      by (subst suminf_split_head, subst (asm) summable_Suc_iff) simp_all
    also have "(n. of_nat (Suc n) * ?f (Suc n) * z^n) + (n. of_nat n * ?f n * z^n) =
                 (n. a * ?f n * z^n)"
      by (subst gbinomial_mult_1, subst suminf_add)
         (insert summable'[OF z] summable2,
          simp_all add: summable_powser_split_head algebra_simps diffs_def)
    also have " = a * f z" unfolding f_f'_def
      by (subst suminf_mult[symmetric]) (simp_all add: summable[OF z] mult_ac)
    finally have "a * f z = (1 + z) * f' z" by simp
  } note deriv = this

  have [derivative_intros]: "(f has_field_derivative f' z) (at z)" if "norm z < of_real K" for z
    unfolding f_f'_def using K that
    by (intro termdiffs_strong[of "?f" K z] summable_strong) simp_all
  have "f 0 = (n. if n = 0 then 1 else 0)" unfolding f_f'_def by (intro suminf_cong) simp
  also have " = 1" using sums_single[of 0 "λ_. 1::real"] unfolding sums_iff by simp
  finally have [simp]: "f 0 = 1" .

  have "f z * (1 + z) powr (-a) = f 0 * (1 + 0) powr (-a)"
  proof (rule DERIV_isconst3[where ?f = "λx. f x * (1 + x) powr (-a)"])
    fix z :: real assume z': "z  {-K<..<K}"
    hence z: "norm z < K" using K by auto
    with K have nz: "1 + z  0" by (auto dest!: minus_unique)
    from z K have "norm z < 1" by simp
    hence "((λz. f z * (1 + z) powr (-a)) has_field_derivative
              f' z * (1 + z) powr (-a) - a * f z * (1 + z) powr (-a-1)) (at z)" using z
      by (auto intro!: derivative_eq_intros)
    also from z have "a * f z = (1 + z) * f' z" by (rule deriv)
    also have "f' z * (1 + z) powr - a - (1 + z) * f' z * (1 + z) powr (- a - 1) = 0"
      using norm z < 1 by (auto simp add: field_simps powr_diff)
    finally show "((λz. f z * (1 + z) powr (-a)) has_field_derivative 0) (at z)" .
  qed (insert K, auto)
  also have "f 0 * (1 + 0) powr -a = 1" by simp
  finally have "f z = (1 + z) powr a" using K
    by (simp add: field_simps dist_real_def powr_minus)
  with summable K show ?thesis unfolding f_f'_def by (simp add: sums_iff)
qed

lemma convergent_powser'_gbinomial:
  "convergent_powser' (gbinomial_series abort p) (λx. (1 + x) powr p)"
proof -
  have "(λn. lcoeff (gbinomial_series abort p) n * x ^ n) sums (1 + x) powr p" if "abs x < 1" for x
    using that gen_binomial_real[of x p] by simp
  thus ?thesis unfolding convergent_powser'_def by (force intro!: exI[of _ 1])
qed

lemma convergent_powser'_gbinomial':
  "convergent_powser' (gbinomial_series abort (real n)) (λx. (1 + x) ^ n)"
proof -
  {
    fix x :: real
    have "(λk. if k  {0..n} then real (n choose k) * x ^ k else 0) sums (x + 1) ^ n"
      using sums_If_finite_set[of "{..n}" "λk. real (n choose k) * x ^ k"]
      by (subst binomial_ring) simp
    also have "(λk. if k  {0..n} then real (n choose k) * x ^ k else 0) = 
                 (λk. lcoeff (gbinomial_series abort (real n)) k * x ^ k)"
      by (simp add: fun_eq_iff binomial_gbinomial [symmetric])
    finally have " sums (1 + x) ^ n" by (simp add: add_ac)
  }
  thus ?thesis unfolding convergent_powser'_def
    by (auto intro!: exI[of _ 1])
qed


subsubsection ‹Sine and cosine›

primcorec sin_series_stream_aux :: "bool  real  real  real msstream" where
  "sin_series_stream_aux b acc m = 
    MSSCons (if b then -inverse acc else inverse acc)
          (sin_series_stream_aux (¬b) (acc * (m + 1) * (m + 2)) (m + 2))"

lemma mssnth_sin_series_stream_aux:
  "mssnth (sin_series_stream_aux b (fact m) m) n = 
     (if b then -1 else 1) * (-1) ^ n / (fact (2 * n + m))"
proof (induction n arbitrary: b m) 
  case (0 b m)
  thus ?case by (subst sin_series_stream_aux.code) (simp add: field_simps)
next
  case (Suc n b m)
  show ?case using Suc.IH[of "¬b" "m + 2"]
    by (subst sin_series_stream_aux.code) (auto simp: field_simps)
qed

definition sin_series_stream where
  "sin_series_stream = sin_series_stream_aux False 1 1"

lemma mssnth_sin_series_stream:
  "mssnth sin_series_stream n = (-1) ^ n / fact (2 * n + 1)"
  using mssnth_sin_series_stream_aux[of False 1 n] unfolding sin_series_stream_def by simp

definition cos_series_stream where
  "cos_series_stream = sin_series_stream_aux False 1 0"

lemma mssnth_cos_series_stream:
  "mssnth cos_series_stream n = (-1) ^ n / fact (2 * n)"
  using mssnth_sin_series_stream_aux[of False 0 n] unfolding cos_series_stream_def by simp

lemma summable_pre_sin: "summable (λn. mssnth sin_series_stream n * (x::real) ^ n)"
proof -
  have *: "summable (λn. 1 / fact n * abs x ^ n)"
    using exp_converges[of "abs x"] by (simp add: sums_iff field_simps)
  {
    fix n :: nat
    have "¦x¦ ^ n / fact (2 * n + 1)  ¦x¦ ^ n / fact n"
      by (intro divide_left_mono fact_mono) auto
    hence "norm (mssnth sin_series_stream n * x ^ n)  1 / fact n * abs x ^ n"
      by (simp add: mssnth_sin_series_stream abs_mult power_abs field_simps)
  }
  thus ?thesis
    by (intro summable_comparison_test[OF _ *]) (auto intro!: exI[of _ 0])
qed

lemma summable_pre_cos: "summable (λn. mssnth cos_series_stream n * (x::real) ^ n)"
proof -
  have *: "summable (λn. 1 / fact n * abs x ^ n)"
    using exp_converges[of "abs x"] by (simp add: sums_iff field_simps)
  {
    fix n :: nat
    have "¦x¦ ^ n / fact (2 * n)  ¦x¦ ^ n / fact n"
      by (intro divide_left_mono fact_mono) auto
    hence "norm (mssnth cos_series_stream n * x ^ n)  1 / fact n * abs x ^ n"
      by (simp add: mssnth_cos_series_stream abs_mult power_abs field_simps)
  }
  thus ?thesis
    by (intro summable_comparison_test[OF _ *]) (auto intro!: exI[of _ 0])
qed

lemma cos_conv_pre_cos:
  "cos x = powser (msllist_of_msstream cos_series_stream) (x ^ 2)"
proof -
  have "(λn. cos_coeff (2 * n) * x ^ (2 * n)) sums cos x"
    using cos_converges[of x]
    by (subst sums_mono_reindex[of "λn. 2 * n"])
       (auto simp: strict_mono_def cos_coeff_def elim!: evenE)
  also have "(λn. cos_coeff (2 * n) * x ^ (2 * n)) =
               (λn. mssnth cos_series_stream n * (x ^ 2) ^ n)"
    by (simp add: fun_eq_iff mssnth_cos_series_stream cos_coeff_def power_mult)
  finally have sums: "(λn. mssnth cos_series_stream n * x2 ^ n) sums cos x" .
  thus ?thesis by (simp add: powser_def sums_iff)
qed

lemma sin_conv_pre_sin:
  "sin x = x * powser (msllist_of_msstream sin_series_stream) (x ^ 2)"
proof -
  have "(λn. sin_coeff (2 * n + 1) * x ^ (2 * n + 1)) sums sin x"
    using sin_converges[of x]
    by (subst sums_mono_reindex[of "λn. 2 * n + 1"])
       (auto simp: strict_mono_def sin_coeff_def elim!: oddE)
  also have "(λn. sin_coeff (2 * n + 1) * x ^ (2 * n + 1)) =
               (λn. x * mssnth sin_series_stream n * (x ^ 2) ^ n)"
    by (simp add: fun_eq_iff mssnth_sin_series_stream sin_coeff_def power_mult [symmetric] algebra_simps)
  finally have sums: "(λn. x * mssnth sin_series_stream n * x2 ^ n) sums sin x" .
  thus ?thesis using summable_pre_sin[of "x^2"] 
    by (simp add: powser_def sums_iff suminf_mult [symmetric] mult.assoc)
qed

lemma convergent_powser_sin_series_stream:
  "convergent_powser (msllist_of_msstream sin_series_stream)"
  (is "convergent_powser ?cs")
proof -
  show ?thesis using summable_pre_sin unfolding convergent_powser_def
    by (intro exI[of _ 1]) auto
qed

lemma convergent_powser_cos_series_stream:
  "convergent_powser (msllist_of_msstream cos_series_stream)"
  (is "convergent_powser ?cs")
proof -
  show ?thesis using summable_pre_cos unfolding convergent_powser_def
    by (intro exI[of _ 1]) auto
qed
  
  
subsubsection ‹Arc tangent›

definition arctan_coeffs :: "nat  'a :: {real_normed_div_algebra, banach}" where
  "arctan_coeffs n = (if odd n then (-1) ^ (n div 2) / of_nat n else 0)"

lemma arctan_converges:
  assumes "norm x < 1"
  shows   "(λn. arctan_coeffs n * x ^ n) sums arctan x"
proof -
  have A: "(λn. diffs arctan_coeffs n * x ^ n) sums (1 / (1 + x^2))" if "norm x < 1" for x :: real
  proof -
    let ?f = "λn. (if odd n then 0 else (-1) ^ (n div 2)) * x ^ n"
    from that have "norm x ^ 2 < 1 ^ 2" by (intro power_strict_mono) simp_all
    hence "(λn. (-(x^2))^n) sums (1 / (1 - (-(x^2))))" by (intro geometric_sums) simp_all
    also have "1 - (-(x^2)) = 1 + x^2" by simp
    also have "(λn. (-(x^2))^n) = (λn. ?f (2*n))" by (auto simp: fun_eq_iff power_minus' power_mult)
    also have "range ((*) (2::nat)) = {n. even n}"
      by (auto elim!: evenE)
    hence "(λn. ?f (2*n)) sums (1 / (1 + x^2))  ?f sums (1 / (1 + x^2))"
      by (intro sums_mono_reindex) (simp_all add: strict_mono_Suc_iff)
    also have "?f = (λn. diffs arctan_coeffs n * x ^ n)"
      by (simp add: fun_eq_iff arctan_coeffs_def diffs_def)
    finally show "(λn. diffs arctan_coeffs n * x ^ n) sums (1 / (1 + x^2))" .
  qed
  
  have B: "summable (λn. arctan_coeffs n * x ^ n)" if "norm x < 1" for x :: real
  proof (rule summable_comparison_test)
    show "N. nN. norm (arctan_coeffs n * x ^ n)  1 * norm x ^ n"
      unfolding norm_mult norm_power
      by (intro exI[of _ "0::nat"] allI impI mult_right_mono)
         (simp_all add: arctan_coeffs_def field_split_simps)
    from that show "summable (λn. 1 * norm x ^ n)"
      by (intro summable_mult summable_geometric) simp_all
  qed
  
  define F :: "real  real" where "F = (λx. n. arctan_coeffs n * x ^ n)"
  have [derivative_intros]:
    "(F has_field_derivative (1 / (1 + x^2))) (at x)" if "norm x < 1" for x :: real
  proof -
    define K :: real where "K = (1 + norm x) / 2"
    from that have K: "norm x < K" "K < 1" by (simp_all add: K_def field_simps)
    have "(F has_field_derivative (n. diffs arctan_coeffs n * x ^ n)) (at x)"
      unfolding F_def using K by (intro termdiffs_strong[where K = K] B) simp_all
    also from A[OF that] have "(n. diffs arctan_coeffs n * x ^ n) = 1 / (1 + x^2)"
      by (simp add: sums_iff)
    finally show ?thesis .
  qed
  from assms have "arctan x - F x = arctan 0 - F 0"
    by (intro DERIV_isconst3[of "-1" 1 _ _ "λx. arctan x - F x"])
       (auto intro!: derivative_eq_intros simp: field_simps)
  hence "F x = arctan x" by (simp add: F_def arctan_coeffs_def)
  with B[OF assms] show ?thesis by (simp add: sums_iff F_def)
qed

primcorec arctan_series_stream_aux :: "bool  real  real msstream" where
  "arctan_series_stream_aux b n = 
     MSSCons (if b then -inverse n else inverse n) (arctan_series_stream_aux (¬b) (n + 2))"

lemma mssnth_arctan_series_stream_aux: 
  "mssnth (arctan_series_stream_aux b n) m = (-1) ^ (if b then Suc m else m) / (2*m + n)"
  by (induction m arbitrary: b n; subst arctan_series_stream_aux.code)
     (auto simp: field_simps split: if_splits)

definition arctan_series_stream where
  "arctan_series_stream = arctan_series_stream_aux False 1"

lemma mssnth_arctan_series_stream:
  "mssnth arctan_series_stream n = (-1) ^ n / (2 * n + 1)"
  by (simp add: arctan_series_stream_def mssnth_arctan_series_stream_aux)

lemma summable_pre_arctan:
  assumes "norm x < 1"
  shows   "summable (λn. mssnth arctan_series_stream n * x ^ n)" (is "summable ?f")
proof (rule summable_comparison_test)
  show "summable (λn. abs x ^ n)" using assms by (intro summable_geometric) auto
  show "N. nN. norm (?f n)  abs x ^ n"
  proof (intro exI[of _ 0] allI impI)
    fix n :: nat
    have "norm (?f n) = ¦x¦ ^ n / (2 * real n + 1)"
      by (simp add: mssnth_arctan_series_stream abs_mult power_abs)
    also have "  ¦x¦ ^ n / 1" by (intro divide_left_mono) auto
    finally show "norm (?f n)  abs x ^ n" by simp
  qed
qed

lemma arctan_conv_pre_arctan:
  assumes "norm x < 1"
  shows   "arctan x = x * powser (msllist_of_msstream arctan_series_stream) (x ^ 2)"
proof -
  from assms have "norm x * norm x < 1 * 1" by (intro mult_strict_mono) auto
  hence norm_less: "norm (x ^ 2) < 1" by (simp add: power2_eq_square)
  have "(λn. arctan_coeffs n * x ^ n) sums arctan x"
    by (intro arctan_converges assms)
  also have "?this  (λn. arctan_coeffs (2 * n + 1) * x ^ (2 * n + 1)) sums arctan x"
    by (intro sums_mono_reindex [symmetric])
       (auto simp: arctan_coeffs_def strict_mono_def elim!: oddE)
  also have "(λn. arctan_coeffs (2 * n + 1) * x ^ (2 * n + 1)) =
               (λn. x * mssnth arctan_series_stream n * (x ^ 2) ^ n)"
    by (simp add: fun_eq_iff mssnth_arctan_series_stream 
                  power_mult [symmetric] arctan_coeffs_def mult_ac)
  finally have "(λn. x * mssnth arctan_series_stream n * x2 ^ n) sums arctan x" .
  thus ?thesis using summable_pre_arctan[OF norm_less] assms
    by (simp add: powser_def sums_iff suminf_mult [symmetric] mult.assoc)
qed

lemma convergent_powser_arctan: 
  "convergent_powser (msllist_of_msstream arctan_series_stream)"
  unfolding convergent_powser_def
  using summable_pre_arctan by (intro exI[of _ 1]) auto

lemma arctan_inverse_pos: "x > 0  arctan x = pi / 2 - arctan (inverse x)"
  using arctan_inverse[of x] by (simp add: field_simps)
    
lemma arctan_inverse_neg: "x < 0  arctan x = -pi / 2 - arctan (inverse x)"
  using arctan_inverse[of x] by (simp add: field_simps)



subsection ‹Multiseries›

subsubsection ‹Asymptotic bases›

type_synonym basis = "(real  real) list"

lemma transp_ln_smallo_ln: "transp (λf g. (λx::real. ln (g x))  o(λx. ln (f x)))"
  by (rule transpI, erule landau_o.small.trans)

definition basis_wf :: "basis  bool" where
  "basis_wf basis  (fset basis. filterlim f at_top at_top)  
                      sorted_wrt (λf g. (λx. ln (g x))  o(λx. ln (f x))) basis"

lemma basis_wf_Nil [simp]: "basis_wf []"
  by (simp add: basis_wf_def)

lemma basis_wf_Cons: 
  "basis_wf (f # basis)  
     filterlim f at_top at_top  (gset basis. (λx. ln (g x))  o(λx. ln (f x)))  basis_wf basis"
  unfolding basis_wf_def by auto

(* TODO: Move *)
lemma sorted_wrt_snoc:
  assumes trans: "transp P" and "sorted_wrt P xs" "P (last xs) y"
  shows   "sorted_wrt P (xs @ [y])"
proof (subst sorted_wrt_append, intro conjI ballI)
  fix x y' assume x: "x  set xs" and y': "y'  set [y]"
  from y' have [simp]: "y' = y" by simp
  show "P x y'"
  proof (cases "x = last xs")
    case False
    from x have eq: "xs = butlast xs @ [last xs]"
      by (subst append_butlast_last_id) auto
    from x and False have x': "x  set (butlast xs)" by (subst (asm) eq) auto
    have "sorted_wrt P xs" by fact
    also note eq
    finally have "P x (last xs)" using x'
      by (subst (asm) sorted_wrt_append) auto
    with P (last xs) y have "P x y" using transpD[OF trans] by blast
    thus ?thesis by simp
  qed (insert assms, auto)
qed (insert assms, auto)  

lemma basis_wf_snoc:
  assumes "bs  []"
  assumes "basis_wf bs" "filterlim b at_top at_top"
  assumes "(λx. ln (b x))  o(λx. ln (last bs x))"
  shows   "basis_wf (bs @ [b])"
proof -
  have "transp (λf g. (λx. ln (g x))  o(λx. ln (f x)))"
    by (auto simp: transp_def intro: landau_o.small.trans)
  thus ?thesis using assms
    by (auto simp add: basis_wf_def sorted_wrt_snoc[OF transp_ln_smallo_ln])
qed

lemma basis_wf_singleton: "basis_wf [b]  filterlim b at_top at_top"
  by (simp add: basis_wf_Cons)

lemma basis_wf_many: "basis_wf (b # b' # bs) 
  filterlim b at_top at_top  (λx. ln (b' x))  o(λx. ln (b x))  basis_wf (b' # bs)"
  unfolding basis_wf_def by (subst sorted_wrt2[OF transp_ln_smallo_ln]) auto


subsubsection ‹Monomials›

type_synonym monom = "real × real list"

definition eval_monom :: "monom  basis  real  real" where
  "eval_monom = (λ(c,es) basis x. c * prod_list (map (λ(b,e). b x powr e) (zip basis es)))"
  
lemma eval_monom_Nil1 [simp]: "eval_monom (c, []) basis = (λ_. c)"
  by (simp add: eval_monom_def split: prod.split)

lemma eval_monom_Nil2 [simp]: "eval_monom monom [] = (λ_. fst monom)"
  by (simp add: eval_monom_def split: prod.split)

lemma eval_monom_Cons: 
  "eval_monom (c, e # es) (b # basis) = (λx. eval_monom (c, es) basis x * b x powr e)"
  by (simp add: eval_monom_def fun_eq_iff split: prod.split)

definition inverse_monom :: "monom  monom" where
  "inverse_monom monom = (case monom of (c, es)  (inverse c, map uminus es))"

lemma length_snd_inverse_monom [simp]: 
  "length (snd (inverse_monom monom)) = length (snd monom)"
  by (simp add: inverse_monom_def split: prod.split)

lemma eval_monom_pos:
  assumes "basis_wf basis" "fst monom > 0"
  shows   "eventually (λx. eval_monom monom basis x > 0) at_top"
proof (cases monom)
  case (Pair c es)
  with assms have "c > 0" by simp
  with assms(1) show ?thesis unfolding Pair
  proof (induction es arbitrary: basis)
    case (Cons e es)
    note A = this
    thus ?case
    proof (cases basis)
      case (Cons b basis')
      with A(1)[of basis'] A(2,3) 
        have A: "F x in at_top. eval_monom (c, es) basis' x > 0" 
         and B: "eventually (λx. b x > 0) at_top"
        by (simp_all add: basis_wf_Cons filterlim_at_top_dense)
      thus ?thesis by eventually_elim (simp add: Cons eval_monom_Cons)
    qed simp
  qed simp
qed

lemma eval_monom_uminus: "eval_monom (-c, es) basis x = -eval_monom (c, es) basis x"
  by (simp add: eval_monom_def)

lemma eval_monom_neg:
  assumes "basis_wf basis" "fst monom < 0"
  shows   "eventually (λx. eval_monom monom basis x < 0) at_top"    
proof -
  from assms have "eventually (λx. eval_monom (-fst monom, snd monom) basis x > 0) at_top"
    by (intro eval_monom_pos) simp_all
  thus ?thesis by (simp add: eval_monom_uminus)
qed

lemma eval_monom_nonzero:
  assumes "basis_wf basis" "fst monom  0"
  shows   "eventually (λx. eval_monom monom basis x  0) at_top"
proof (cases "fst monom" "0 :: real" rule: linorder_cases)
  case greater
  with assms have "eventually (λx. eval_monom monom basis x > 0) at_top"
    by (simp add: eval_monom_pos)
  thus ?thesis by eventually_elim simp
next
  case less
  with assms have "eventually (λx. eval_monom (-fst monom, snd monom) basis x > 0) at_top"
    by (simp add: eval_monom_pos)
  thus ?thesis by eventually_elim (simp add: eval_monom_uminus)
qed (insert assms, simp_all)


subsubsection ‹Typeclass for multiseries›

class multiseries = plus + minus + times + uminus + inverse + 
  fixes is_expansion :: "'a  basis  bool"
    and expansion_level :: "'a itself  nat"
    and eval :: "'a  real  real"
    and zero_expansion :: 'a
    and const_expansion :: "real  'a"
    and powr_expansion :: "bool  'a  real  'a"
    and power_expansion :: "bool  'a  nat  'a"
    and trimmed :: "'a  bool"
    and dominant_term :: "'a  monom"

  assumes is_expansion_length:
    "is_expansion F basis  length basis = expansion_level (TYPE('a))"
  assumes is_expansion_zero:
    "basis_wf basis  length basis = expansion_level (TYPE('a))  
       is_expansion zero_expansion basis"
  assumes is_expansion_const:
    "basis_wf basis  length basis = expansion_level (TYPE('a))  
       is_expansion (const_expansion c) basis"
  assumes is_expansion_uminus:
    "basis_wf basis  is_expansion F basis  is_expansion (-F) basis"
  assumes is_expansion_add: 
    "basis_wf basis  is_expansion F basis  is_expansion G basis  
       is_expansion (F + G) basis"
  assumes is_expansion_minus: 
    "basis_wf basis  is_expansion F basis  is_expansion G basis  
       is_expansion (F - G) basis"
  assumes is_expansion_mult: 
    "basis_wf basis  is_expansion F basis  is_expansion G basis  
       is_expansion (F * G) basis"
  assumes is_expansion_inverse:
    "basis_wf basis  trimmed F  is_expansion F basis  
       is_expansion (inverse F) basis"
  assumes is_expansion_divide:
    "basis_wf basis  trimmed G  is_expansion F basis  is_expansion G basis  
       is_expansion (F / G) basis"
  assumes is_expansion_powr:
    "basis_wf basis  trimmed F  fst (dominant_term F) > 0  is_expansion F basis 
       is_expansion (powr_expansion abort F p) basis"
  assumes is_expansion_power:
    "basis_wf basis  trimmed F  is_expansion F basis 
       is_expansion (power_expansion abort F n) basis"
  
  assumes is_expansion_imp_smallo:
    "basis_wf basis  is_expansion F basis  filterlim b at_top at_top 
       (gset basis. (λx. ln (g x))  o(λx. ln (b x)))  e > 0  eval F  o(λx. b x powr e)"
  assumes is_expansion_imp_smallomega:
    "basis_wf basis  is_expansion F basis  filterlim b at_top at_top  trimmed F  
       (gset basis. (λx. ln (g x))  o(λx. ln (b x)))  e < 0  eval F  ω(λx. b x powr e)"
  assumes trimmed_imp_eventually_sgn:
    "basis_wf basis  is_expansion F basis  trimmed F 
       eventually (λx. sgn (eval F x) = sgn (fst (dominant_term F))) at_top"
  assumes trimmed_imp_eventually_nz: 
    "basis_wf basis  is_expansion F basis  trimmed F  
       eventually (λx. eval F x  0) at_top"
  assumes trimmed_imp_dominant_term_nz: "trimmed F  fst (dominant_term F)  0"
  
  assumes dominant_term: 
    "basis_wf basis  is_expansion F basis  trimmed F 
       eval F ∼[at_top] eval_monom (dominant_term F) basis"
  assumes dominant_term_bigo:
    "basis_wf basis  is_expansion F basis 
       eval F  O(eval_monom (1, snd (dominant_term F)) basis)"
  assumes length_dominant_term:
    "length (snd (dominant_term F)) = expansion_level TYPE('a)"
  assumes fst_dominant_term_uminus [simp]: "fst (dominant_term (- F)) = -fst (dominant_term F)"
  assumes trimmed_uminus_iff [simp]: "trimmed (-F)  trimmed F"
  
  assumes add_zero_expansion_left [simp]: "zero_expansion + F = F"
  assumes add_zero_expansion_right [simp]: "F + zero_expansion = F"
  
  assumes eval_zero [simp]: "eval zero_expansion x = 0"
  assumes eval_const [simp]: "eval (const_expansion c) x = c"
  assumes eval_uminus [simp]: "eval (-F) = (λx. -eval F x)"
  assumes eval_plus [simp]: "eval (F + G) = (λx. eval F x + eval G x)"
  assumes eval_minus [simp]: "eval (F - G) = (λx. eval F x - eval G x)"
  assumes eval_times [simp]: "eval (F * G) = (λx. eval F x * eval G x)"
  assumes eval_inverse [simp]: "eval (inverse F) = (λx. inverse (eval F x))"
  assumes eval_divide [simp]: "eval (F / G) = (λx. eval F x / eval G x)"
  assumes eval_powr [simp]: "eval (powr_expansion abort F p) = (λx. eval F x powr p)"
  assumes eval_power [simp]: "eval (power_expansion abort F n) = (λx. eval F x ^ n)"
  assumes minus_eq_plus_uminus: "F - G = F + (-G)"
  assumes times_const_expansion_1: "const_expansion 1 * F = F"
  assumes trimmed_const_expansion: "trimmed (const_expansion c)  c  0"
begin

definition trimmed_pos where
  "trimmed_pos F  trimmed F  fst (dominant_term F) > 0"

definition trimmed_neg where
  "trimmed_neg F  trimmed F  fst (dominant_term F) < 0"

lemma trimmed_pos_uminus: "trimmed_neg F  trimmed_pos (-F)"
  by (simp add: trimmed_neg_def trimmed_pos_def)

lemma trimmed_pos_imp_trimmed: "trimmed_pos x  trimmed x"
  by (simp add: trimmed_pos_def)

lemma trimmed_neg_imp_trimmed: "trimmed_neg x  trimmed x"
  by (simp add: trimmed_neg_def)

end


subsubsection ‹Zero-rank expansions›

instantiation real :: multiseries
begin

inductive is_expansion_real :: "real  basis  bool" where
  "is_expansion_real c []"
  
definition expansion_level_real :: "real itself  nat" where
  expansion_level_real_def [simp]: "expansion_level_real _ = 0"

definition eval_real :: "real  real  real" where
  eval_real_def [simp]: "eval_real c = (λ_. c)"

definition zero_expansion_real :: "real" where
  zero_expansion_real_def [simp]: "zero_expansion_real = 0"
  
definition const_expansion_real :: "real  real" where
  const_expansion_real_def [simp]: "const_expansion_real x = x"

definition powr_expansion_real :: "bool  real  real  real" where
  powr_expansion_real_def [simp]: "powr_expansion_real _ x p = x powr p"

definition power_expansion_real :: "bool  real  nat  real" where
  power_expansion_real_def [simp]: "power_expansion_real _ x n = x ^ n"

definition trimmed_real :: "real  bool" where
  trimmed_real_def [simp]: "trimmed_real x  x  0"

definition dominant_term_real :: "real  monom" where
  dominant_term_real_def [simp]: "dominant_term_real c = (c, [])"

instance proof
  fix basis :: basis and b :: "real  real" and e F :: real
  assume *: "basis_wf basis" "filterlim b at_top at_top" "is_expansion F basis" "0 < e"
  have "(λx. b x powr e)  ω(λ_. 1)"
    by (intro smallomegaI_filterlim_at_infinity filterlim_at_top_imp_at_infinity) 
       (auto intro!: filterlim_compose[OF real_powr_at_top] * )
  with * show "eval F  o(λx. b x powr e)"
    by (cases "F = 0") (auto elim!: is_expansion_real.cases simp: smallomega_iff_smallo)
next
  fix basis :: basis and b :: "real  real" and e F :: real
  assume *: "basis_wf basis" "filterlim b at_top at_top" "is_expansion F basis" 
            "e < 0" "trimmed F"
  from * have pos: "eventually (λx. b x > 0) at_top" by (simp add: filterlim_at_top_dense)
  have "(λx. b x powr -e)  ω(λ_. 1)"
    by (intro smallomegaI_filterlim_at_infinity filterlim_at_top_imp_at_infinity) 
       (auto intro!: filterlim_compose[OF real_powr_at_top] * )
  with pos have "(λx. b x powr e)  o(λ_. 1)" unfolding powr_minus
    by (subst (asm) landau_omega.small.inverse_eq2)
       (auto elim: eventually_mono simp: smallomega_iff_smallo)
  with * show "eval F  ω(λx. b x powr e)"
    by (auto elim!: is_expansion_real.cases simp: smallomega_iff_smallo)
qed (auto intro!: is_expansion_real.intros elim!: is_expansion_real.cases)

end

lemma eval_real: "eval (c :: real) x = c" by simp


subsubsection ‹Operations on multiseries›

lemma ms_aux_cases [case_names MSLNil MSLCons]:
  fixes xs :: "('a × real) msllist"
  obtains "xs = MSLNil" | c e xs' where "xs = MSLCons (c, e) xs'"
proof (cases xs)
  case (MSLCons x xs')
  with that(2)[of "fst x" "snd x" xs'] show ?thesis by auto
qed auto

definition ms_aux_hd_exp :: "('a × real) msllist  real option" where
  "ms_aux_hd_exp xs = (case xs of MSLNil  None | MSLCons (_, e) _  Some e)"

primrec ms_exp_gt :: "real  real option  bool" where
  "ms_exp_gt x None = True"
| "ms_exp_gt x (Some y)  x > y"

lemma ms_aux_hd_exp_MSLNil [simp]: "ms_aux_hd_exp MSLNil = None"
  by (simp add: ms_aux_hd_exp_def split: prod.split)
  
lemma ms_aux_hd_exp_MSLCons [simp]: "ms_aux_hd_exp (MSLCons x xs) = Some (snd x)"
  by (simp add: ms_aux_hd_exp_def split: prod.split)

coinductive is_expansion_aux :: 
    "('a :: multiseries × real) msllist  (real  real)  basis  bool" where
  is_expansion_MSLNil: 
    "eventually (λx. f x = 0) at_top  length basis = Suc (expansion_level TYPE('a)) 
       is_expansion_aux MSLNil f basis"
| is_expansion_MSLCons: 
    "is_expansion_aux xs (λx. f x - eval C x * b x powr e) (b # basis) 
       is_expansion C basis 
       (e'. e' > e  f  o(λx. b x powr e'))  ms_exp_gt e (ms_aux_hd_exp xs) 
       is_expansion_aux (MSLCons (C, e) xs) f (b # basis)"

inductive_cases is_expansion_aux_MSLCons: "is_expansion_aux (MSLCons (c, e) xs) f basis"
 
lemma is_expansion_aux_basis_nonempty: "is_expansion_aux F f basis  basis  []"
  by (erule is_expansion_aux.cases) auto

lemma is_expansion_aux_expansion_level:
  assumes "is_expansion_aux (G :: ('a::multiseries × real) msllist) g basis"
  shows   "length basis = Suc (expansion_level TYPE('a))"
  using assms by (cases rule: is_expansion_aux.cases) (auto dest: is_expansion_length)  

lemma is_expansion_aux_imp_smallo:
  assumes "is_expansion_aux xs f basis" "ms_exp_gt p (ms_aux_hd_exp xs)" 
  shows   "f  o(λx. hd basis x powr p)"
  using assms
proof (cases rule: is_expansion_aux.cases)
  case is_expansion_MSLNil
  show ?thesis by (simp add: landau_o.small.in_cong[OF is_expansion_MSLNil(2)])
next
  case (is_expansion_MSLCons xs C b e basis)
  with assms have "f  o(λx. b x powr p)"
    by (intro is_expansion_MSLCons) (simp_all add: ms_aux_hd_exp_def)
  thus ?thesis by (simp add: is_expansion_MSLCons)
qed

lemma is_expansion_aux_imp_smallo':
  assumes "is_expansion_aux xs f basis"
  obtains e where "f  o(λx. hd basis x powr e)"
proof -
  define e where "e = (case ms_aux_hd_exp xs of None  0 | Some e  e + 1)"
  have "ms_exp_gt e (ms_aux_hd_exp xs)"
    by (auto simp add: e_def ms_aux_hd_exp_def split: msllist.splits)
  from assms this have "f  o(λx. hd basis x powr e)" by (rule is_expansion_aux_imp_smallo)
  from that[OF this] show ?thesis .
qed

lemma is_expansion_aux_imp_smallo'':
  assumes "is_expansion_aux xs f basis" "ms_exp_gt e' (ms_aux_hd_exp xs)"
  obtains e where "e < e'" "f  o(λx. hd basis x powr e)"
proof -
  define e where "e = (case ms_aux_hd_exp xs of None  e' - 1 | Some e  (e' + e) / 2)"
  from assms have "ms_exp_gt e (ms_aux_hd_exp xs)" "e < e'"
    by (cases xs; simp add: e_def field_simps)+
  from assms(1) this(1) have "f  o(λx. hd basis x powr e)" by (rule is_expansion_aux_imp_smallo)
  from that[OF e < e' this] show ?thesis .
qed

definition trimmed_ms_aux :: "('a :: multiseries × real) msllist  bool" where
  "trimmed_ms_aux xs = (case xs of MSLNil  False | MSLCons (C, _) _  trimmed C)"
 
lemma not_trimmed_ms_aux_MSLNil [simp]: "¬trimmed_ms_aux MSLNil"
  by (simp add: trimmed_ms_aux_def)

lemma trimmed_ms_aux_MSLCons: "trimmed_ms_aux (MSLCons x xs) = trimmed (fst x)"
  by (simp add: trimmed_ms_aux_def split: prod.split)

lemma trimmed_ms_aux_imp_nz:
  assumes "basis_wf basis" "is_expansion_aux xs f basis" "trimmed_ms_aux xs"
  shows   "eventually (λx. f x  0) at_top"
proof (cases xs rule: ms_aux_cases)
  case (MSLCons C e xs')
  from assms this obtain b basis' where *: "basis = b # basis'"
    "is_expansion_aux xs' (λx. f x - eval C x * b x powr e) (b # basis')" 
    "ms_exp_gt