Quelle  Pell.thy

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

  Basic facts about the solutions of Pell's equation
*)
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²"
    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'² = 1 + D * y'²" .
  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²›.
 ›
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² = 1 + D * y²"
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² - u² * 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² - u² * 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² - u² * D∣ ≤ 1 + 2 * sqrt D" by simp
    hence "t² - u² * D ≥ -M" "t² - u² * D ≤ M" unfolding M_def by linarith+
    thus "t² - u² * 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² - u1² * D) = k" "(t2² - u2² * D) = k"
    by (auto simp: S'_def)
  hence "(t1² - u1² * D) * (t2² - u2² * D) = k ^ 2"
    unfolding power2_eq_square by simp
  also have "(t1² - u1² * D) * (t2² - u2² * 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)² - D * (t1 * u2 - t2 * u1)² = k²" .

  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² dvd k² + D * (t1 * u2 - t2 * u1)²"
    using dvd1 by (intro dvd_add) auto
  also from eq have "… = (t1 * t2 - D * u1 * u2)²"
    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)² - D * (k * y)² = k² * (x² - D * y²)"
    by (simp add: power_mult_distrib algebra_simps)
  finally have eq': "x² - D * y² = 1"
    using k_nz by simp
  hence "x² = 1 + D * y²" by simp
  also have "x² = int (nat ∣x∣ ^ 2)" by simp
  also have "1 + D * y² = int (1 + D * nat ∣y∣ ^ 2)" by simp
  also note of_nat_eq_iff
  finally have eq'': "(nat ∣x∣)² = 1 + D * (nat ∣y∣)²" .

  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² = 1 + D * y²" 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² = 1 + of_nat D * b²)"

definition nontriv_solution :: "('a × 'a :: comm_semiring_1) ==> bool" where
  "nontriv_solution = (λ(a, b). (a, b) ≠ (1, 0) ∧ a² = 1 + of_nat D * b²)"

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² > 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² ≤ y²"
      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² = Suc (D * y²)"
      by (simp add: solution_def)
    also from that and assms have "… = x²"
      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²))"
    have "∃y>0. is_square (1 + D * y²)"
      using pell_solution_exists by (auto simp: eq_commute[of _ "Suc _"])
    hence y: "y > 0 ∧ is_square (1 + D * y²)"
      unfolding y_def by (rule LeastI_ex)
    have y_le: "y ≤ y'" if "y' > 0" "is_square (1 + D * y'²)" for y'
      unfolding y_def using that by (intro Least_le) auto
    from y obtain x where x: "x² = 1 + D * y²"
      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²) = a²"
        by (auto simp: nontriv_solution_def intro!: Nat.gr0I)
      hence "is_square (1 + D * b²)"
        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² = 1 + D * y²"
    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²"
      using D_pos by (intro add_mono mult_mono) auto
    also from eq have "… = x²" ..
    finally have "x² > 1²" 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² - 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² - Suc 0)"
proof
  show "¬is_square (a² - Suc 0)"
  proof
    assume "is_square (a² - Suc 0)"
    then obtain k where "k² = a² - 1" by (auto elim: is_nth_powerE)
    with a have "a² = Suc (k²)" 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² - 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²"
 obtains x₀ y₀ :: nat
 where   "∀(x::int) (y::int).
            x² - D * y² = 1 ⟷
            (∃n::nat. nat ∣x∣ + sqrt D * nat ∣y∣ = (x₀ + sqrt D * y₀) ^ 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² - int D * y² = 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²"
 shows   "infinite {(x :: int, y :: int). x² - D * y² = 1}"
proof -
  from assms interpret pell
    by unfold_locales (auto simp: is_nth_power_def)
  have "{(x :: int, y :: int). x² - D * y² = 1} = {z. solution z}"
    by (auto simp: solution_def)
  also have "infinite …" by (rule infinite_solutions')
  finally show ?thesis .
qed

end