Theory Liouville_Numbers

(*
File:    Liouville_Numbers.thy
Author:  Manuel Eberl <manuel@pruvisto.org>

The definition of Liouville numbers and their standard construction, plus the proof
that any Liouville number is transcendental.
*)

theory Liouville_Numbers
imports
Complex_Main
"HOL-Computational_Algebra.Polynomial"
Liouville_Numbers_Misc
begin

(*
TODO: Move definition of algebraic numbers out of Algebraic_Numbers to reduce unnecessary
dependencies.
*)

text ‹
A Liouville number is a real number that can be approximated well -- but not perfectly --
by a sequence of rational numbers. Well, in this context, means that the error of the
$n$-th rational in the sequence is bounded by the $n$-th power of its denominator.

Our approach will be the following:
\begin{itemize}
\item Liouville numbers cannot be rational.
\item Any irrational algebraic number cannot be approximated in the Liouville sense
\item Therefore, all Liouville numbers are transcendental.
\item The standard construction fulfils all the properties of Liouville numbers.
\end{itemize}
›

subsection ‹Definition of Liouville numbers›

text ‹
The following definitions and proofs are largely adapted from those in the Wikipedia
article on Liouville numbers.~\<^cite>‹"wikipedia"›
›

text ‹
A Liouville number is a real number that can be approximated well -- but not perfectly --
by a sequence of rational numbers. The error of the $n$-th term $\frac{p_n}{q_n}$ is at most
$q_n^{-n}$, where $p_n\in\isasymint$ and $q_n \in\isasymint_{\geq 2}$.

We will say that such a number can be approximated in the Liouville sense.
›
locale liouville =
fixes x :: real and p q :: "nat ⇒ int"
assumes approx_int_pos: "abs (x - p n / q n) > 0"
and denom_gt_1:     "q n > 1"
and approx_int:     "abs (x - p n / q n) < 1 / of_int (q n) ^ n"

