Theory Pi_Transcendental

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

   A proof of the transcendence of pi
*)
section ‹The Transcendence of $\pi$›
theory Pi_Transcendental
imports
  "E_Transcendental.E_Transcendental"
  "Symmetric_Polynomials.Symmetric_Polynomials"
  "HOL-Real_Asymp.Real_Asymp" 
begin

lemma ring_homomorphism_to_poly [intro]: "ring_homomorphism (λi. [:i:])"
  by standard auto

lemma (in ring_closed) coeff_power_closed:
  "(⋀m. coeff p m ∈ A) ⟹ coeff (p ^ n) m ∈ A"
  by (induction n arbitrary: m)
     (auto simp: mpoly_coeff_1 coeff_mpoly_times intro!: prod_fun_closed)

lemma (in ring_closed) coeff_prod_closed:
  "(⋀x m. x ∈ X ⟹ coeff (f x) m ∈ A) ⟹ coeff (prod f X) m ∈ A"
  by (induction X arbitrary: m rule: infinite_finite_induct)
     (auto simp: mpoly_coeff_1 coeff_mpoly_times intro!: prod_fun_closed)

lemma map_of_rat_of_int_poly [simp]: "map_poly of_rat (of_int_poly p) = of_int_poly p"
  by (intro poly_eqI) (auto simp: coeff_map_poly)

text ‹
  Given a polynomial with rational coefficients, we can obtain an integer polynomial that
  differs from it only by a nonzero constant by clearing the denominators.
›
lemma ratpoly_to_intpoly:
  assumes "∀i. poly.coeff p i ∈ ℚ"
  obtains q c where "c ≠ 0" "p = Polynomial.smult (inverse (of_nat c)) (of_int_poly q)"
proof (cases "p = 0")
  case True
  with that[of 1 0] show ?thesis by auto
