Theory Pell
section ‹Pell's equation›
theory Pell
imports
Complex_Main
"HOL-Computational_Algebra.Computational_Algebra"
begin
text ‹
Pell's equation has the general form $x^2 = 1 + D y^2$ where ‹D ∈ ℕ› is a parameter
and ‹x›, ‹y› are ‹ℤ›-valued variables. As we will see, that case where ‹D› is a
perfect square is trivial and therefore uninteresting; we will therefore assume that ‹D› is
not a perfect square for the most part.
Furthermore, it is obvious that the solutions to Pell's equation are symmetric around the
origin in the sense that ‹(x, y)› is a solution iff ‹(±x, ±y)› is a solution. We will
therefore mostly look at solutions ‹(x, y)› where both ‹x› and ‹y› are non-negative, since
the remaining solutions are a trivial consequence of these.
Information on the material treated in this formalisation can be found in many textbooks and
lecture notes, e.\,g.\ \<^cite>‹"jacobson2008solving" and "auckland_pell"›.
›
subsection ‹Preliminary facts›
lemma gcd_int_nonpos_iff [simp]: "gcd x (y :: int) ≤ 0 ⟷ x = 0 ∧ y = 0"
proof
assume "gcd x y ≤ 0"
with gcd_ge_0_int[of x y] have "gcd x y = 0" by linarith
thus "x = 0 ∧ y = 0" by auto
qed auto
lemma minus_in_Ints_iff [simp]:
"-x ∈ ℤ ⟷ x ∈ ℤ"
using Ints_minus[of x] Ints_minus[of "-x"] by auto
text ‹
A (positive) square root of a natural number is either a natural number or irrational.
›
lemma nonneg_sqrt_nat_or_irrat:
assumes "x ^ 2 = real a" and "x ≥ 0"
shows "x ∈ ℕ ∨ x ∉ ℚ"
proof safe
assume "x ∉ ℕ" and "x ∈ ℚ"
from Rats_abs_nat_div_natE[OF this(2)]
obtain p q :: nat where q_nz [simp]: "q ≠ 0" and "abs x = p / q" and coprime: "coprime p q" .
with ‹x ≥ 0› have x: "x = p / q"
by simp
with assms have "real (q ^ 2) * real a = real (p ^ 2)"
by (simp add: field_simps)
also have "real (q ^ 2) * real a = real (q ^ 2 * a)"
by simp
finally have "p ^ 2 = q ^ 2 * a"
by (subst (asm) of_nat_eq_iff) auto
hence "q ^ 2 dvd p ^ 2"
by simp
hence "q dvd p"
by simp
with coprime have "q = 1"
by auto
with x and ‹x ∉ ℕ› show False
by simp
qed
text ‹
A square root of a natural number is either an integer or irrational.
›
corollary sqrt_nat_or_irrat:
assumes "x ^ 2 = real a"
shows "x ∈ ℤ ∨ x ∉ ℚ"
proof (cases "x ≥ 0")
case True
with nonneg_sqrt_nat_or_irrat[OF assms this]
show ?thesis by (auto simp: Nats_altdef2)
next
case False
from assms have "(-x) ^ 2 = real a"
by simp
moreover from False have "-x ≥ 0"
by simp
ultimately have "-x ∈ ℕ ∨ -x ∉ ℚ"
by (rule nonneg_sqrt_nat_or_irrat)
thus ?thesis
by (auto simp: Nats_altdef2)
qed
corollary sqrt_nat_or_irrat':
"sqrt (real a) ∈ ℕ ∨ sqrt (real a) ∉ ℚ"
using nonneg_sqrt_nat_or_irrat[of "sqrt a" a] by auto
text ‹
The square root of a natural number ‹n› is again a natural number iff ‹n is a perfect square.›
›
corollary sqrt_nat_iff_is_square:
"sqrt (real n) ∈ ℕ ⟷ is_square n"
proof
assume "sqrt (real n) ∈ ℕ"
then obtain k where "sqrt (real n) = real k" by (auto elim!: Nats_cases)
hence "sqrt (real n) ^ 2 = real (k ^ 2)" by (simp only: of_nat_power)
also have "sqrt (real n) ^ 2 = real n" by simp
finally have "n = k ^ 2" by (simp only: of_nat_eq_iff)
thus "is_square n" by blast
qed (auto elim!: is_nth_powerE)
corollary irrat_sqrt_nonsquare: "¬is_square n ⟹ sqrt (real n) ∉ ℚ"
using sqrt_nat_or_irrat'[of n] by (auto simp: sqrt_nat_iff_is_square)
subsection ‹The case of a perfect square›
text ‹
As we have noted, the case where ‹D› is a perfect square is trivial: In fact, we will
show that the only solutions in this case are the trivial solutions ‹(x, y) = (±1, 0)› if
‹D› is a non-zero perfect square, or ‹(±1, y)› for arbitrary ‹y ∈ ℤ› if ‹D = 0›.
›
context
fixes D :: nat
assumes square_D: "is_square D"
begin
lemma pell_square_solution_nat_aux:
fixes x y :: nat
assumes "D > 0" and "x ^ 2 = 1 + D * y ^ 2"
shows "(x, y) = (1, 0)"
proof -
from assms have x_nz: "x > 0" by (auto intro!: Nat.gr0I)
from square_D obtain d where [simp]: "D = d⇧2"
by (auto elim: is_nth_powerE)
have "int x ^ 2 = int (x ^ 2)" by simp
also note assms(2)
also have "int (1 + D * y ^ 2) = 1 + int D * int y ^ 2" by simp
finally have "(int x + int d * int y) * (int x - int d * int y) = 1"
by (simp add: algebra_simps power2_eq_square)
hence *: "int x + int d * int y = 1 ∧ int x - int d * int y = 1"
using x_nz by (subst (asm) pos_zmult_eq_1_iff) (auto intro: add_pos_nonneg)
from * have [simp]: "x = 1" by simp
moreover from * and assms(1) have "y = 0" by auto
ultimately show ?thesis by simp
qed
lemma pell_square_solution_int_aux:
fixes x y :: int
assumes "D > 0" and "x ^ 2 = 1 + D * y ^ 2"
shows "x ∈ {-1, 1} ∧ y = 0"
proof -
define x' y' where "x' = nat ¦x¦" and "y' = nat ¦y¦"
have x: "x = sgn x * x'" and y: "y = sgn y * y'"
by (auto simp: sgn_if x'_def y'_def)
have zero_iff: "x = 0 ⟷ x' = 0" "y = 0 ⟷ y' = 0"
by (auto simp: x'_def y'_def)
note assms(2)
also have "x ^ 2 = int (x' ^ 2)"
by (subst x) (auto simp: sgn_if zero_iff)
also have "1 + D * y ^ 2 = int (1 + D * y' ^ 2)"
by (subst y) (auto simp: sgn_if zero_iff)
also note of_nat_eq_iff
finally have "x'⇧2 = 1 + D * y'⇧2" .
from ‹D > 0› and this have "(x', y') = (1, 0)"
by (rule pell_square_solution_nat_aux)
thus ?thesis by (auto simp: x'_def y'_def)
qed
lemma pell_square_solution_nat_iff:
fixes x y :: nat
shows "x ^ 2 = 1 + D * y ^ 2 ⟷ x = 1 ∧ (D = 0 ∨ y = 0)"
using pell_square_solution_nat_aux[of x y] by (cases "D = 0") auto
lemma pell_square_solution_int_iff:
fixes x y :: int
shows "x ^ 2 = 1 + D * y ^ 2 ⟷ x ∈ {-1, 1} ∧ (D = 0 ∨ y = 0)"
using pell_square_solution_int_aux[of x y] by (cases "D = 0") (auto simp: power2_eq_1_iff)
end
subsection ‹Existence of a non-trivial solution›
text ‹
Let us now turn to the case where ‹D› is not a perfect square.
We first show that Pell's equation always has at least one non-trivial solution (apart
from the trivial solution ‹(1, 0)›). For this, we first need a lemma about the existence
of rational approximations of real numbers.
The following lemma states that for any positive integer ‹s› and real number ‹x›, we can find a
rational approximation ‹t / u› to ‹x› with an error of most ‹1 / (u * s)› where the denominator
‹u› is at most ‹s›.
›
lemma pell_approximation_lemma:
fixes s :: nat and x :: real
assumes s: "s > 0"
shows "∃u::nat. ∃t::int. u > 0 ∧ coprime u t ∧ 1 / s ∈ {¦t - u * x¦<..1 / u}"
proof -
define f where "f = (λu. ⌈u * x⌉)"
define g :: "nat ⇒ int" where "g = (λu. ⌊frac (u * x) * s⌋)"
{
fix u :: nat assume u: "u ≤ s"
hence "frac (u * x) * real s < 1 * real s"
using s by (intro mult_strict_right_mono) (auto simp: frac_lt_1)
hence "g u < int s" by (auto simp: floor_less_iff g_def)
}
hence "g ` {..s} ⊆ {0..<s}"
by (auto simp: g_def floor_less_iff)
hence "card (g ` {..s}) ≤ card {0..<int s}"
by (intro card_mono) auto
also have "… < card {..s}" by simp
finally have "¬inj_on g {..s}" by (rule pigeonhole)
then obtain a b where ab: "a ≤ s" "b ≤ s" "a ≠ b" "g a = g b"
by (auto simp: inj_on_def)
define u1 and u2 where "u1 = max a b" and "u2 = min a b"
have u12: "u1 ≤ s" "u2 ≤ s" "u2 < u1" "g u1 = g u2"
using ab by (auto simp: u1_def u2_def)
define u t where "u = u1 - u2" and "t = ⌊u1 * x⌋ - ⌊u2 * x⌋"
have u: "u > 0" "¦u¦ ≤ s"
using u12 by (simp_all add: u_def)
from ‹g u1 = g u2› have "¦frac (u2 * x) * s - frac (u1 * x) * s¦ < 1"
unfolding g_def by linarith
also have "¦frac (u2 * x) * s - frac (u1 * x) * s¦ =
¦real s¦ * ¦frac (u2 * x) - frac (u1 * x)¦"
by (subst abs_mult [symmetric]) (simp add: algebra_simps)
finally have "¦t - u * x¦ * s < 1" using ‹u1 > u2›
by (simp add: g_def u_def t_def frac_def algebra_simps of_nat_diff)
with ‹s > 0› have less: "¦t - u * x¦ < 1 / s" by (simp add: divide_simps)
define d where "d = gcd (nat ¦t¦) u"
define t' :: int and u' :: nat where "t'= t div d" and "u' = u div d"
from u have "d ≠ 0"
by (intro notI) (auto simp: d_def)
have "int (gcd (nat ¦t¦) u) = gcd ¦t¦ (int u)"
by simp
hence t'_u': "t = t' * d" "u = u' * d"
by (auto simp: t'_def u'_def d_def nat_dvd_iff)
from ‹d ≠ 0› have "¦t' - u' * x¦ * 1 ≤ ¦t' - u' * x¦ * ¦real d¦"
by (intro mult_left_mono) auto
also have "… = ¦t - u * x¦" by (subst abs_mult [symmetric]) (simp add: algebra_simps t'_u')
also note less
finally have "¦t' - u' * x¦ < 1 / s" by simp
moreover {
from ‹s > 0› and u have "1 / s ≤ 1 / u"
by (simp add: divide_simps u_def)
also have "… = 1 / u' / d" by (simp add: t'_u' divide_simps)
also have "… ≤ 1 / u' / 1" using ‹d ≠ 0› by (intro divide_left_mono) auto
finally have "1 / s ≤ 1 / u'" by simp
}
ultimately have "1 / s ∈ {¦t' - u' * x¦<..1 / u'}" by auto
moreover from ‹u > 0› have "u' > 0" by (auto simp: t'_u')
moreover {
have "gcd u t = gcd t' u' * int d"
by (simp add: t'_u' gcd_mult_right gcd.commute)
also have "int d = gcd u t"
by (simp add: d_def gcd.commute)
finally have "gcd u' t' = 1" using u by (simp add: gcd.commute)
}
ultimately show ?thesis by blast
qed
text ‹
As a simple corollary of this, we can show that for irrational ‹x›, there is an infinite
number of rational approximations ‹t / u› to ‹x› whose error is less that ‹1 / u⇧2›.
›
corollary pell_approximation_corollary:
fixes x :: real
assumes "x ∉ ℚ"
shows "infinite {(t :: int, u :: nat). u > 0 ∧ coprime u t ∧ ¦t - u * x¦ < 1 / u}"
(is "infinite ?A")
proof
assume fin: "finite ?A"
let ?f = "λ(t :: int, u :: nat). ¦t - u * x¦"
from fin have fin': "finite (insert 1 (?f ` ?A))" by blast
have "Min (insert 1 (?f ` ?A)) > 0"
proof (subst Min_gr_iff)
have "a ≠ b * x" if "b > 0" for a :: int and b :: nat
proof
assume "a = b * x"
with ‹b > 0› have "x = a / b" by (simp add: field_simps)
with ‹x ∉ ℚ› and ‹b > 0› show False by (auto simp: Rats_eq_int_div_nat)
qed
thus "∀x∈insert 1 (?f ` ?A). x > 0" by auto
qed (insert fin', simp_all)
also note real_arch_inverse
finally obtain M :: nat where M: "M ≠ 0" "inverse M < Min (insert 1 (?f ` ?A))"
by blast
hence "M > 0" by simp
from pell_approximation_lemma[OF this, of x] obtain u :: nat and t :: int
where ut: "u > 0" "coprime u t" "1 / real M ∈ {?f (t, u)<..1 / u}" by auto
from ut have "?f (t, u) < 1 / real M" by simp
also from M have "… < Min (insert 1 (?f ` ?A))"
by (simp add: divide_simps)
also from ut have "Min (insert 1 (?f ` ?A)) ≤ ?f (t, u)"
by (intro Min.coboundedI fin') auto
finally show False by simp
qed
locale pell =
fixes D :: nat
assumes nonsquare_D: "¬is_square D"
begin
lemma D_gt_1: "D > 1"
proof -
from nonsquare_D have "D ≠ 0" "D ≠ 1" by (auto intro!: Nat.gr0I)
thus ?thesis by simp
qed
lemma D_pos: "D > 0"
using nonsquare_D by (intro Nat.gr0I) auto
text ‹
With the above corollary, we can show the existence of a non-trivial solution. We restrict our
attention to solutions ‹(x, y)› where both ‹x› and ‹y› are non-negative.
›
theorem pell_solution_exists: "∃(x::nat) (y::nat). y ≠ 0 ∧ x⇧2 = 1 + D * y⇧2"
proof -
define S where "S = {(t :: int, u :: nat). u > 0 ∧ coprime u t ∧ ¦t - u * sqrt D¦ < 1 / u}"
let ?f = "λ(t :: int, u :: nat). t⇧2 - u⇧2 * D"
define M where "M = ⌊1 + 2 * sqrt D⌋"
have infinite: "¬finite S" unfolding S_def
by (intro pell_approximation_corollary irrat_sqrt_nonsquare nonsquare_D)
have subset: "?f ` S ⊆ {-M..M}"
proof safe
fix u :: nat and t :: int
assume tu: "(t, u) ∈ S"
from tu have [simp]: "u > 0" by (auto simp: S_def)
have "¦t + u * sqrt D¦ = ¦t - u * sqrt D + 2 * u * sqrt D¦" by simp
also have "… ≤ ¦t - u * sqrt D¦ + ¦2 * u * sqrt D¦"
by (rule abs_triangle_ineq)
also have "¦2 * u * sqrt D¦ = 2 * u * sqrt D" by simp
also have "¦t - u * sqrt D¦ ≤ 1 / u"
using tu by (simp add: S_def)
finally have le: "¦t + u * sqrt D¦ ≤ 1 / u + 2 * u * sqrt D" by simp
have "¦t⇧2 - u⇧2 * D¦ = ¦t - u * sqrt D¦ * ¦t + u * sqrt D¦"
by (subst abs_mult [symmetric]) (simp add: algebra_simps power2_eq_square)
also have "… ≤ 1 / u * (1 / u + 2 * u * sqrt D)"
using tu by (intro mult_mono le) (auto simp: S_def)
also have "… = 1 / real u ^ 2 + 2 * sqrt D"
by (simp add: algebra_simps power2_eq_square)
also from ‹u > 0› have "real u ≥ 1" by linarith
hence "1 / real u ^ 2 ≤ 1 / 1 ^ 2"
by (intro divide_left_mono power_mono) auto
finally have "¦t⇧2 - u⇧2 * D¦ ≤ 1 + 2 * sqrt D" by simp
hence "t⇧2 - u⇧2 * D ≥ -M" "t⇧2 - u⇧2 * D ≤ M" unfolding M_def by linarith+
thus "t⇧2 - u⇧2 * D ∈ {-M..M}" by simp
qed
hence fin: "finite (?f ` S)" by (rule finite_subset) auto
from pigeonhole_infinite[OF infinite fin]
obtain z where z: "z ∈ S" "infinite {z' ∈ S. ?f z' = ?f z}" by blast
define k where "k = ?f z"
with subset and z have k: "k ∈ {-M..M}" "infinite {z∈S. ?f z = k}"
by (auto simp: k_def)
have k_nz: "k ≠ 0"
proof
assume [simp]: "k = 0"
note k(2)
also have "?f z ≠ 0" if "z ∈ S" for z
proof
assume *: "?f z = 0"
obtain t u where [simp]: "z = (t, u)" by (cases z)
from * have "t ^ 2 = int u ^ 2 * int D" by simp
hence "int u ^ 2 dvd t ^ 2" by simp
hence "int u dvd t" by simp
then obtain k where [simp]: "t = int u * k" by (auto elim!: dvdE)
from * and ‹z ∈ S› have "k ^ 2 = int D"
by (auto simp: power_mult_distrib S_def)
also have "k ^ 2 = int (nat ¦k¦ ^ 2)" by simp
finally have "D = nat ¦k¦ ^ 2" by (simp only: of_nat_eq_iff)
hence "is_square D" by auto
with nonsquare_D show False by contradiction
qed
hence "{z∈S. ?f z = k} = {}" by auto
finally show False by simp
qed
let ?h = "λ(t :: int, u :: nat). (t mod (abs k), u mod (abs k))"
have "?h ` {z∈S. ?f z = k} ⊆ {0..<abs k} × {0..< abs k}"
using k_nz by (auto simp: case_prod_unfold)
hence "finite (?h ` {z∈S. ?f z = k})" by (rule finite_subset) auto
from pigeonhole_infinite[OF k(2) this] obtain z'
where z': "z' ∈ S" "?f z' = k" "infinite {z''∈{z∈S. ?f z = k}. ?h z'' = ?h z'}"
by blast
define l1 and l2 where "l1 = fst (?h z')" and "l2 = snd (?h z')"
define S' where "S' = {(t,u) ∈ S. ?f (t,u) = k ∧ t mod abs k = l1 ∧ u mod abs k = l2}"
note z'(3)
also have "{z''∈{z∈S. ?f z = k}. ?h z'' = ?h z'} = S'"
by (auto simp: l1_def l2_def case_prod_unfold S'_def)
finally have infinite: "infinite S'" .
from z'(1) and k_nz have l12: "l1 ∈ {0..<abs k}" "l2 ∈ {0..<abs k}"
by (auto simp: l1_def l2_def case_prod_unfold)
from infinite_arbitrarily_large[OF infinite]
obtain X where X: "finite X" "card X = 2" "X ⊆ S'" by blast
from finite_distinct_list[OF this(1)] obtain xs where xs: "set xs = X" "distinct xs"
by blast
with X have "length xs = 2" using distinct_card[of xs] by simp
then obtain z1 z2 where [simp]: "xs = [z1, z2]"
by (auto simp: length_Suc_conv eval_nat_numeral)
from X xs have S': "z1 ∈ S'" "z2 ∈ S'" and neq: "z1 ≠ z2" by auto
define t1 u1 t2 u2 where "t1 = fst z1" and "u1 = snd z1" and "t2 = fst z2" and "u2 = snd z2"
have [simp]: "z1 = (t1, u1)" "z2 = (t2, u2)"
by (simp_all add: t1_def u1_def t2_def u2_def)
from S' have * [simp]: "t1 mod abs k = l1" "t2 mod abs k = l1" "u1 mod abs k = l2" "u2 mod abs k = l2"
by (simp_all add: S'_def)
define x where "x = (t1 * t2 - D * u1 * u2) div k"
define y where "y = (t1 * u2 - t2 * u1) div k"
from S' have "(t1⇧2 - u1⇧2 * D) = k" "(t2⇧2 - u2⇧2 * D) = k"
by (auto simp: S'_def)
hence "(t1⇧2 - u1⇧2 * D) * (t2⇧2 - u2⇧2 * D) = k ^ 2"
unfolding power2_eq_square by simp
also have "(t1⇧2 - u1⇧2 * D) * (t2⇧2 - u2⇧2 * D) =
(t1 * t2 - D * u1 * u2) ^ 2 - D * (t1 * u2 - t2 * u1) ^ 2"
by (simp add: power2_eq_square algebra_simps)
finally have eq: "(t1 * t2 - D * u1 * u2)⇧2 - D * (t1 * u2 - t2 * u1)⇧2 = k⇧2" .
have "(t1 * u2 - t2 * u1) mod abs k = (l1 * l2 - l1 * l2) mod abs k"
using l12 by (intro mod_diff_cong mod_mult_cong) (auto simp: mod_pos_pos_trivial)
hence dvd1: "k dvd t1 * u2 - t2 * u1" by (simp add: mod_eq_0_iff_dvd)
have "k⇧2 dvd k⇧2 + D * (t1 * u2 - t2 * u1)⇧2"
using dvd1 by (intro dvd_add) auto
also from eq have "… = (t1 * t2 - D * u1 * u2)⇧2"
by (simp add: algebra_simps)
finally have dvd2: "k dvd t1 * t2 - D * u1 * u2"
by simp
note eq
also from dvd2 have "t1 * t2 - D * u1 * u2 = k * x"
by (simp add: x_def)
also from dvd1 have "t1 * u2 - t2 * u1 = k * y"
by (simp add: y_def)
also have "(k * x)⇧2 - D * (k * y)⇧2 = k⇧2 * (x⇧2 - D * y⇧2)"
by (simp add: power_mult_distrib algebra_simps)
finally have eq': "x⇧2 - D * y⇧2 = 1"
using k_nz by simp
hence "x⇧2 = 1 + D * y⇧2" by simp
also have "x⇧2 = int (nat ¦x¦ ^ 2)" by simp
also have "1 + D * y⇧2 = int (1 + D * nat ¦y¦ ^ 2)" by simp
also note of_nat_eq_iff
finally have eq'': "(nat ¦x¦)⇧2 = 1 + D * (nat ¦y¦)⇧2" .
have "t1 * u2 ≠ t2 * u1"
proof
assume *: "t1 * u2 = t2 * u1"
hence "¦t1¦ * ¦u2¦ = ¦t2¦ * ¦u1¦" by (simp only: abs_mult [symmetric])
moreover from S' have "coprime u1 t1" "coprime u2 t2"
by (auto simp: S'_def S_def)
ultimately have eq: "¦t1¦ = ¦t2¦ ∧ u1 = u2"
by (subst (asm) coprime_crossproduct_int) (auto simp: S'_def S_def gcd.commute coprime_commute)
moreover from S' have "u1 ≠ 0" "u2 ≠ 0" by (auto simp: S'_def S_def)
ultimately have "t1 = t2" using * by auto
with eq and neq show False by auto
qed
with dvd1 have "y ≠ 0"
by (auto simp add: y_def dvd_div_eq_0_iff)
hence "nat ¦y¦ ≠ 0" by auto
with eq'' show "∃x y. y ≠ 0 ∧ x⇧2 = 1 + D * y⇧2" by blast
qed
subsection ‹Definition of solutions›
text ‹
We define some abbreviations for the concepts of a solution and a non-trivial solution.
›
definition solution :: "('a × 'a :: comm_semiring_1) ⇒ bool" where
"solution = (λ(a, b). a⇧2 = 1 + of_nat D * b⇧2)"
definition nontriv_solution :: "('a × 'a :: comm_semiring_1) ⇒ bool" where
"nontriv_solution = (λ(a, b). (a, b) ≠ (1, 0) ∧ a⇧2 = 1 + of_nat D * b⇧2)"
lemma nontriv_solution_altdef: "nontriv_solution z ⟷ solution z ∧ z ≠ (1, 0)"
by (auto simp: solution_def nontriv_solution_def)
lemma solution_trivial_nat [simp, intro]: "solution (Suc 0, 0)"
by (simp add: solution_def)
lemma solution_trivial [simp, intro]: "solution (1, 0)"
by (simp add: solution_def)
lemma solution_uminus_left [simp]: "solution (-x, y :: 'a :: comm_ring_1) ⟷ solution (x, y)"
by (simp add: solution_def)
lemma solution_uminus_right [simp]: "solution (x, -y :: 'a :: comm_ring_1) ⟷ solution (x, y)"
by (simp add: solution_def)
lemma solution_0_snd_nat_iff [simp]: "solution (a :: nat, 0) ⟷ a = 1"
by (auto simp: solution_def)
lemma solution_0_snd_iff [simp]: "solution (a :: 'a :: idom, 0) ⟷ a ∈ {1, -1}"
by (auto simp: solution_def power2_eq_1_iff)
lemma no_solution_0_fst_nat [simp]: "¬solution (0, b :: nat)"
by (auto simp: solution_def)
lemma no_solution_0_fst_int [simp]: "¬solution (0, b :: int)"
proof -
have "1 + int D * b⇧2 > 0" by (intro add_pos_nonneg) auto
thus ?thesis by (auto simp add: solution_def)
qed
lemma solution_of_nat_of_nat [simp]:
"solution (of_nat a, of_nat b :: 'a :: {comm_ring_1, ring_char_0}) ⟷ solution (a, b)"
by (simp only: solution_def prod.case of_nat_power [symmetric]
of_nat_1 [symmetric, where ?'a = 'a] of_nat_add [symmetric]
of_nat_mult [symmetric] of_nat_eq_iff of_nat_id)
lemma solution_of_nat_of_nat' [simp]:
"solution (case z of (a, b) ⇒ (of_nat a, of_nat b :: 'a :: {comm_ring_1, ring_char_0})) ⟷
solution z"
by (auto simp: case_prod_unfold)
lemma solution_nat_abs_nat_abs [simp]:
"solution (nat ¦x¦, nat ¦y¦) ⟷ solution (x, y)"
proof -
define x' and y' where "x' = nat ¦x¦" and "y' = nat ¦y¦"
have x: "x = sgn x * x'" and y: "y = sgn y * y'"
by (auto simp: x'_def y'_def sgn_if)
have [simp]: "x = 0 ⟷ x' = 0" "y = 0 ⟷ y' = 0"
by (auto simp: x'_def y'_def)
show "solution (x', y') ⟷ solution (x, y)"
by (subst x, subst y) (auto simp: sgn_if)
qed
lemma nontriv_solution_of_nat_of_nat [simp]:
"nontriv_solution (of_nat a, of_nat b :: 'a :: {comm_ring_1, ring_char_0}) ⟷ nontriv_solution (a, b)"
by (auto simp: nontriv_solution_altdef)
lemma nontriv_solution_of_nat_of_nat' [simp]:
"nontriv_solution (case z of (a, b) ⇒ (of_nat a, of_nat b :: 'a :: {comm_ring_1, ring_char_0})) ⟷
nontriv_solution z"
by (auto simp: case_prod_unfold)
lemma nontriv_solution_imp_solution [dest]: "nontriv_solution z ⟹ solution z"
by (auto simp: nontriv_solution_altdef)
subsection ‹The Pell valuation function›
text ‹
Solutions ‹(x,y)› have an interesting correspondence to the ring $\mathbb{Z}[\sqrt{D}]$ via
the map $(x,y) \mapsto x + y \sqrt{D}$. We call this map the ∗‹Pell valuation function›.
It is obvious that this map is injective, since $\sqrt{D}$ is irrational.
›
definition pell_valuation :: "int × int ⇒ real" where
"pell_valuation = (λ(a,b). a + b * sqrt D)"
lemma pell_valuation_nonneg [simp]: "fst z ≥ 0 ⟹ snd z ≥ 0 ⟹ pell_valuation z ≥ 0"
by (auto simp: pell_valuation_def case_prod_unfold)
lemma pell_valuation_uminus_uminus [simp]: "pell_valuation (-x, -y) = -pell_valuation (x, y)"
by (simp add: pell_valuation_def)
lemma pell_valuation_eq_iff [simp]:
"pell_valuation z1 = pell_valuation z2 ⟷ z1 = z2"
proof
assume *: "pell_valuation z1 = pell_valuation z2"
obtain a b where [simp]: "z1 = (a, b)" by (cases z1)
obtain u v where [simp]: "z2 = (u, v)" by (cases z2)
have "b = v"
proof (rule ccontr)
assume "b ≠ v"
with * have "sqrt D = (u - a) / (b - v)"
by (simp add: field_simps pell_valuation_def)
also have "… ∈ ℚ" by auto
finally show False using irrat_sqrt_nonsquare nonsquare_D by blast
qed
moreover from this and * have "a = u"
by (simp add: pell_valuation_def)
ultimately show "z1 = z2" by simp
qed auto
subsection ‹Linear ordering of solutions›
text ‹
Next, we show that solutions are linearly ordered w.\,r.\,t.\ the pointwise order on products.
This means thatfor two different solutions ‹(a, b)› and ‹(x, y)›, we always either have
‹a < x› and ‹b < y› or ‹a > x› and ‹b > y›.
›
lemma solutions_linorder:
fixes a b x y :: nat
assumes "solution (a, b)" "solution (x, y)"
shows "a ≤ x ∧ b ≤ y ∨ a ≥ x ∧ b ≥ y"
proof -
have "b ≤ y" if "a ≤ x" "solution (a, b)" "solution (x, y)" for a b x y :: nat
proof -
from that have "a ^ 2 ≤ x ^ 2" by (intro power_mono) auto
with that and D_gt_1 have "b⇧2 ≤ y⇧2"
by (simp add: solution_def)
thus "b ≤ y"
by (simp add: power2_nat_le_eq_le)
qed
from this[of a x b y] and this[of x a y b] and assms show ?thesis
by (cases "a ≤ x") auto
qed
lemma solutions_linorder_strict:
fixes a b x y :: nat
assumes "solution (a, b)" "solution (x, y)"
shows "(a, b) = (x, y) ∨ a < x ∧ b < y ∨ a > x ∧ b > y"
proof -
have "b = y" if "a = x"
using that assms and D_gt_1 by (simp add: solution_def)
moreover have "a = x" if "b = y"
proof -
from that and assms have "a⇧2 = Suc (D * y⇧2)"
by (simp add: solution_def)
also from that and assms have "… = x⇧2"
by (simp add: solution_def)
finally show "a = x" by simp
qed
ultimately have [simp]: "a = x ⟷ b = y" ..
show ?thesis using solutions_linorder[OF assms]
by (cases a x rule: linorder_cases; cases b y rule: linorder_cases) simp_all
qed
lemma solutions_le_iff_pell_valuation_le:
fixes a b x y :: nat
assumes "solution (a, b)" "solution (x, y)"
shows "a ≤ x ∧ b ≤ y ⟷ pell_valuation (a, b) ≤ pell_valuation (x, y)"
proof
assume "a ≤ x ∧ b ≤ y"
thus "pell_valuation (a, b) ≤ pell_valuation (x, y)"
unfolding pell_valuation_def prod.case using D_gt_1
by (intro add_mono mult_right_mono) auto
next
assume *: "pell_valuation (a, b) ≤ pell_valuation (x, y)"
from assms have "a ≤ x ∧ b ≤ y ∨ x ≤ a ∧ y ≤ b"
by (rule solutions_linorder)
thus "a ≤ x ∧ b ≤ y"
proof
assume "x ≤ a ∧ y ≤ b"
hence "pell_valuation (a, b) ≥ pell_valuation (x, y)"
unfolding pell_valuation_def prod.case using D_gt_1
by (intro add_mono mult_right_mono) auto
with * have "pell_valuation (a, b) = pell_valuation (x, y)" by linarith
hence "(a, b) = (x, y)" by simp
thus "a ≤ x ∧ b ≤ y" by simp
qed auto
qed
lemma solutions_less_iff_pell_valuation_less:
fixes a b x y :: nat
assumes "solution (a, b)" "solution (x, y)"
shows "a < x ∧ b < y ⟷ pell_valuation (a, b) < pell_valuation (x, y)"
proof
assume "a < x ∧ b < y"
thus "pell_valuation (a, b) < pell_valuation (x, y)"
unfolding pell_valuation_def prod.case using D_gt_1
by (intro add_strict_mono mult_strict_right_mono) auto
next
assume *: "pell_valuation (a, b) < pell_valuation (x, y)"
from assms have "(a, b) = (x, y) ∨ a < x ∧ b < y ∨ x < a ∧ y < b"
by (rule solutions_linorder_strict)
thus "a < x ∧ b < y"
proof (elim disjE)
assume "x < a ∧ y < b"
hence "pell_valuation (a, b) > pell_valuation (x, y)"
unfolding pell_valuation_def prod.case using D_gt_1
by (intro add_strict_mono mult_strict_right_mono) auto
with * have False by linarith
thus ?thesis ..
qed (insert *, auto)
qed
subsection ‹The fundamental solution›
text ‹
The ∗‹fundamental solution› is the non-trivial solution ‹(x, y)› with non-negative ‹x› and ‹y›
for which the Pell valuation $x + y\sqrt{D}$ is minimal, or, equivalently, for which ‹x› and ‹y›
are minimal.
›
definition fund_sol :: "nat × nat" where
"fund_sol = (THE z::nat×nat. is_arg_min (pell_valuation :: nat × nat ⇒ real) nontriv_solution z)"
text ‹
The well-definedness of this follows from the injectivity of the Pell valuation and the fact
that smaller Pell valuation of a solution is smaller than that of another iff the components
are both smaller.
›
theorem fund_sol_is_arg_min:
"is_arg_min (pell_valuation :: nat × nat ⇒ real) nontriv_solution fund_sol"
unfolding fund_sol_def
proof (rule theI')
show "∃!z::nat×nat. is_arg_min (pell_valuation :: nat × nat ⇒ real) nontriv_solution z"
proof (rule ex_ex1I)
fix z1 z2 :: "nat × nat"
assume "is_arg_min (pell_valuation :: nat × nat ⇒ real) nontriv_solution z1"
"is_arg_min (pell_valuation :: nat × nat ⇒ real) nontriv_solution z2"
hence "pell_valuation z1 = pell_valuation z2"
by (cases z1, cases z2, intro antisym) (auto simp: is_arg_min_def not_less)
thus "z1 = z2" by (auto split: prod.splits)
next
define y where "y = (LEAST y. y > 0 ∧ is_square (1 + D * y⇧2))"
have "∃y>0. is_square (1 + D * y⇧2)"
using pell_solution_exists by (auto simp: eq_commute[of _ "Suc _"])
hence y: "y > 0 ∧ is_square (1 + D * y⇧2)"
unfolding y_def by (rule LeastI_ex)
have y_le: "y ≤ y'" if "y' > 0" "is_square (1 + D * y'⇧2)" for y'
unfolding y_def using that by (intro Least_le) auto
from y obtain x where x: "x⇧2 = 1 + D * y⇧2"
by (auto elim: is_nth_powerE)
with y have "nontriv_solution (x, y)"
by (auto simp: nontriv_solution_def)
have "is_arg_min (pell_valuation :: nat × nat ⇒ real) nontriv_solution (x, y)"
unfolding is_arg_min_linorder
proof safe
fix a b :: nat
assume *: "nontriv_solution (a, b)"
hence "b > 0" and "Suc (D * b⇧2) = a⇧2"
by (auto simp: nontriv_solution_def intro!: Nat.gr0I)
hence "is_square (1 + D * b⇧2)"
by (auto simp: nontriv_solution_def)
from ‹b > 0› and this have "y ≤ b" by (rule y_le)
with ‹nontriv_solution (x, y)› and * have "x ≤ a"
using solutions_linorder_strict[of x y a b] by (auto simp: nontriv_solution_altdef)
with ‹y ≤ b› show "pell_valuation (int x, int y) ≤ pell_valuation (int a, int b)"
unfolding pell_valuation_def prod.case by (intro add_mono mult_right_mono) auto
qed fact+
thus "∃z. is_arg_min (pell_valuation :: nat × nat ⇒ real) nontriv_solution z" ..
qed
qed
corollary
fund_sol_is_nontriv_solution: "nontriv_solution fund_sol"
and fund_sol_minimal:
"nontriv_solution (a, b) ⟹ pell_valuation fund_sol ≤ pell_valuation (int a, int b)"
and fund_sol_minimal':
"nontriv_solution (z :: nat × nat) ⟹ pell_valuation fund_sol ≤ pell_valuation z"
using fund_sol_is_arg_min by (auto simp: is_arg_min_linorder case_prod_unfold)
lemma fund_sol_minimal'':
assumes "nontriv_solution z"
shows "fst fund_sol ≤ fst z" "snd fund_sol ≤ snd z"
proof -
have "pell_valuation (fst fund_sol, snd fund_sol) ≤ pell_valuation (fst z, snd z)"
using fund_sol_minimal'[OF assms] by (simp add: case_prod_unfold)
hence "fst fund_sol ≤ fst z ∧ snd fund_sol ≤ snd z"
using assms fund_sol_is_nontriv_solution
by (subst solutions_le_iff_pell_valuation_le) (auto simp: case_prod_unfold)
thus "fst fund_sol ≤ fst z" "snd fund_sol ≤ snd z" by blast+
qed
subsection ‹Group structure on solutions›
text ‹
As was mentioned already, the Pell valuation function provides an injective map from
solutions of Pell's equation into the ring $\mathbb{Z}[\sqrt{D}]$. We shall see now that
the solutions are actually a subgroup of the multiplicative group of $\mathbb{Z}[\sqrt{D}]$ via
the valuation function as a homomorphism:
▪ The trivial solution ‹(1, 0)› has valuation ‹1›, which is the neutral element of
$\mathbb{Z}[\sqrt{D}]^*$
▪ Multiplication of two solutions $a + b \sqrt D$ and
$x + y \sqrt D$ leads to $\bar x + \bar y\sqrt D$
with $\bar x = xa + ybD$ and $\bar y = xb + ya$, which is again a solution.
▪ The conjugate ‹(x, -y)› of a solution ‹(x, y)› is an inverse element to this
multiplication operation, since $(x + y \sqrt D) (x - y \sqrt D) = 1$.
›
definition pell_mul :: "('a :: comm_semiring_1 × 'a) ⇒ ('a × 'a) ⇒ ('a × 'a)" where
"pell_mul = (λ(a,b) (x,y). (x * a + y * b * of_nat D, x * b + y * a))"
definition pell_cnj :: "('a :: comm_ring_1 × 'a) ⇒ 'a × 'a" where
"pell_cnj = (λ(a,b). (a, -b))"
lemma pell_cnj_snd_0 [simp]: "snd z = 0 ⟹ pell_cnj z = z"
by (cases z) (simp_all add: pell_cnj_def)
lemma pell_mul_commutes: "pell_mul z1 z2 = pell_mul z2 z1"
by (auto simp: pell_mul_def algebra_simps case_prod_unfold)
lemma pell_mul_assoc: "pell_mul z1 (pell_mul z2 z3) = pell_mul (pell_mul z1 z2) z3"
by (auto simp: pell_mul_def algebra_simps case_prod_unfold)
lemma pell_mul_trivial_left [simp]: "pell_mul (1, 0) z = z"
by (auto simp: pell_mul_def algebra_simps case_prod_unfold)
lemma pell_mul_trivial_right [simp]: "pell_mul z (1, 0) = z"
by (auto simp: pell_mul_def algebra_simps case_prod_unfold)
lemma pell_mul_trivial_left_nat [simp]: "pell_mul (Suc 0, 0) z = z"
by (auto simp: pell_mul_def algebra_simps case_prod_unfold)
lemma pell_mul_trivial_right_nat [simp]: "pell_mul z (Suc 0, 0) = z"
by (auto simp: pell_mul_def algebra_simps case_prod_unfold)
definition pell_power :: "('a :: comm_semiring_1 × 'a) ⇒ nat ⇒ ('a × 'a)" where
"pell_power z n = ((λz'. pell_mul z' z) ^^ n) (1, 0)"
lemma pell_power_0 [simp]: "pell_power z 0 = (1, 0)"
by (simp add: pell_power_def)
lemma pell_power_one [simp]: "pell_power (1, 0) n = (1, 0)"
by (induction n) (auto simp: pell_power_def)
lemma pell_power_one_right [simp]: "pell_power z 1 = z"
by (simp add: pell_power_def)
lemma pell_power_Suc: "pell_power z (Suc n) = pell_mul z (pell_power z n)"
by (simp add: pell_power_def pell_mul_commutes)
lemma pell_power_add: "pell_power z (m + n) = pell_mul (pell_power z m) (pell_power z n)"
by (induction m arbitrary: z )
(simp_all add: funpow_add o_def pell_power_Suc pell_mul_assoc)
lemma pell_valuation_mult [simp]:
"pell_valuation (pell_mul z1 z2) = pell_valuation z1 * pell_valuation z2"
by (simp add: pell_valuation_def pell_mul_def case_prod_unfold algebra_simps)
lemma pell_valuation_mult_nat [simp]:
"pell_valuation (case pell_mul z1 z2 of (a, b) ⇒ (int a, int b)) =
pell_valuation z1 * pell_valuation z2"
by (simp add: pell_valuation_def pell_mul_def case_prod_unfold algebra_simps)
lemma pell_valuation_trivial [simp]: "pell_valuation (1, 0) = 1"
by (simp add: pell_valuation_def)
lemma pell_valuation_trivial_nat [simp]: "pell_valuation (Suc 0, 0) = 1"
by (simp add: pell_valuation_def)
lemma pell_valuation_cnj: "pell_valuation (pell_cnj z) = fst z - snd z * sqrt D"
by (simp add: pell_valuation_def pell_cnj_def case_prod_unfold)
lemma pell_valuation_snd_0 [simp]: "pell_valuation (a, 0) = of_int a"
by (simp add: pell_valuation_def)
lemma pell_valuation_0_iff [simp]: "pell_valuation z = 0 ⟷ z = (0, 0)"
proof
assume *: "pell_valuation z = 0"
have "snd z = 0"
proof (rule ccontr)
assume "snd z ≠ 0"
with * have "sqrt D = -fst z / snd z"
by (simp add: pell_valuation_def case_prod_unfold field_simps)
also have "… ∈ ℚ" by auto
finally show False using nonsquare_D irrat_sqrt_nonsquare by blast
qed
with * have "fst z = 0" by (simp add: pell_valuation_def case_prod_unfold)
with ‹snd z = 0› show "z = (0, 0)" by (cases z) auto
qed (auto simp: pell_valuation_def)
lemma pell_valuation_solution_pos_nat:
fixes z :: "nat × nat"
assumes "solution z"
shows "pell_valuation z > 0"
proof -
from assms have "z ≠ (0, 0)" by (intro notI) auto
hence "pell_valuation z ≠ 0" by (auto split: prod.splits)
moreover have "pell_valuation z ≥ 0" by (intro pell_valuation_nonneg) (auto split: prod.splits)
ultimately show ?thesis by linarith
qed
lemma
assumes "solution z"
shows pell_mul_cnj_right: "pell_mul z (pell_cnj z) = (1, 0)"
and pell_mul_cnj_left: "pell_mul (pell_cnj z) z = (1, 0)"
using assms by (auto simp: pell_mul_def pell_cnj_def solution_def power2_eq_square)
lemma pell_valuation_cnj_solution:
fixes z :: "nat × nat"
assumes "solution z"
shows "pell_valuation (pell_cnj z) = 1 / pell_valuation z"
proof -
have "pell_valuation (pell_cnj z) * pell_valuation z = pell_valuation (pell_mul (pell_cnj z) z)"
by simp
also from assms have "pell_mul (pell_cnj z) z = (1, 0)"
by (subst pell_mul_cnj_left) (auto simp: case_prod_unfold)
finally show ?thesis using pell_valuation_solution_pos_nat[OF assms]
by (auto simp: divide_simps)
qed
lemma pell_valuation_power [simp]: "pell_valuation (pell_power z n) = pell_valuation z ^ n"
by (induction n) (simp_all add: pell_power_Suc)
lemma pell_valuation_power_nat [simp]:
"pell_valuation (case pell_power z n of (a, b) ⇒ (int a, int b)) = pell_valuation z ^ n"
by (induction n) (simp_all add: pell_power_Suc)
lemma pell_valuation_fund_sol_ge_2: "pell_valuation fund_sol ≥ 2"
proof -
obtain x y where [simp]: "fund_sol = (x, y)" by (cases fund_sol)
from fund_sol_is_nontriv_solution have eq: "x⇧2 = 1 + D * y⇧2"
by (auto simp: nontriv_solution_def)
consider "y > 0" | "y = 0" "x ≠ 1"
using fund_sol_is_nontriv_solution by (force simp: nontriv_solution_def)
thus ?thesis
proof cases
assume "y > 0"
hence "1 + 1 * 1 ≤ 1 + D * y⇧2"
using D_pos by (intro add_mono mult_mono) auto
also from eq have "… = x⇧2" ..
finally have "x⇧2 > 1⇧2" by simp
hence "x > 1" by (rule power2_less_imp_less) auto
with ‹y > 0› have "x + y * sqrt D ≥ 2 + 1 * 1"
using D_pos by (intro add_mono mult_mono) auto
thus ?thesis by (simp add: pell_valuation_def)
next
assume [simp]: "y = 0" and "x ≠ 1"
with eq have "x ≠ 0" by (intro notI) auto
with ‹x ≠ 1› have "x ≥ 2" by simp
thus ?thesis by (auto simp: pell_valuation_def)
qed
qed
lemma solution_pell_mul [intro]:
assumes "solution z1" "solution z2"
shows "solution (pell_mul z1 z2)"
proof -
obtain a b where [simp]: "z1 = (a, b)" by (cases z1)
obtain c d where [simp]: "z2 = (c, d)" by (cases z2)
from assms show ?thesis
by (simp add: solution_def pell_mul_def case_prod_unfold power2_eq_square algebra_simps)
qed
lemma solution_pell_cnj [intro]:
assumes "solution z"
shows "solution (pell_cnj z)"
using assms by (auto simp: solution_def pell_cnj_def)
lemma solution_pell_power [simp, intro]: "solution z ⟹ solution (pell_power z n)"
by (induction n) (auto simp: pell_power_Suc)
lemma pell_mul_eq_trivial_nat_iff:
"pell_mul z1 z2 = (Suc 0, 0) ⟷ z1 = (Suc 0, 0) ∧ z2 = (Suc 0, 0)"
using D_gt_1 by (cases z1; cases z2) (auto simp: pell_mul_def)
lemma nontriv_solution_pell_nat_mul1:
"solution (z1 :: nat × nat) ⟹ nontriv_solution z2 ⟹ nontriv_solution (pell_mul z1 z2)"
by (auto simp: nontriv_solution_altdef pell_mul_eq_trivial_nat_iff)
lemma nontriv_solution_pell_nat_mul2:
"nontriv_solution (z1 :: nat × nat) ⟹ solution z2 ⟹ nontriv_solution (pell_mul z1 z2)"
by (auto simp: nontriv_solution_altdef pell_mul_eq_trivial_nat_iff)
lemma nontriv_solution_power_nat [intro]:
assumes "nontriv_solution (z :: nat × nat)" "n > 0"
shows "nontriv_solution (pell_power z n)"
proof -
have "nontriv_solution (pell_power z n) ∨ n = 0"
by (induction n)
(insert assms(1), auto intro: nontriv_solution_pell_nat_mul1 simp: pell_power_Suc)
with assms(2) show ?thesis by auto
qed
subsection ‹The different regions of the valuation function›
text ‹
Next, we shall explore what happens to the valuation function for solutions ‹(x, y)› for
different signs of ‹x› and ‹y›:
▪ If ‹x > 0› and ‹y > 0›, we have $x + y \sqrt D > 1$.
▪ If ‹x > 0› and ‹y < 0›, we have $0 < x + y \sqrt D < 1$.
▪ If ‹x < 0› and ‹y > 0›, we have $-1 < x + y \sqrt D < 0$.
▪ If ‹x < 0› and ‹y < 0›, we have $x + y \sqrt D < -1$.
In particular, this means that we can deduce the sign of ‹x› and ‹y› if we know in
which of these four regions the valuation lies.
›
lemma
assumes "x > 0" "y > 0" "solution (x, y)"
shows pell_valuation_pos_pos: "pell_valuation (x, y) > 1"
and pell_valuation_pos_neg_aux: "pell_valuation (x, -y) ∈ {0<..<1}"
proof -
from D_gt_1 assms have "x + y * sqrt D ≥ 1 + 1 * 1"
by (intro add_mono mult_mono) auto
hence gt_1: "x + y * sqrt D > 1" by simp
thus "pell_valuation (x, y) > 1" by (simp add: pell_valuation_def)
from assms have "1 = x^2 - D * y^2" by (simp add: solution_def)
also have "of_int … = (x - y * sqrt D) * (x + y * sqrt D)"
by (simp add: field_simps power2_eq_square)
finally have eq: "(x - y * sqrt D) = 1 / (x + y * sqrt D)"
using gt_1 by (simp add: field_simps)
note eq
also from gt_1 have "1 / (x + y * sqrt D) < 1 / 1"
by (intro divide_strict_left_mono) auto
finally have "x - y * sqrt D < 1" by simp
note eq
also from gt_1 have "1 / (x + y * sqrt D) > 0"
by (intro divide_pos_pos) auto
finally have "x - y * sqrt D > 0" .
with ‹x - y * sqrt D < 1› show "pell_valuation (x, -y) ∈ {0<..<1}"
by (simp add: pell_valuation_def)
qed
lemma pell_valuation_pos_neg:
assumes "x > 0" "y < 0" "solution (x, y)"
shows "pell_valuation (x, y) ∈ {0<..<1}"
using pell_valuation_pos_neg_aux[of x "-y"] assms by simp
lemma pell_valuation_neg_neg:
assumes "x < 0" "y < 0" "solution (x, y)"
shows "pell_valuation (x, y) < -1"
using pell_valuation_pos_pos[of "-x" "-y"] assms by simp
lemma pell_valuation_neg_pos:
assumes "x < 0" "y > 0" "solution (x, y)"
shows "pell_valuation (x, y) ∈ {-1<..<0}"
using pell_valuation_pos_neg[of "-x" "-y"] assms by simp
lemma pell_valuation_solution_gt1D:
assumes "solution z" "pell_valuation z > 1"
shows "fst z > 0 ∧ snd z > 0"
using pell_valuation_pos_pos[of "fst z" "snd z"] pell_valuation_pos_neg[of "fst z" "snd z"]
pell_valuation_neg_pos[of "fst z" "snd z"] pell_valuation_neg_neg[of "fst z" "snd z"]
assms
by (cases "fst z" "0 :: int" rule: linorder_cases;
cases "snd z" "0 :: int" rule: linorder_cases;
cases z) auto
subsection ‹Generating property of the fundamental solution›
text ‹
We now show that the fundamental solution generates the set of the (non-negative) solutions
in the sense that each solution is a power of the fundamental solution. Combined with the
symmetry property that ‹(x,y)› is a solution iff ‹(±x, ±y)› is a solution, this gives us
a complete characterisation of all solutions of Pell's equation.
›
definition nth_solution :: "nat ⇒ nat × nat" where
"nth_solution n = pell_power fund_sol n"
lemma pell_valuation_nth_solution [simp]:
"pell_valuation (nth_solution n) = pell_valuation fund_sol ^ n"
by (simp add: nth_solution_def)
theorem nth_solution_inj: "inj nth_solution"
proof
fix m n :: nat
assume "nth_solution m = nth_solution n"
hence "pell_valuation (nth_solution m) = pell_valuation (nth_solution n)"
by (simp only: )
also have "pell_valuation (nth_solution m) = pell_valuation fund_sol ^ m"
by simp
also have "pell_valuation (nth_solution n) = pell_valuation fund_sol ^ n"
by simp
finally show "m = n"
using pell_valuation_fund_sol_ge_2 by (subst (asm) power_inject_exp) auto
qed
theorem nth_solution_sound [intro]: "solution (nth_solution n)"
using fund_sol_is_nontriv_solution by (auto simp: nth_solution_def)
theorem nth_solution_sound' [intro]: "n > 0 ⟹ nontriv_solution (nth_solution n)"
using fund_sol_is_nontriv_solution by (auto simp: nth_solution_def)
theorem nth_solution_complete:
fixes z :: "nat × nat"
assumes "solution z"
shows "z ∈ range nth_solution"
proof (cases "z = (1, 0)")
case True
hence "z = nth_solution 0" by (simp add: nth_solution_def)
thus ?thesis by auto
next
case False
with assms have "nontriv_solution z" by (auto simp: nontriv_solution_altdef)
show ?thesis
proof (rule ccontr)
assume "¬?thesis"
hence *: "pell_power fund_sol n ≠ z" for n unfolding nth_solution_def by blast
define u where "u = pell_valuation fund_sol"
define v where "v = pell_valuation z"
define n where "n = nat ⌊log u v⌋"
have u_ge_2: "u ≥ 2" using pell_valuation_fund_sol_ge_2 by (auto simp: u_def)
have v_pos: "v > 0" unfolding v_def using assms
by (intro pell_valuation_solution_pos_nat) auto
have u_le_v: "u ≤ v" unfolding u_def v_def by (rule fund_sol_minimal') fact
have u_power_neq_v: "u ^ k ≠ v" for k
proof
assume "u ^ k = v"
also have "u ^ k = pell_valuation (pell_power fund_sol k)"
by (simp add: u_def)
also have "… = v ⟷ pell_power fund_sol k = z"
unfolding v_def by (subst pell_valuation_eq_iff) (auto split: prod.splits)
finally show False using * by blast
qed
from u_le_v v_pos u_ge_2 have log_ge_1: "log u v ≥ 1"
by (subst one_le_log_cancel_iff) auto
define z' where "z' = pell_mul z (pell_power (pell_cnj fund_sol) n)"
define x and y where "x = nat ¦fst z'¦" and "y = nat ¦snd z'¦"
have "solution z'" using assms fund_sol_is_nontriv_solution unfolding z'_def
by (intro solution_pell_mul solution_pell_power solution_pell_cnj) (auto simp: case_prod_unfold)
have "u ^ n < v"
proof -
from u_ge_2 have "u ^ n = u powr real n" by (subst powr_realpow) auto
also have "… ≤ u powr log u v" using u_ge_2 log_ge_1
by (intro powr_mono) (auto simp: n_def)
also have "… = v"
using u_ge_2 v_pos by (subst powr_log_cancel) auto
finally have "u ^ n ≤ v" .
with u_power_neq_v[of n] show ?thesis by linarith
qed
moreover have "v < u ^ Suc n"
proof -
have "v = u powr log u v"
using u_ge_2 v_pos by (subst powr_log_cancel) auto
also have "log u v ≤ 1 + real_of_int ⌊log u v⌋" by linarith
hence "u powr log u v ≤ u powr real (Suc n)" using u_ge_2 log_ge_1
by (intro powr_mono) (auto simp: n_def)
also have "… = u ^ Suc n" using u_ge_2 by (subst powr_realpow) auto
finally have "u ^ Suc n ≥ v" .
with u_power_neq_v[of "Suc n"] show ?thesis by linarith
qed
ultimately have "v / u ^ n ∈ {1<..<u}"
using u_ge_2 by (simp add: field_simps)
also have "v / u ^ n = pell_valuation z'"
using fund_sol_is_nontriv_solution
by (auto simp add: z'_def u_def v_def pell_valuation_cnj_solution field_simps)
finally have val: "pell_valuation z' ∈ {1<..<u}" .
from val and ‹solution z'› have "nontriv_solution z'"
by (auto simp: nontriv_solution_altdef)
from ‹solution z'› and val have "fst z' > 0 ∧ snd z' > 0"
by (intro pell_valuation_solution_gt1D) auto
hence [simp]: "z' = (int x, int y)"
by (auto simp: x_def y_def)
from ‹nontriv_solution z'› have "pell_valuation (int x, int y) ≥ u"
unfolding u_def by (intro fund_sol_minimal) auto
with val show False by simp
qed
qed
corollary solution_iff_nth_solution:
fixes z :: "nat × nat"
shows "solution z ⟷ z ∈ range nth_solution"
using nth_solution_sound nth_solution_complete by blast
corollary solution_iff_nth_solution':
fixes z :: "int × int"
shows "solution (a, b) ⟷ (nat ¦a¦, nat ¦b¦) ∈ range nth_solution"
proof -
have "solution (a, b) ⟷ solution (nat ¦a¦, nat ¦b¦)"
by simp
also have "… ⟷ (nat ¦a¦, nat ¦b¦) ∈ range nth_solution"
by (rule solution_iff_nth_solution)
finally show ?thesis .
qed
corollary infinite_solutions: "infinite {z :: nat × nat. solution z}"
proof -
have "infinite (range nth_solution)"
by (intro range_inj_infinite nth_solution_inj)
also have "range nth_solution = {z :: nat × nat. solution z}"
by (auto simp: solution_iff_nth_solution)
finally show ?thesis .
qed
corollary infinite_solutions': "infinite {z :: int × int. solution z}"
proof
assume "finite {z :: int × int. solution z}"
hence "finite (map_prod (nat ∘ abs) (nat ∘ abs) ` {z :: int × int. solution z})"
by (rule finite_imageI)
also have "(map_prod (nat ∘ abs) (nat ∘ abs) ` {z :: int × int. solution z}) =
{z :: nat × nat. solution z}"
by (auto simp: map_prod_def image_iff intro!: exI[of _ "int x" for x])
finally show False using infinite_solutions by contradiction
qed
lemma strict_mono_pell_valuation_nth_solution: "strict_mono (pell_valuation ∘ nth_solution)"
using pell_valuation_fund_sol_ge_2
by (auto simp: strict_mono_def intro!: power_strict_increasing)
lemma strict_mono_nth_solution:
"strict_mono (fst ∘ nth_solution)" "strict_mono (snd ∘ nth_solution)"
proof -
let ?g = nth_solution
have "fst (?g m) < fst (?g n) ∧ snd (?g m) < snd (?g n)" if "m < n" for m n
using pell_valuation_fund_sol_ge_2 that
by (subst solutions_less_iff_pell_valuation_less) auto
thus "strict_mono (fst ∘ nth_solution)" "strict_mono (snd ∘ nth_solution)"
by (auto simp: strict_mono_def)
qed
end
subsection ‹The case of an ``almost square'' parameter›
text ‹
If ‹D› is equal to ‹a⇧2 - 1› for some ‹a > 1›, we have a particularly simple case
where the fundamental solution is simply ‹(1, a)›.
›
context
fixes a :: nat
assumes a: "a > 1"
begin
lemma pell_square_minus1: "pell (a⇧2 - Suc 0)"
proof
show "¬is_square (a⇧2 - Suc 0)"
proof
assume "is_square (a⇧2 - Suc 0)"
then obtain k where "k⇧2 = a⇧2 - 1" by (auto elim: is_nth_powerE)
with a have "a⇧2 = Suc (k⇧2)" by simp
hence "a = 1" using pell_square_solution_nat_iff[of "1" a k] by simp
with a show False by simp
qed
qed
interpretation pell "a⇧2 - Suc 0"
by (rule pell_square_minus1)
lemma fund_sol_square_minus1: "fund_sol = (a, 1)"
proof -
from a have sol: "nontriv_solution (a, 1)"
by (simp add: nontriv_solution_def)
from sol have "snd fund_sol ≤ 1"
using fund_sol_minimal''[of "(a, 1)"] by auto
with solutions_linorder_strict[of a 1 "fst fund_sol" "snd fund_sol"]
fund_sol_is_nontriv_solution sol
show "fund_sol = (a, 1)"
by (cases fund_sol) (auto simp: nontriv_solution_altdef)
qed
end
subsection ‹Alternative presentation of the main results›
theorem pell_solutions:
fixes D :: nat
assumes "∄k. D = k⇧2"
obtains x⇩0 y⇩0 :: nat
where "∀(x::int) (y::int).
x⇧2 - D * y⇧2 = 1 ⟷
(∃n::nat. nat ¦x¦ + sqrt D * nat ¦y¦ = (x⇩0 + sqrt D * y⇩0) ^ n)"
proof -
from assms interpret pell
by unfold_locales (auto simp: is_nth_power_def)
show ?thesis
proof (rule that[of "fst fund_sol" "snd fund_sol"], intro allI, goal_cases)
case (1 x y)
have "(x⇧2 - int D * y⇧2 = 1) ⟷ solution (x, y)"
by (auto simp: solution_def)
also have "… ⟷ (∃n. (nat ¦x¦, nat ¦y¦) = nth_solution n)"
by (subst solution_iff_nth_solution') blast
also have "(λn. (nat ¦x¦, nat ¦y¦) = nth_solution n) =
(λn. pell_valuation (nat ¦x¦, nat ¦y¦) = pell_valuation (nth_solution n))"
by (subst pell_valuation_eq_iff) (auto simp add: case_prod_unfold prod_eq_iff fun_eq_iff)
also have "… = (λn. nat ¦x¦ + sqrt D * nat ¦y¦ = (fst fund_sol + sqrt D * snd fund_sol) ^ n)"
by (subst pell_valuation_nth_solution)
(simp add: pell_valuation_def case_prod_unfold mult_ac)
finally show ?case .
qed
qed
corollary pell_solutions_infinite:
fixes D :: nat
assumes "∄k. D = k⇧2"
shows "infinite {(x :: int, y :: int). x⇧2 - D * y⇧2 = 1}"
proof -
from assms interpret pell
by unfold_locales (auto simp: is_nth_power_def)
have "{(x :: int, y :: int). x⇧2 - D * y⇧2 = 1} = {z. solution z}"
by (auto simp: solution_def)
also have "infinite …" by (rule infinite_solutions')
finally show ?thesis .
qed
end