Theory Pell.Pell

(*
  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 = d2"
    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 / u2.
›
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 "xinsert 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  x2 = 1 + D * y2"
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). t2 - u2 * 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 "¦t2 - u2 * 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 "¦t2 - u2 * D¦  1 + 2 * sqrt D" by simp
    hence "t2 - u2 * D  -M" "t2 - u2 * D  M" unfolding M_def by linarith+
    thus "t2 - u2 * 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 {zS. ?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 "{zS. ?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 ` {zS. ?f z = k}  {0..<abs k} × {0..< abs k}"
    using k_nz by (auto simp: case_prod_unfold)
  hence "finite (?h ` {zS. ?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''{zS. ?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''{zS. ?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 "(t12 - u12 * D) = k" "(t22 - u22 * D) = k"
    by (auto simp: S'_def)
  hence "(t12 - u12 * D) * (t22 - u22 * D) = k ^ 2"
    unfolding power2_eq_square by simp
  also have "(t12 - u12 * D) * (t22 - u22 * 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 = k2" .

  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 "k2 dvd k2 + 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 = k2 * (x2 - D * y2)"
    by (simp add: power_mult_distrib algebra_simps)
  finally have eq': "x2 - D * y2 = 1"
    using k_nz by simp
  hence "x2 = 1 + D * y2" by simp
  also have "x2 = int (nat ¦x¦ ^ 2)" by simp
  also have "1 + D * y2 = 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  x2 = 1 + D * y2" 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). a2 = 1 + of_nat D * b2)"

definition nontriv_solution :: "('a × 'a :: comm_semiring_1)  bool" where
  "nontriv_solution = (λ(a, b). (a, b)  (1, 0)  a2 = 1 + of_nat D * b2)"

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 * b2 > 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 "b2  y2"
      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 "a2 = Suc (D * y2)"
      by (simp add: solution_def)
    also from that and assms have " = x2"
      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 * y2))"
    have "y>0. is_square (1 + D * y2)"
      using pell_solution_exists by (auto simp: eq_commute[of _ "Suc _"])
    hence y: "y > 0  is_square (1 + D * y2)"
      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: "x2 = 1 + D * y2"
      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 * b2) = a2"
        by (auto simp: nontriv_solution_def intro!: Nat.gr0I)
      hence "is_square (1 + D * b2)"
        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: "x2 = 1 + D * y2"
    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 * y2"
      using D_pos by (intro add_mono mult_mono) auto
    also from eq have " = x2" ..
    finally have "x2 > 12" 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 a2 - 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 (a2 - Suc 0)"
proof
  show "¬is_square (a2 - Suc 0)"
  proof
    assume "is_square (a2 - Suc 0)"
    then obtain k where "k2 = a2 - 1" by (auto elim: is_nth_powerE)
    with a have "a2 = Suc (k2)" 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 "a2 - 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 = k2"
 obtains x0 y0 :: nat
 where   "(x::int) (y::int).
            x2 - D * y2 = 1 
            (n::nat. nat ¦x¦ + sqrt D * nat ¦y¦ = (x0 + sqrt D * y0) ^ 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 "(x2 - int D * y2 = 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 = k2"
 shows   "infinite {(x :: int, y :: int). x2 - D * y2 = 1}"
proof -
  from assms interpret pell
    by unfold_locales (auto simp: is_nth_power_def)
  have "{(x :: int, y :: int). x2 - D * y2 = 1} = {z. solution z}"
    by (auto simp: solution_def)
  also have "infinite " by (rule infinite_solutions')
  finally show ?thesis .
qed

end