next
  case False
  from assms obtain p' where p': "p = map_poly of_rat p'"
    using ratpolyE by auto
  define c where "c = Lcm ((nat ∘ snd ∘ quotient_of ∘ poly.coeff p') ` {..Polynomial.degree p'})"
  have "¬snd (quotient_of x) ≤ 0" for x
    using quotient_of_denom_pos[of x, OF surjective_pairing] by auto
  hence "c ≠ 0" by (auto simp: c_def)
  define q where "q = Polynomial.smult (of_nat c) p'"

  have "poly.coeff q i ∈ ℤ" for i
  proof (cases "i > Polynomial.degree p'")
    case False
    define m n 
      where "m = fst (quotient_of (poly.coeff p' i))"
        and "n = nat (snd (quotient_of (poly.coeff p' i)))"
    have mn: "n > 0" "poly.coeff p' i = of_int m / of_nat n"
      using quotient_of_denom_pos[of "poly.coeff p' i", OF surjective_pairing]
            quotient_of_div[of "poly.coeff p' i", OF surjective_pairing]
      by (auto simp: m_def n_def)
    from False have "n dvd c" unfolding c_def
      by (intro dvd_Lcm) (auto simp: c_def n_def o_def not_less)
    hence "of_nat c * (of_int m / of_nat n) = (of_nat (c div n) * of_int m :: rat)"
      by (auto simp: of_nat_div)
    also have "… ∈ ℤ" by auto
    finally show ?thesis using mn by (auto simp: q_def)
  qed (auto simp: q_def coeff_eq_0)
  with intpolyE obtain q' where q': "q = of_int_poly q'" by auto
  moreover have "p = Polynomial.smult (inverse (of_nat c)) (map_poly of_rat (of_int_poly q'))"
    unfolding smult_conv_map_poly q'[symmetric] p' using ‹c ≠ 0›
    by (intro poly_eqI) (auto simp: coeff_map_poly q_def of_rat_mult)
  ultimately show ?thesis
    using q' p' ‹c ≠ 0› by (auto intro!: that[of c q'])
qed

lemma symmetric_mpoly_symmetric_sum:
  assumes "⋀π. π permutes A ⟹ g π permutes X"
  assumes "⋀x π. x ∈ X ⟹ π permutes A ⟹ mpoly_map_vars π (f x) = f (g π x)"
  shows "symmetric_mpoly A (∑x∈X. f x)"
  unfolding symmetric_mpoly_def
proof safe
  fix π assume π: "π permutes A"
  have "mpoly_map_vars π (sum f X) = (∑x∈X. mpoly_map_vars π (f x))"
    by simp
  also have "… = (∑x∈X. f (g π x))"
    by (intro sum.cong assms π refl)
  also have "… = (∑x∈g π`X. f x)"
    using assms(1)[OF π] by (subst sum.reindex) (auto simp: permutes_inj_on)
  also have "g π ` X = X"
    using assms(1)[OF π] by (simp add: permutes_image)
  finally show "mpoly_map_vars π (sum f X) = sum f X" .
qed

(* TODO: The version of this theorem in the AFP is to weak and should be replaced by this one. *)
lemma symmetric_mpoly_symmetric_prod:
  assumes "g permutes X"
  assumes "⋀x π. x ∈ X ⟹ π permutes A ⟹ mpoly_map_vars π (f x) = f (g x)"
  shows "symmetric_mpoly A (∏x∈X. f x)"
  unfolding symmetric_mpoly_def
proof safe
  fix π assume π: "π permutes A"
  have "mpoly_map_vars π (prod f X) = (∏x∈X. mpoly_map_vars π (f x))"
    by simp
  also have "… = (∏x∈X. f (g x))"
    by (intro prod.cong assms π refl)
  also have "… = (∏x∈g`X. f x)"
    using assms by (subst prod.reindex) (auto simp: permutes_inj_on)
  also have "g ` X = X"
    using assms by (simp add: permutes_image)
  finally show "mpoly_map_vars π (prod f X) = prod f X" .
qed


text ‹
  We now prove the transcendence of $i\pi$, from which the transcendence of $\pi$ will follow
  as a trivial corollary. The first proof of this was given by von Lindemann~\<^cite>‹"lindemann_pi82"›.
  The central ingredient is the fundamental theorem of symmetric functions.

  The proof can, by now, be considered folklore and one can easily find many similar variants of
  it, but we mostly follows the nice exposition given by Niven~\<^cite>‹"niven_pi39"›.

  An independent previous formalisation in Coq that uses the same basic techniques was given by
  Bernard et al.~\<^cite>‹"bernard_pi16"›. They later also formalised the much stronger
  Lindemann--Weierstra{\ss} theorem~\<^cite>‹"bernard_lw17"›.
›
lemma transcendental_i_pi: "¬algebraic (𝗂 * pi)"
proof
  ― ‹Suppose $i\pi$ were algebraic.›
  assume "algebraic (𝗂 * pi)"
  ― ‹We obtain some nonzero integer polynomial that has $i\pi$ as a root. We can assume
      w.\,l.\,o.\,g.\ that the constant coefficient of this polynomial is nonzero.›
  then obtain p
    where p: "poly (of_int_poly p) (𝗂 * pi) = 0" "p ≠ 0" "poly.coeff p 0 ≠ 0"
    by (elim algebraicE'_nonzero) auto
  define n where "n = Polynomial.degree p"

  ― ‹We define the sequence of the roots of this polynomial:›
  obtain root where "Polynomial.smult (Polynomial.lead_coeff (of_int_poly p))
                       (∏i<n. [:-root i :: complex, 1:]) = of_int_poly p"
    using complex_poly_decompose'[of "of_int_poly p"] unfolding n_def by auto
  note root = this [symmetric]

  ― ‹We note that $i\pi$ is, of course, among these roots.›
  from p and root obtain idx where idx: "idx < n" "root idx = 𝗂 * pi"
    by (auto simp: poly_prod)

  ― ‹We now define a new polynomial ‹P'›, whose roots are all numbers that arise as a
      sum of any subset of roots of ‹p›. We also count all those subsets that sum up to 0
      and call their number ‹A›.›
  define root' where "root' = (λX. (∑j∈X. root j))"
  define P where "P = (λi. ∏X | X ⊆ {..<n} ∧ card X = i. [:-root' X, 1:])"
  define P' where "P' = (∏i∈{0<..n}. P i)"
  define A where "A = card {X∈Pow {..<n}. root' X = 0}"
  have [simp]: "P' ≠ 0" by (auto simp: P'_def P_def)

  ― ‹We give the name ‹Roots'› to those subsets that do not sum to zero and note that there
      is at least one, namely $\{i\pi\}$.›
  define Roots' where "Roots' = {X. X ⊆ {..<n} ∧ root' X ≠ 0}"
  have [intro]: "finite Roots'" by (auto simp: Roots'_def)
  have "{idx} ∈ Roots'" using idx by (auto simp: Roots'_def root'_def)
  hence "Roots' ≠ {}" by auto
  hence card_Roots': "card Roots' > 0" by (auto simp: card_eq_0_iff)

  have P'_altdef: "P' = (∏X∈Pow {..<n} - {{}}. [:-root' X, 1:])"
  proof -
    have "P' = (∏(i, X)∈(SIGMA x:{0<..n}. {X. X ⊆ {..<n} ∧ card X = x}). [:- root' X, 1:])"
      unfolding P'_def P_def by (subst prod.Sigma) auto
    also have "… = (∏X∈Pow{..<n} - {{}}. [:- root' X, 1:])"
      using card_mono[of "{..<n}"]
      by (intro prod.reindex_bij_witness[of _ "λX. (card X, X)" "λ(_, X). X"])
         (auto simp: case_prod_unfold card_gt_0_iff intro: finite_subset[of _ "{..<n}"])
    finally show ?thesis .
  qed

  ― ‹Clearly, @{term A} is nonzero, since the empty set sums to 0.›
  have "A > 0"
  proof -
    have "{} ∈ {X∈Pow {..<n}. root' X = 0}"
      by (auto simp: root'_def)
    thus ?thesis by (auto simp: A_def card_gt_0_iff)
  qed

  ― ‹Since $e^{i\pi} + 1 = 0$, we know the following:›
  have "0 = (∏i<n. exp (root i) + 1)"
    using idx by force
  ― ‹We rearrange this product of sums into a sum of products and collect all summands that are
      1 into a separate sum, which we call @{term A}:›
  also have "… = (∑X∈Pow {..<n}. ∏i∈X. exp (root i))"
    by (subst prod_add) auto
  also have "… = (∑X∈Pow {..<n}. exp (root' X))"
    by (intro sum.cong refl, subst exp_sum [symmetric])
       (auto simp: root'_def intro: finite_subset[of _ "{..<n}"])
  also have "Pow {..<n} = {X∈Pow {..<n}. root' X ≠ 0} ∪ {X∈Pow {..<n}. root' X = 0}"
    by auto
  also have "(∑X∈…. exp (root' X)) = (∑X | X ⊆ {..<n} ∧ root' X ≠ 0. exp (root' X)) +
                                       (∑X | X ⊆ {..<n} ∧ root' X = 0. exp (root' X))"
    by (subst sum.union_disjoint) auto
  also have "(∑X | X ⊆ {..<n} ∧ root' X = 0. exp (root' X)) = of_nat A"
    by (simp add: A_def)
  ― ‹Finally, we obtain the fact that the sum of $\exp(u)$ with $u$ ranging over all the non-zero
      roots of @{term P'} is a negative integer.›
  finally have eq: "(∑X | X ⊆ {..<n} ∧ root' X ≠ 0. exp (root' X)) = -of_nat A"
    by (simp add: add_eq_0_iff2)

  ― ‹Next, we show that ‹P'› is a rational polynomial since it can be written as a symmetric
      polynomial expression (with rational coefficients) in the roots of ‹p›.›
  define ratpolys where "ratpolys = {p::complex poly. ∀i. poly.coeff p i ∈ ℚ}"
  have ratpolysI: "p ∈ ratpolys" if "⋀i. poly.coeff p i ∈ ℚ" for p
    using that by (auto simp: ratpolys_def)

  have "P' ∈ ratpolys"
  proof -
    define Pmv :: "nat ⇒ complex poly mpoly"
      where "Pmv = (λi. ∏X | X ⊆ {..<n} ∧ card X = i. Const ([:0,1:]) -
                          (∑i∈X. monom (Poly_Mapping.single i 1) 1))"
    define P'mv where "P'mv = (∏i∈{0<..n}. Pmv i)"
    have "insertion (λi. [:root i:]) P'mv ∈ ratpolys"
    proof (rule symmetric_poly_of_roots_in_subring[where l = "λx. [:x:]"])
      show "ring_closed ratpolys"
        by standard (auto simp: ratpolys_def coeff_mult)
      then interpret ring_closed ratpolys .
      show "∀m. coeff P'mv m ∈ ratpolys"
        by (auto simp: P'mv_def Pmv_def coeff_monom when_def mpoly_coeff_Const coeff_pCons' ratpolysI
                 intro!: coeff_prod_closed minus_closed sum_closed uminus_closed)
      show "∀i. [:poly.coeff (of_int_poly p) i:] ∈ ratpolys"
        by (intro ratpolysI allI) (auto simp: coeff_pCons')
      show "[:inverse (of_int (Polynomial.lead_coeff p)):] *
              [:of_int (Polynomial.lead_coeff p) :: complex:] = 1"
        using ‹p ≠ 0› by (auto intro!: poly_eqI simp: field_simps)
    next
      have "symmetric_mpoly {..<n} (Pmv k)" for k
        unfolding symmetric_mpoly_def
      proof safe
        fix π :: "nat ⇒ nat" assume π: "π permutes {..<n}"
        hence "mpoly_map_vars π (Pmv k) =
                 (∏X | X ⊆ {..<n} ∧ card X = k. Const [:0, 1:] -
                    (∑x∈X. MPoly_Type.monom (Poly_Mapping.single (π x) (Suc 0)) 1))"
          by (simp add: Pmv_def permutes_bij)
        also have "… = (∏X | X ⊆ {..<n} ∧ card X = k. Const [:0, 1:] -
                          (∑x∈π`X. MPoly_Type.monom (Poly_Mapping.single x (Suc 0)) 1))"
          using π by (subst sum.reindex) (auto simp: permutes_inj_on)
        also have "… = (∏X ∈ (λX. π`X)`{X. X ⊆ {..<n} ∧ card X = k}. Const [:0, 1:] -
                          (∑x∈X. MPoly_Type.monom (Poly_Mapping.single x (Suc 0)) 1))"
          by (subst prod.reindex) (auto intro!: inj_on_image permutes_inj_on[OF π])
        also have "(λX. π`X)`{X. X ⊆ {..<n} ∧ card X = k} = {X. X ⊆ π ` {..<n} ∧ card X = k}"
          using π by (subst image_image_fixed_card_subset) (auto simp: permutes_inj_on)
        also have "π ` {..<n} = {..<n}"
          by (intro permutes_image π)
        finally show "mpoly_map_vars π (Pmv k) = Pmv k" by (simp add: Pmv_def)
      qed
      thus "symmetric_mpoly {..<n} P'mv"
        unfolding P'mv_def by (intro symmetric_mpoly_prod) auto
    next
      show vars_P'mv: "vars P'mv ⊆ {..<n}"
        unfolding P'mv_def Pmv_def
        by (intro order.trans[OF vars_prod] UN_least order.trans[OF vars_diff]
                   Un_least order.trans[OF vars_sum] order.trans[OF vars_monom_subset]) auto
    qed (insert root, auto intro!: ratpolysI simp: coeff_pCons')
    also have "insertion (λi. [:root i:]) (Pmv k) = P k" for k
      by (simp add: Pmv_def insertion_prod insertion_diff insertion_sum root'_def P_def
                    sum_to_poly del: insertion_monom)
      (* TODO: insertion_monom should not be a simp rule *)
    hence "insertion (λi. [:root i:]) P'mv = P'"
      by (simp add: P'mv_def insertion_prod P'_def)
    finally show "P' ∈ ratpolys" .
  qed

  ― ‹We clear the denominators and remove all powers of $X$ from @{term P'} to obtain a new
      integer polynomial ‹Q›.›
  define Q' where "Q' = (∏X∈Roots'. [:- root' X, 1:])"
  have "P' = (∏X∈Pow {..<n}-{{}}. [:-root' X, 1:])"
    by (simp add: P'_altdef)
  also have "Pow {..<n}-{{}} = Roots' ∪
               {X. X ∈ Pow {..<n} - {{}} ∧ root' X = 0}" by (auto simp: root'_def Roots'_def)
  also have "(∏X∈…. [:-root' X, 1:]) =
               Q' * [:0, 1:] ^ card {X. X ⊆ {..<n} ∧ X ≠ {} ∧ root' X = 0}"
    by (subst prod.union_disjoint) (auto simp: Q'_def Roots'_def)
  also have "{X. X ⊆ {..<n} ∧ X ≠ {} ∧ root' X = 0} = {X. X ⊆ {..<n} ∧ root' X = 0} - {{}}"
    by auto
  also have "card … = A - 1"  unfolding A_def
    by (subst card_Diff_singleton) (auto simp: root'_def)
  finally have Q': "P' = Polynomial.monom 1 (A - 1) * Q'"
    by (simp add: Polynomial.monom_altdef)
  have degree_Q': "Polynomial.degree P' = Polynomial.degree Q' + (A - 1)"
    by (subst Q')
       (auto simp: Q'_def Roots'_def degree_mult_eq Polynomial.degree_monom_eq degree_prod_sum_eq)

  have "∀i. poly.coeff Q' i ∈ ℚ"
  proof
    fix i :: nat
    have "poly.coeff Q' i = Polynomial.coeff P' (i + (A - 1))"
      by (simp add: Q' Polynomial.coeff_monom_mult)
    also have "… ∈ ℚ" using ‹P' ∈ ratpolys› by (auto simp: ratpolys_def)
    finally show "poly.coeff Q' i ∈ ℚ" .
  qed
  from ratpoly_to_intpoly[OF this] obtain c Q
    where [simp]: "c ≠ 0" and Q: "Q' = Polynomial.smult (inverse (of_nat c)) (of_int_poly Q)"
    by metis
  have [simp]: "Q ≠ 0" using Q Q' by auto
  have Q': "of_int_poly Q = Polynomial.smult (of_nat c) Q'"
    using Q by simp
  have degree_Q: "Polynomial.degree Q = Polynomial.degree Q'"
    by (subst Q) auto
  have "Polynomial.lead_coeff (of_int_poly Q :: complex poly) = c"
    by (subst Q') (simp_all add: degree_Q Q'_def Polynomial.lead_coeff_prod)
  hence lead_coeff_Q: "Polynomial.lead_coeff Q = int c"
    using of_int_eq_iff[of "Polynomial.lead_coeff Q" "of_nat c"] by (auto simp del: of_int_eq_iff)  
  have Q_decompose: "of_int_poly Q =
               Polynomial.smult (of_nat c) (∏X∈Roots'. [:- root' X, 1:])"
    by (subst Q') (auto simp: Q'_def lead_coeff_Q)
  have "poly (of_int_poly Q) (𝗂 * pi) = 0"
    using ‹{idx} ∈ Roots'› ‹finite Roots'› idx 
    by (force simp: root'_def Q_decompose poly_prod)
  have degree_Q: "Polynomial.degree (of_int_poly Q :: complex poly) = card Roots'"
    by (subst Q') (auto simp: Q'_def degree_prod_sum_eq)
  have "poly (of_int_poly Q) (0 :: complex) ≠ 0"
    by (subst Q') (auto simp: Q'_def Roots'_def poly_prod)
  hence [simp]: "poly Q 0 ≠ 0" by simp
  have [simp]: "poly (of_int_poly Q) (root' Y) = 0" if "Y ∈ Roots'" for Y
    using that ‹finite Roots'› by (auto simp: Q' Q'_def poly_prod)

  ― ‹We find some closed ball that contains all the roots of @{term Q}.›
  define r where "r = Polynomial.degree Q"
  have "r > 0" using degree_Q card_Roots' by (auto simp: r_def)
  define Radius where "Radius = Max ((λY. norm (root' Y)) ` Roots')"
  have Radius: "norm (root' Y) ≤ Radius" if "Y ∈ Roots'" for Y
    using ‹finite Roots'› that by (auto simp: Radius_def)
  from Radius[of "{idx}"] have "Radius ≥ pi"
    using idx by (auto simp: Roots'_def norm_mult root'_def)
  hence Radius_nonneg: "Radius ≥ 0" and "Radius > 0" using pi_gt3 by linarith+

  ― ‹Since this ball is compact, @{term Q} is bounded on it. We obtain such a bound.›
  have "compact (poly (of_int_poly Q :: complex poly) ` cball 0 Radius)"
    by (intro compact_continuous_image continuous_intros) auto
  then obtain Q_ub
    where Q_ub: "Q_ub > 0" 
                "⋀u :: complex. u ∈ cball 0 Radius ⟹ norm (poly (of_int_poly Q) u) ≤ Q_ub"
    by (auto dest!: compact_imp_bounded simp: bounded_pos cball_def)

  ― ‹Using this, define another upper bound that we will need later.›
  define fp_ub
    where "fp_ub = (λp. ¦c¦ ^ (r * p - 1) / fact (p - 1) * (Radius ^ (p - 1) * Q_ub ^ p))"
  have fp_ub_nonneg: "fp_ub p ≥ 0" for p
    unfolding fp_ub_def using ‹Radius ≥ 0› Q_ub
    by (intro mult_nonneg_nonneg divide_nonneg_pos zero_le_power) auto
  define C where "C = card Roots' * Radius * exp Radius"

  ― ‹We will now show that any sufficiently large prime number leads to
      ‹C * fp_ub p ≥ 1›, from which we will then derive a contradiction.›
  define primes_at_top where "primes_at_top = inf_class.inf sequentially (principal {p. prime p})"
  have "eventually (λp. ∀x∈{nat ¦poly Q 0¦, c, A}. p > x) sequentially"
    by (intro eventually_ball_finite ballI eventually_gt_at_top) auto
  hence "eventually (λp. ∀x∈{nat ¦poly Q 0¦, c, A}. p > x) primes_at_top"
    unfolding primes_at_top_def eventually_inf_principal by eventually_elim auto
  moreover have "eventually (λp. prime p) primes_at_top"
    by (auto simp: primes_at_top_def eventually_inf_principal)
  ultimately have "eventually (λp. C * fp_ub p ≥ 1) primes_at_top"
  proof eventually_elim
    case (elim p)
    hence p: "prime p" "p > nat ¦poly Q 0¦" "p > c" "p > A" by auto
    hence "p > 1" by (auto dest: prime_gt_1_nat)

    ― ‹We define the polynomial $f(X) = \frac{c^s}{(p-1)!} X^{p-1} Q(X)^p$, where $c$ is
        the leading coefficient of $Q$. We also define $F(X)$ to be the sum of all its
        derivatives.›
    define s where "s = r * p - 1"
    define fp :: "complex poly"
      where "fp = Polynomial.smult (of_nat c ^ s / fact (p - 1))
                    (Polynomial.monom 1 (p - 1) * of_int_poly Q ^ p)"
    define Fp where "Fp = (∑i≤s+p. (pderiv ^^ i) fp)"
    define f F where "f = poly fp" and "F = poly Fp"
    have degree_fp: "Polynomial.degree fp = s + p" using degree_Q card_Roots' ‹p > 1›
      by (simp add: fp_def s_def degree_mult_eq degree_monom_eq
                    degree_power_eq r_def algebra_simps)

    ― ‹Using the same argument as in the case of the transcendence of $e$, we now 
        consider the function
        \[I(u) := e^u F(0) - F(u) = u \int\limits_0^1 e^{(1-t)x} f(tx)\,\textrm{d}t\]
        whose absolute value can be bounded with a standard ``maximum times length'' estimate
        using our upper bound on $f$. All of this can be reused from the proof for $e$, so
        there is not much to do here. In particular, we will look at $\sum I(x_i)$ with the
        $x_i$ ranging over the roots of $Q$ and bound this sum in two different ways.›
    interpret lindemann_weierstrass_aux fp .
    have I_altdef: "I = (λu. exp u * F 0 - F u)"
      by (intro ext) (simp add: I_def degree_fp F_def Fp_def poly_sum)

    ― ‹We show that @{term fp_ub} is indeed an upper bound for $f$.›
    have fp_ub:  "norm (poly fp u) ≤ fp_ub p" if "u ∈ cball 0 Radius" for u
    proof -
      have "norm (poly fp u) = ¦c¦ ^ (r * p - 1) / fact (p - 1) * (norm u ^ (p - 1) *
                                 norm (poly (of_int_poly Q) u) ^ p)"
        by (simp add: fp_def f_def s_def norm_mult poly_monom norm_divide norm_power)
      also have "… ≤ fp_ub p"
        unfolding fp_ub_def using that Q_ub ‹Radius ≥ 0›
        by (intro mult_left_mono[OF mult_mono] power_mono zero_le_power) auto
      finally show ?thesis .
    qed    

    ― ‹We now show that the following sum is an integer multiple of $p$. This argument again
        uses the fundamental theorem of symmetric functions, exploiting that the inner sums are
        symmetric over the roots of $Q$.›
    have "(∑i=p..s+p. ∑Y∈Roots'. poly ((pderiv ^^ i) fp) (root' Y)) / p ∈ ℤ"
    proof (subst sum_divide_distrib, intro Ints_sum[of "{a..b}" for a b])
      fix i assume i: "i ∈ {p..s+p}"
      then obtain roots' where roots': "distinct roots'" "set roots' = Roots'"
        using finite_distinct_list ‹finite Roots'› by metis
      define l where "l = length roots'"
      define fp' where "fp' = (pderiv ^^ i) fp"
      define d where "d = Polynomial.degree fp'"
      ― ‹We define a multivariate polynomial for the inner sum $\sum f(x_i)/p$ in order
          to show that it is indeed a symmetric function over the $x_i$.›
      define R where "R = (smult (1 / of_nat p) (∑k≤d. ∑i<l. smult (poly.coeff fp' k)
                             (monom (Poly_Mapping.single i k) (1 / of_int (c ^ k)))) :: complex mpoly)"

      ― ‹The $j$-th coefficient of the $i$-th derivative of $f$ are integer multiples
          of $c^j p$ since $i \geq p$.›
      have integer: "poly.coeff fp' j / (of_nat c ^ j * of_nat p) ∈ ℤ" if "j ≤ d" for j
      proof -
        define fp'' where "fp'' = Polynomial.monom 1 (p - 1) * Q ^ p"
        define x
          where "x = c ^ s * poly.coeff ((pderiv ^^ i) (Polynomial.monom 1 (p - 1) * Q ^ p)) j"
        have "[:fact p:] dvd ([:fact i:] :: int poly)" using i
          by (auto intro: fact_dvd)
        also have "[:fact i:] dvd ((pderiv ^^ i) (Polynomial.monom 1 (p - 1) * Q ^ p))"
          by (rule fact_dvd_higher_pderiv)
        finally have "c ^ j * fact p dvd x" unfolding x_def of_nat_mult using that i
          by (intro mult_dvd_mono)
             (auto intro!: le_imp_power_dvd simp: s_def d_def fp'_def degree_higher_pderiv degree_fp)
        hence "of_int x / (of_int (c ^ j * fact p) :: complex) ∈ ℤ"
          by (intro of_int_divide_in_Ints) auto
        also have "of_int x / (of_int (c ^ j * fact p) :: complex) =
                     poly.coeff fp' j / (of_nat c ^ j * of_nat p)" using ‹p > 1›
          unfolding x_def fp'_def fp_def
          by (simp add: fact_reduce[of p] field_simps hom_distribs higher_pderiv_smult
                   flip: of_int_hom.coeff_map_poly_hom)
        finally show ?thesis .
      qed

      ― ‹Evaluating $R$ yields is an integer since it is symmetric.›
      have "insertion (λi. c * root' (roots' ! i)) R ∈ ℤ"
      proof (intro symmetric_poly_of_roots_in_subring_monic allI)
        define Q' where "Q' = of_int_poly Q ∘⇩p [:0, 1 / of_nat c :: complex:]"
        show "symmetric_mpoly {..<l} R" unfolding R_def
          by (intro symmetric_mpoly_smult symmetric_mpoly_sum[of "{..d}"] symmetric_mpoly_symmetric_sum)
             (simp_all add: mpoly_map_vars_monom permutes_bij permutep_single bij_imp_bij_inv permutes_inv_inv)
        show "MPoly_Type.coeff R m ∈ ℤ" for m unfolding R_def coeff_sum coeff_smult sum_distrib_left
          using integer by (auto simp: R_def coeff_monom when_def intro!: Ints_sum)
        show "vars R ⊆ {..<l}" unfolding R_def
          by (intro order.trans[OF vars_smult] order.trans[OF vars_sum] UN_least
                    order.trans[OF vars_monom_subset]) auto
        show "ring_closed ℤ" by standard auto
  
        have "(∏i<l. [:- (of_nat c * root' (roots' ! i)), 1:]) =
                (∏Y←roots'. [:- (of_nat c * root' Y), 1:])"
          by (subst prod.list_conv_set_nth) (auto simp: atLeast0LessThan l_def)
        also have "… = (∏Y∈Roots'. [:- (of_nat c * root' Y), 1:])"
          using roots' by (subst prod.distinct_set_conv_list [symmetric]) auto
        also have "… = (∏Y∈Roots'. Polynomial.smult (of_nat c) ([:-root' Y, 1:])) ∘⇩p [:0, 1 / c:]"
          by (simp add: pcompose_prod pcompose_pCons)
        also have "(∏Y∈Roots'. Polynomial.smult (of_nat c) ([:-root' Y, 1:])) =
                     Polynomial.smult (of_nat c ^ card Roots') (∏Y∈Roots'. [:-root' Y, 1:])"
          by (subst prod_smult) auto
        also have "… = Polynomial.smult (of_nat c ^ (card Roots' - 1))
                          (Polynomial.smult c (∏Y∈Roots'. [:-root' Y, 1:]))"
          using ‹finite Roots'› and ‹Roots' ≠ {}›
          by (subst power_diff) (auto simp: Suc_le_eq card_gt_0_iff)
        also have "Polynomial.smult c (∏Y∈Roots'. [:-root' Y, 1:]) = of_int_poly Q"
          using Q_decompose by simp
        finally show "Polynomial.smult (of_nat c ^ (card Roots' - 1)) Q' =
                        (∏i<l. [:- (of_nat c * root' (roots' ! i)), 1:])"
          by (simp add: pcompose_smult Q'_def)
        fix i :: nat
        show "poly.coeff (Polynomial.smult (of_nat c ^ (card Roots' - 1)) Q') i ∈ ℤ"
        proof (cases i "Polynomial.degree Q" rule: linorder_cases)
          case greater
          thus ?thesis by (auto simp: Q'_def coeff_pcompose_linear coeff_eq_0)
        next
          case equal
          thus ?thesis using ‹Roots' ≠ {}› degree_Q card_Roots' lead_coeff_Q
            by (auto simp: Q'_def coeff_pcompose_linear lead_coeff_Q power_divide power_diff)
        next
          case less
          have "poly.coeff (Polynomial.smult (of_nat c ^ (card Roots' - 1)) Q') i =
                  of_int (poly.coeff Q i) * (of_int (c ^ (card Roots' - 1)) / of_int (c ^ i))"
            by (auto simp: Q'_def coeff_pcompose_linear power_divide)
          also have "… ∈ ℤ" using less degree_Q
            by (intro Ints_mult of_int_divide_in_Ints) (auto intro!: le_imp_power_dvd)
          finally show ?thesis .
        qed
      qed auto
      ― ‹Moreover, by definition, evaluating @{term R} gives us $\sum f(x_i)/p$.›
      also have "insertion (λi. c * root' (roots' ! i)) R =
                   (∑Y←roots'. poly fp' (root' Y)) / of_nat p"
        by (simp add: insertion_sum R_def poly_altdef d_def sum_list_sum_nth atLeast0LessThan
                    l_def power_mult_distrib algebra_simps 
                    sum.swap[of _ "{..Polynomial.degree fp'}"] del: insertion_monom)
      also have "… = (∑Y∈Roots'. poly ((pderiv ^^ i) fp) (root' Y)) / of_nat p"
        using roots' by (subst sum_list_distinct_conv_sum_set) (auto simp: fp'_def poly_pcompose)
      finally show "… ∈ ℤ" .
    qed
    then obtain K where K: "(∑i=p..s+p. ∑Y∈Roots'. 
                               poly ((pderiv ^^ i) fp) (root' Y)) = of_int K * p"
      using ‹p > 1› by (auto elim!: Ints_cases simp: field_simps)

    ― ‹Next, we show that $F(0)$ is an integer and coprime to $p$.›
    obtain F0 :: int where F0: "F 0 = of_int F0" "coprime (int p) F0"
    proof -
      have "(∑i=p..s + p. poly ((pderiv ^^ i) fp) 0) / of_nat p ∈ ℤ"
        unfolding sum_divide_distrib
      proof (intro Ints_sum)
        fix i assume i: "i ∈ {p..s+p}"
        hence "fact p dvd poly ((pderiv ^^ i) ([:0, 1:] ^ (p - 1) * Q ^ p)) 0"
          by (intro fact_dvd_poly_higher_pderiv_aux') auto
        then obtain k where k: "poly ((pderiv ^^ i) ([:0, 1:] ^ (p - 1) * Q ^ p)) 0 = k * fact p"
          by auto
  
        have "(pderiv ^^ i) fp = Polynomial.smult (of_nat c ^ s / fact (p - 1))
                (of_int_poly ((pderiv ^^ i) ([:0, 1:] ^ (p - 1) * Q ^ p)))"
          by (simp add: fp_def higher_pderiv_smult Polynomial.monom_altdef hom_distribs)
        also have "poly … 0 / of_nat p = of_int (c ^ s * k)"
          using k ‹p > 1› by (simp add: fact_reduce[of p])
        also have "… ∈ ℤ" by simp
        finally show "poly ((pderiv ^^ i) fp) 0 / of_nat p ∈ ℤ" .
      qed
      then obtain S where S: "(∑i=p..s + p. poly ((pderiv ^^ i) fp) 0) = of_int S * p"
        using ‹p > 1› by (auto elim!: Ints_cases simp: field_simps)
  
      have "F 0 = (∑i≤s + p. poly ((pderiv ^^ i) fp) 0)"
        by (auto simp: F_def Fp_def poly_sum)
      also have "… = (∑i∈insert (p - 1) {p..s + p}. poly ((pderiv ^^ i) fp) 0)"
      proof (intro sum.mono_neutral_right ballI)
        fix i assume i: "i ∈ {..s + p} - insert (p - 1) {p..s + p}"
        hence "i < p - 1" by auto
        have "Polynomial.monom 1 (p - 1) dvd fp"
          by (auto simp: fp_def intro: dvd_smult)
        with i show "poly ((pderiv ^^ i) fp) 0 = 0"
          by (intro poly_higher_pderiv_aux1'[of _ "p - 1"]) (auto simp: Polynomial.monom_altdef)
      qed auto
      also have "… = poly ((pderiv ^^ (p - 1)) fp) 0 + of_int S * of_nat p"
        using ‹p > 1› S by (subst sum.insert) auto
      also have "poly ((pderiv ^^ (p - 1)) fp) 0 = of_int (c ^ s * poly Q 0 ^ p)"
        using poly_higher_pderiv_aux2[of "p - 1" 0 "of_int_poly Q ^ p :: complex poly"]
        by (simp add: fp_def higher_pderiv_smult Polynomial.monom_altdef)
      finally have "F 0 = of_int (S * int p + c ^ s * poly Q 0 ^ p)"
        by simp
      moreover have "coprime p c" "coprime (int p) (poly Q 0)"
        using p by (auto intro!: prime_imp_coprime dest: dvd_imp_le_int[rotated])
      hence "coprime (int p) (c ^ s * poly Q 0 ^ p)"
        by auto
      hence "coprime (int p) (S * int p + c ^ s * poly Q 0 ^ p)"
        unfolding coprime_iff_gcd_eq_1 gcd_add_mult by auto
      ultimately show ?thesis using that[of "S * int p + c ^ s * poly Q 0 ^ p"] by blast
    qed

    ― ‹Putting everything together, we have shown that $\sum I(x_i)$ is an integer coprime
        to $p$, and therefore a nonzero integer, and therefore has an absolute value of at least 1.›
    have "(∑Y∈Roots'. I (root' Y)) = F 0 * (∑Y∈Roots'. exp (root' Y)) - (∑Y∈Roots'. F (root' Y))"
      by (simp add: I_altdef sum_subtractf sum_distrib_left sum_distrib_right algebra_simps)
    also have "… = -(of_int (F0 * int A) +
                      (∑i≤s+p. ∑Y∈Roots'. poly ((pderiv ^^ i) fp) (root' Y)))"
      using F0 by (simp add: Roots'_def eq F_def Fp_def poly_sum sum.swap[of _ "{..s+p}"])
    also have "(∑i≤s+p. ∑Y∈Roots'. poly ((pderiv ^^ i) fp) (root' Y)) =
                 (∑i=p..s+p. ∑Y∈Roots'. poly ((pderiv ^^ i) fp) (root' Y))"
    proof (intro sum.mono_neutral_right ballI sum.neutral)
      fix i Y assume i: "i ∈ {..s+p} - {p..s+p}" and Y: "Y ∈ Roots'"
      have "[:-root' Y, 1:] ^ p dvd of_int_poly Q ^ p"
        by (intro dvd_power_same) (auto simp: dvd_iff_poly_eq_0 Y)
      hence "[:-root' Y, 1:] ^ p dvd fp"
        by (auto simp: fp_def intro!: dvd_smult)
      thus "poly ((pderiv ^^ i) fp) (root' Y) = 0"
        using i by (intro poly_higher_pderiv_aux1') auto
    qed auto
    also have "… = of_int (K * int p)" using K by simp
    finally have "(∑Y∈Roots'. I (root' Y)) = -of_int (K * int p + F0 * int A)"
      by simp
    moreover have "coprime p A"
      using p ‹A > 0› by (intro prime_imp_coprime) (auto dest!: dvd_imp_le)
    hence "coprime (int p) (F0 * int A)"
      using F0 by auto
    hence "coprime (int p) (K * int p + F0 * int A)"
      using F0 unfolding coprime_iff_gcd_eq_1 gcd_add_mult by auto
    hence "K * int p + F0 * int A ≠ 0"
      using p by (intro notI) auto
    hence "norm (-of_int (K * int p + F0 * int A) :: complex) ≥ 1"
      unfolding norm_minus_cancel norm_of_int by linarith
    ultimately have "1 ≤ norm (∑Y∈Roots'. I (root' Y))" by metis

    ― ‹The M--L bound on the integral gives us an upper bound:›
    also have "norm (∑Y∈Roots'. I (root' Y)) ≤
                 (∑Y∈Roots'. norm (root' Y) * exp (norm (root' Y)) * fp_ub p)"
    proof (intro sum_norm_le lindemann_weierstrass_integral_bound fp_ub fp_ub_nonneg)
      fix Y u assume *: "Y ∈ Roots'" "u ∈ closed_segment 0 (root' Y)"
      hence "closed_segment 0 (root' Y) ⊆ cball 0 Radius"
        using ‹Radius ≥ 0› Radius[of Y] by (intro closed_segment_subset) auto
      with * show "u ∈ cball 0 Radius" by auto
    qed
    also have "… ≤ (∑Y∈Roots'. Radius * exp (Radius) * fp_ub p)"
      using Radius by (intro sum_mono mult_right_mono mult_mono fp_ub_nonneg ‹Radius ≥ 0›) auto
    also have "… = C * fp_ub p" by (simp add: C_def)
    finally show "1 ≤ C * fp_ub p" .
  qed

  ― ‹It now only remains to show that this inequality is inconsistent for large $p$.
      This is obvious, since the upper bound is an exponential divided by a factorial and 
      therefore clearly tends to zero.›
  have "(λp. C * fp_ub p) ∈ Θ(λp. (C / (Radius * ¦c¦)) * (p / 2 ^ p) *
                                     ((2 * ¦c¦ ^ r * Radius * Q_ub) ^ p / fact p))"
    (is "_ ∈ Θ(?f)") using degree_Q card_Roots' ‹Radius > 0›
    by (intro bigthetaI_cong eventually_mono[OF eventually_gt_at_top[of 0]])
       (auto simp: fact_reduce power_mult [symmetric] r_def
                   fp_ub_def power_diff power_mult_distrib)
  also have "?f ∈ o(λp. 1 * 1 * 1)"
  proof (intro landau_o.big_small_mult landau_o.big_mult)
    have "(λx. (real_of_int (2 * ¦c¦ ^ r) * Radius * Q_ub) ^ x / fact x) ⇢ 0"
      by (intro power_over_fact_tendsto_0)
    thus "(λx. (real_of_int (2 * ¦c¦ ^ r) * Radius * Q_ub) ^ x / fact x) ∈ o(λx. 1)"
      by (intro smalloI_tendsto) auto
  qed real_asymp+
  finally have "(λp. C * fp_ub p) ∈ o(λ_. 1)" by simp
  from smalloD_tendsto[OF this] have "(λp. C * fp_ub p) ⇢ 0" by simp
  hence "eventually (λp. C * fp_ub p < 1) at_top"
    by (intro order_tendstoD) auto
  hence "eventually (λp. C * fp_ub p < 1) primes_at_top"
    unfolding primes_at_top_def eventually_inf_principal by eventually_elim auto
  moreover note ‹eventually (λp. C * fp_ub p ≥ 1) primes_at_top›
  ― ‹We therefore have a contradiction for any sufficiently large prime.›
  ultimately have "eventually (λp. False) primes_at_top"
    by eventually_elim auto

  ― ‹Since sufficiently large primes always exist, this concludes the theorem.›
  moreover have "frequently (λp. prime p) sequentially"
    using primes_infinite by (simp add: cofinite_eq_sequentially[symmetric] Inf_many_def)
  ultimately show False
    by (auto simp: frequently_def eventually_inf_principal primes_at_top_def)
qed

lemma pcompose_conjugates_integer:
  assumes "⋀i. poly.coeff p i ∈ ℤ"
  shows   "poly.coeff (pcompose p [:0, 𝗂:] * pcompose p [:0, -𝗂:]) i ∈ ℤ"
proof -
  let ?c = "λi. poly.coeff p i :: complex"
  have "poly.coeff (pcompose p [:0, 𝗂:] * pcompose p [:0, -𝗂:]) i =
          𝗂 ^ i * (∑k≤i. (-1) ^ (i - k) * ?c k * ?c (i - k))"
    unfolding coeff_mult sum_distrib_left
    by (intro sum.cong) (auto simp: coeff_mult coeff_pcompose_linear power_minus' 
                                    power_diff field_simps intro!: Ints_sum)
  also have "(∑k≤i. (-1) ^ (i - k) * ?c k * ?c (i - k)) =
          (∑k≤i. (-1) ^ k * ?c k * ?c (i - k))" (is "?S1 = ?S2")
    by (intro sum.reindex_bij_witness[of _ "λk. i - k" "λk. i - k"]) (auto simp: mult_ac)
  hence "?S1 = (?S1 + ?S2) / 2" by simp
  also have "… = (∑k≤i. ((-1) ^ k + (-1) ^ (i - k)) / 2 * ?c k * ?c (i - k))"
    by (simp add: ring_distribs sum.distrib sum_divide_distrib [symmetric])
  also have "… = (∑k≤i. (1 + (-1) ^ i) / 2 * (-1) ^ k * ?c k * ?c (i - k))"
    by (intro sum.cong) (auto simp: power_add power_diff field_simps)
  also have "𝗂 ^ i * … ∈ ℤ"
  proof (cases "even i")
    case True
    thus ?thesis
      by (intro Ints_mult Ints_sum assms) (auto elim!: evenE simp: power_mult)
  next
    case False
    hence "1 + (-1) ^ i = (0 :: complex)" by (auto elim!: oddE simp: power_mult)
    thus ?thesis by simp
  qed
  finally show ?thesis .
qed

lemma algebraic_times_i:
  assumes "algebraic x"
  shows   "algebraic (𝗂 * x)" "algebraic (-𝗂 * x)"
proof -
  from assms obtain p where p: "poly p x = 0" "∀i. poly.coeff p i ∈ ℤ" "p ≠ 0"
    by (auto elim!: algebraicE)
  define p' where "p' = pcompose p [:0, 𝗂:] * pcompose p [:0, -𝗂:]"
  have p': "poly p' (𝗂 * x) = 0" "poly p' (-𝗂 * x) = 0" "p' ≠ 0"
    by (auto simp: p'_def poly_pcompose algebra_simps p pcompose_eq_0_iff dest: pcompose_eq_0)
  moreover have "∀i. poly.coeff p' i ∈ ℤ"
    using p unfolding p'_def by (intro allI pcompose_conjugates_integer) auto
  ultimately show "algebraic (𝗂 * x)" "algebraic (-𝗂 * x)" by (intro algebraicI[of p']; simp)+
qed

lemma algebraic_times_i_iff: "algebraic (𝗂 * x) ⟷ algebraic x"
  using algebraic_times_i[of x] algebraic_times_i[of "𝗂 * x"] by auto

theorem transcendental_pi: "¬algebraic pi"
  using transcendental_i_pi by (simp add: algebraic_times_i_iff)

end