text ‹
First, we show that any Liouville number is irrational.
›
lemma (in liouville) irrational: "x ∉ ℚ"
proof
assume "x ∈ ℚ"
then obtain c d :: int where d: "x = of_int c / of_int d" "d > 0"
by (elim Rats_cases') simp
define n where "n = Suc (nat ⌈log 2 d⌉)"
note q_gt_1 = denom_gt_1[of n]

have neq: "c * q n ≠ d * p n"
proof
assume "c * q n = d * p n"
hence "of_int (c * q n) = of_int (d * p n)" by (simp only: )
with approx_int_pos[of n] d q_gt_1 show False by (auto simp: field_simps)
qed
hence abs_pos: "0 < abs (c * q n - d * p n)" by simp

from q_gt_1 d have "of_int d = 2 powr log 2 d" by (subst powr_log_cancel) simp_all
also from d have "log 2 (of_int d) ≥ log 2 1" by (subst log_le_cancel_iff) simp_all
hence "2 powr log 2 d ≤ 2 powr of_nat (nat ⌈log 2 d⌉)"
by (intro powr_mono) simp_all
also have "… = of_int (2 ^ nat ⌈log 2 d⌉)"
by (subst powr_realpow) simp_all
finally have "d ≤ 2 ^ nat ⌈log 2 (of_int d)⌉"
by (subst (asm) of_int_le_iff)
also have "… * q n ≤ q n ^ Suc (nat ⌈log 2 (of_int d)⌉)"
by (subst power_Suc, subst mult.commute, intro mult_left_mono power_mono,
insert q_gt_1) simp_all
also from q_gt_1 have "… = q n ^ n" by (simp add: n_def)
finally have "1 / of_int (q n ^ n) ≤ 1 / real_of_int (d * q n)" using q_gt_1 d
by (intro divide_left_mono Rings.mult_pos_pos of_int_pos, subst of_int_le_iff) simp_all
also have "… ≤ of_int (abs (c * q n - d * p n)) / real_of_int (d * q n)" using q_gt_1 d abs_pos
by (intro divide_right_mono) (linarith, simp)
also have "… = abs (x - of_int (p n) / of_int (q n))" using q_gt_1 d(2)
by (simp_all add: divide_simps d(1) mult_ac)
finally show False using approx_int[of n] by simp
qed

text ‹
Next, any irrational algebraic number cannot be approximated with rational
numbers in the Liouville sense.
›
lemma liouville_irrational_algebraic:
fixes x :: real
assumes irrationsl: "x ∉ ℚ" and "algebraic x"
obtains c :: real and n :: nat
where "c > 0" and "⋀(p::int) (q::int). q > 0 ⟹ abs (x - p / q) > c / of_int q ^ n"
proof -
from ‹algebraic x› obtain p where p: "⋀i. coeff p i ∈ ℤ" "p ≠ 0" "poly p x = 0"
by (elim algebraicE) blast
define n where "n = degree p"

― ‹The derivative of @{term p} is bounded within @{term "{x-1..x+1}"}.›
let ?f = "λt. ¦poly (pderiv p) t¦"
define M where "M = (SUP t∈{x-1..x+1}. ?f t)"
define roots where "roots = {x. poly p x = 0} - {x}"

define A_set where "A_set = {1, 1/M} ∪ {abs (x' - x) |x'. x' ∈ roots}"
define A' where "A' = Min A_set"
define A  where "A = A' / 2"
― ‹We define @{term A} to be a positive real number that is less than @{term 1},
@{term "1/M"} and any distance from @{term x} to another root of @{term p}.›

― ‹Properties of @{term M}, our upper bound for the derivative of @{term p}:›
have "∃s∈{x-1..x+1}. ∀y∈{x-1..x+1}. ?f s ≥ ?f y"
by (intro continuous_attains_sup) (auto intro!: continuous_intros)
hence bdd: "bdd_above (?f  {x-1..x+1})" by auto

have M_pos: "M > 0"
proof -
from p have "pderiv p ≠ 0" by (auto dest!: pderiv_iszero)
hence "finite {x. poly (pderiv p) x = 0}" using poly_roots_finite by blast
moreover have "¬finite {x-1..x+1}" by (simp add: infinite_Icc)
ultimately have "¬finite ({x-1..x+1} - {x. poly (pderiv p) x = 0})" by simp
hence "{x-1..x+1} - {x. poly (pderiv p) x = 0} ≠ {}" by (intro notI) simp
then obtain t where t: "t ∈ {x-1..x+1}" and "poly (pderiv p) t ≠ 0" by blast
hence "0 < ?f t" by simp
also from t and bdd have "… ≤ M" unfolding M_def by (rule cSUP_upper)
finally show "M > 0" .
qed

have M_sup: "?f t ≤ M" if "t ∈ {x-1..x+1}" for t
proof -
from that and bdd show "?f t ≤ M"
unfolding M_def by (rule cSUP_upper)
qed

― ‹Properties of @{term A}:›
from p poly_roots_finite[of p] have "finite A_set"
unfolding A_set_def roots_def by simp
have "x ∉ roots" unfolding roots_def by auto
hence "A' > 0" using Min_gr_iff[OF ‹finite A_set›, folded A'_def, of 0]
by (auto simp: A_set_def M_pos)
hence A_pos: "A > 0" by (simp add: A_def)

from ‹A' > 0› have "A < A'" by (simp add: A_def)
moreover from Min_le[OF ‹finite A_set›, folded A'_def]
have "A' ≤ 1" "A' ≤ 1/M" "⋀x'. x' ≠ x ⟹ poly p x' = 0 ⟹ A' ≤ abs (x' - x)"
unfolding A_set_def roots_def by auto
ultimately have A_less: "A < 1" "A < 1/M" "⋀x'. x' ≠ x ⟹ poly p x' = 0 ⟹ A < abs (x' - x)"
by fastforce+

― ‹Finally: no rational number can approximate @{term x} well enough''.›
have "∀(a::int) (b :: int). b > 0 ⟶ ¦x - of_int a / of_int b¦ > A / of_int b ^ n"
proof (intro allI impI, rule ccontr)
fix a b :: int
assume b: "b > 0" and "¬(A / of_int b ^ n < ¦x - of_int a / of_int b¦)"
hence ab: "abs (x - of_int a / of_int b) ≤ A / of_int b ^ n" by simp
also from A_pos b have "A / of_int b ^ n ≤ A / 1"
by (intro divide_left_mono) simp_all
finally have ab': "abs (x - a / b) ≤ A" by simp
also have "… ≤ 1" using A_less by simp
finally have ab'': "a / b ∈ {x-1..x+1}" by auto

have no_root: "poly p (a / b) ≠ 0"
proof
assume "poly p (a / b) = 0"
moreover from ‹x ∉ ℚ› have "x ≠ a / b" by auto
ultimately have "A < abs (a / b - x)" using A_less(3) by simp
with ab' show False by simp
qed

have "∃x0∈{x-1..x+1}. poly p (a / b) - poly p x = (a / b - x) * poly (pderiv p) x0"
using ab'' by (intro poly_MVT') (auto simp: min_def max_def)
with p obtain x0 :: real where x0:
"x0 ∈ {x-1..x+1}" "poly p (a / b) = (a / b - x) * poly (pderiv p) x0" by auto

from x0(2) no_root have deriv_pos: "poly (pderiv p) x0 ≠ 0" by auto

from b p no_root have p_ab: "abs (poly p (a / b)) ≥ 1 / of_int b ^ n"
unfolding n_def by (intro int_poly_rat_no_root_ge)

note ab
also from A_less b have "A / of_int b ^ n < (1/M) / of_int b ^ n"
by (intro divide_strict_right_mono) simp_all
also have "… = (1 / b ^ n) / M" by simp
also from M_pos ab p_ab have "… ≤ abs (poly p (a / b)) / M"
by (intro divide_right_mono) simp_all
also from deriv_pos M_pos x0(1)
have "… ≤ abs (poly p (a / b)) / abs (poly (pderiv p) x0)"
by (intro divide_left_mono M_sup) simp_all
also have "¦poly p (a / b)¦ = ¦a / b - x¦ * ¦poly (pderiv p) x0¦"
by (subst x0(2)) (simp add: abs_mult)
with deriv_pos have "¦poly p (a / b)¦ / ¦poly (pderiv p) x0¦ = ¦x - a / b¦"
by (subst nonzero_divide_eq_eq) simp_all
finally show False by simp
qed
with A_pos show ?thesis using that[of A n] by blast
qed

text ‹
Since Liouville numbers are irrational, but can be approximated well by rational
numbers in the Liouville sense, they must be transcendental.
›
lemma (in liouville) transcendental: "¬algebraic x"
proof
assume "algebraic x"
from liouville_irrational_algebraic[OF irrational this]
obtain c n where cn:
"c > 0" "⋀p q. q > 0 ⟹ c / real_of_int q ^ n < ¦x - real_of_int p / real_of_int q¦"
by auto

define r where "r = nat ⌈log 2 (1 / c)⌉"
define m where "m = n + r"
from cn(1) have "(1/c) = 2 powr log 2 (1/c)" by (subst powr_log_cancel) simp_all
also have "log 2 (1/c) ≤ of_nat (nat ⌈log 2 (1/c)⌉)" by linarith
hence "2 powr log 2 (1/c) ≤ 2 powr …" by (intro powr_mono) simp_all
also have "… = 2 ^ r" unfolding r_def by (simp add: powr_realpow)
finally have "1 / (2^r) ≤ c" using cn(1) by (simp add: field_simps)

have "abs (x - p m / q m) < 1 / of_int (q m) ^ m" by (rule approx_int)
also have "… = (1 / of_int (q m) ^ r) * (1 / real_of_int (q m) ^ n)"
also from denom_gt_1[of m] have "1 / real_of_int (q m) ^ r ≤ 1 / 2 ^ r"
by (intro divide_left_mono power_mono) simp_all
also have "… ≤ c" by fact
finally have "abs (x - p m / q m) < c / of_int (q m) ^ n"
using denom_gt_1[of m] by - (simp_all add: divide_right_mono)
with cn(2)[of "q m" "p m"] denom_gt_1[of m] show False by simp
qed

subsection ‹Standard construction for Liouville numbers›

text ‹
We now define the standard construction for Liouville numbers.
›
definition standard_liouville :: "(nat ⇒ int) ⇒ int ⇒ real" where
"standard_liouville p q = (∑k. p k / of_int q ^ fact (Suc k))"

lemma standard_liouville_summable:
fixes p :: "nat ⇒ int" and q :: int
assumes "q > 1" "range p ⊆ {0..<q}"
shows   "summable (λk. p k / of_int q ^ fact (Suc k))"
proof (rule summable_comparison_test')
from assms(1) show "summable (λn. (1 / q) ^ n)"
by (intro summable_geometric) simp_all
next
fix n :: nat
from assms have "p n ∈ {0..<q}" by blast
with assms(1) have "norm (p n / of_int q ^ fact (Suc n)) ≤
q / of_int q ^ (fact (Suc n))" by (auto simp: field_simps)
also from assms(1) have "… = 1 / of_int q ^ (fact (Suc n) - 1)"
by (subst power_diff) (auto simp del: fact_Suc)
also have "Suc n ≤ fact (Suc n)" by (rule fact_ge_self)
with assms(1) have "1 / real_of_int q ^ (fact (Suc n) - 1) ≤ 1 / of_int q ^ n"
by (intro divide_left_mono power_increasing)
(auto simp del: fact_Suc simp add: algebra_simps)
finally show "norm (p n / of_int q ^ fact (Suc n)) ≤ (1 / q) ^ n"
qed

lemma standard_liouville_sums:
assumes "q > 1" "range p ⊆ {0..<q}"
shows   "(λk. p k / of_int q ^ fact (Suc k)) sums standard_liouville p q"
using standard_liouville_summable[OF assms] unfolding standard_liouville_def

text ‹
Now we prove that the standard construction indeed yields Liouville numbers.
›
lemma standard_liouville_is_liouville:
assumes "q > 1" "range p ⊆ {0..<q}" "frequently (λn. p n ≠ 0) sequentially"
defines "b ≡ λn. q ^ fact (Suc n)"
defines "a ≡ λn. (∑k≤n. p k * q ^ (fact (Suc n) - fact (Suc k)))"
shows   "liouville (standard_liouville p q) a b"
proof
fix n :: nat
from assms(1) have "1 < q ^ 1" by simp
also from assms(1) have "… ≤ b n" unfolding b_def
by (intro power_increasing) (simp_all del: fact_Suc)
finally show "b n > 1" .

note summable = standard_liouville_summable[OF assms(1,2)]
let ?S = "∑k. p (k + n + 1) / of_int q ^ (fact (Suc (k + n + 1)))"
let ?C = "(q - 1) / of_int q ^ (fact (n+2))"

have "a n / b n = (∑k≤n. p k * (of_int q ^ (fact (n+1) - fact (k+1)) / of_int q ^ (fact (n+1))))"
by (simp add: a_def b_def sum_divide_distrib of_int_sum)
also have "… = (∑k≤n. p k / of_int q ^ (fact (Suc k)))"
by (intro sum.cong refl, subst inverse_divide [symmetric], subst power_diff [symmetric])
(insert assms(1), simp_all add: divide_simps fact_mono_nat del: fact_Suc)
also have "standard_liouville p q - … = ?S" unfolding standard_liouville_def
by (subst diff_eq_eq) (intro suminf_split_initial_segment' summable)
finally have "abs (standard_liouville p q - a n / b n) = abs ?S" by (simp only: )
moreover from assms have "?S ≥ 0"
by (intro suminf_nonneg allI divide_nonneg_pos summable_ignore_initial_segment summable) force+
ultimately have eq: "abs (standard_liouville p q - a n / b n) = ?S" by simp

also have "?S ≤ (∑k. ?C * (1 / q) ^ k)"
proof (intro suminf_le allI summable_ignore_initial_segment summable)
from assms show "summable (λk. ?C * (1 / q) ^ k)"
by (intro summable_mult summable_geometric) simp_all
next
fix k :: nat
from assms have "p (k + n + 1) ≤ q - 1" by force
with ‹q > 1› have "p (k + n + 1) / of_int q ^ fact (Suc (k + n + 1)) ≤
(q - 1) / of_int q ^ (fact (Suc (k + n + 1)))"
by (intro divide_right_mono) (simp_all del: fact_Suc)
also from ‹q > 1› have "… ≤ (q - 1) / of_int q ^ (fact (n+2) + k)"
using fact_ineq[of "n+2" k]
also have "… = ?C * (1 / q) ^ k"
finally show "p (k + n + 1) / of_int q ^ fact (Suc (k + n + 1)) ≤ …" .
qed
also from assms have "… = ?C * (∑k. (1 / q) ^ k)"
by (intro suminf_mult summable_geometric) simp_all
also from assms have "(∑k. (1 / q) ^ k) = 1 / (1 - 1 / q)"
by (intro suminf_geometric) simp_all
also from assms(1) have "?C * … = of_int q ^ 1 / of_int q ^ fact (n + 2)"
also from assms(1) have "… ≤ of_int q ^ fact (n + 1) / of_int q ^ fact (n + 2)"
by (intro divide_right_mono power_increasing) (simp_all add: field_simps del: fact_Suc)
also from assms(1) have "… = 1 / (of_int q ^ (fact (n + 2) - fact (n + 1)))"
by (subst power_diff) simp_all
also have "fact (n + 2) - fact (n + 1) = (n + 1) * fact (n + 1)" by (simp add: algebra_simps)
also from assms(1) have "1 / (of_int q ^ …) < 1 / (real_of_int q ^ (fact (n + 1) * n))"
by (intro divide_strict_left_mono power_increasing mult_right_mono) simp_all
also have "… = 1 / of_int (b n) ^ n"
by (simp add: power_mult b_def power_divide del: fact_Suc)
finally show "¦standard_liouville p q - a n / b n¦ < 1 / of_int (b n) ^ n" .

from assms(3) obtain k where k: "k ≥ n + 1" "p k ≠ 0"
by (auto simp: frequently_def eventually_at_top_linorder)
define k' where "k' = k - n - 1"
from assms k have "p k ≥ 0" by force
with k assms have k': "p (k' + n + 1) > 0" unfolding k'_def by force
with assms(1,2) have "?S > 0"
by (intro suminf_pos2[of _ k'] summable_ignore_initial_segment summable)
(force intro!: divide_pos_pos divide_nonneg_pos)+
with eq show "¦standard_liouville p q - a n / b n¦ > 0" by (simp only: )
qed

text ‹
We can now show our main result: any standard Liouville number is transcendental.
›
theorem transcendental_standard_liouville:
assumes "q > 1" "range p ⊆ {0..<q}" "frequently (λk. p k ≠ 0) sequentially"
shows   "¬algebraic (standard_liouville p q)"
proof -
from assms interpret
liouville "standard_liouville p q"
"λn. (∑k≤n. p k * q ^ (fact (Suc n) - fact (Suc k)))"
"λn. q ^ fact (Suc n)"
by (rule standard_liouville_is_liouville)
from transcendental show ?thesis .
qed

text ‹
In particular: The the standard construction for constant sequences, such as the
classic'' Liouville constant $\sum_{n=1}^\infty 10^{-n!} = 0.11000100\ldots$,
are transcendental.

This shows that Liouville numbers exists and therefore gives a concrete and
elementary proof that transcendental numbers exist.
›
corollary transcendental_standard_standard_liouville:
"a ∈ {0<..<b} ⟹ ¬algebraic (standard_liouville (λ_. int a) (int b))"
by (intro transcendental_standard_liouville) auto

corollary transcendental_liouville_constant:
"¬algebraic (standard_liouville (λ_. 1) 10)"
by (intro transcendental_standard_liouville) auto

end
`