Theory Gaussian_Integers

(*
  File:    Gaussian_Integers.thy
  Author:  Manuel Eberl, TU München
*)
section ‹Gaussian Integers›
theory Gaussian_Integers
imports
  "HOL-Computational_Algebra.Computational_Algebra"
  "HOL-Number_Theory.Number_Theory"
begin

subsection ‹Auxiliary material›

lemma coprime_iff_prime_factors_disjoint:
  fixes x y :: "'a :: factorial_semiring"
  assumes "x  0" "y  0"
  shows "coprime x y  prime_factors x  prime_factors y = {}"
proof 
  assume "coprime x y"
  have False if "p  prime_factors x" "p  prime_factors y" for p
  proof -
    from that assms have "p dvd x" "p dvd y"
      by (auto simp: prime_factors_dvd)
    with coprime x y have "p dvd 1"
      using coprime_common_divisor by auto
    with that assms show False by (auto simp: prime_factors_dvd)
  qed
  thus "prime_factors x  prime_factors y = {}" by auto
next
  assume disjoint: "prime_factors x  prime_factors y = {}"
  show "coprime x y"
  proof (rule coprimeI)
    fix d assume d: "d dvd x" "d dvd y"
    show "is_unit d"
    proof (rule ccontr)
      assume "¬is_unit d"
      moreover from this and d assms have "d  0" by auto
      ultimately obtain p where p: "prime p" "p dvd d"
        using prime_divisor_exists by auto
      with d and assms have "p  prime_factors x  prime_factors y"
        by (auto simp: prime_factors_dvd)
      with disjoint show False by auto
    qed
  qed
qed

lemma product_dvd_irreducibleD:
  fixes a b x :: "'a :: algebraic_semidom"
  assumes "irreducible x"
  assumes "a * b dvd x"
  shows "a dvd 1  b dvd 1"
proof -
  from assms obtain c where "x = a * b * c"
    by auto
  hence "x = a * (b * c)"
    by (simp add: mult_ac)
  from irreducibleD[OF assms(1) this] show "a dvd 1  b dvd 1"
    by (auto simp: is_unit_mult_iff)
qed

lemma prime_elem_mult_dvdI:
  assumes "prime_elem p" "p dvd c" "b dvd c" "¬p dvd b"
  shows   "p * b dvd c"
proof -
  from assms(3) obtain a where c: "c = a * b"
    using mult.commute by blast
  with assms(2) have "p dvd a * b"
    by simp
  with assms have "p dvd a"
    by (subst (asm) prime_elem_dvd_mult_iff) auto
  with c show ?thesis by (auto intro: mult_dvd_mono)
qed

lemma prime_elem_power_mult_dvdI:
  fixes p :: "'a :: algebraic_semidom"
  assumes "prime_elem p" "p ^ n dvd c" "b dvd c" "¬p dvd b"
  shows   "p ^ n * b dvd c"
proof (cases "n = 0")
  case False
  from assms(3) obtain a where c: "c = a * b"
    using mult.commute by blast
  with assms(2) have "p ^ n dvd b * a"
    by (simp add: mult_ac)
  hence "p ^ n dvd a"
    by (rule prime_power_dvd_multD[OF assms(1)]) (use assms False in auto)
  with c show ?thesis by (auto intro: mult_dvd_mono)
qed (use assms in auto)

lemma prime_mod_4_cases:
  fixes p :: nat
  assumes "prime p"
  shows   "p = 2  [p = 1] (mod 4)  [p = 3] (mod 4)"
proof (cases "p = 2")
  case False
  with prime_gt_1_nat[of p] assms have "p > 2" by auto
  have "¬4 dvd p"
    using assms product_dvd_irreducibleD[of p 2 2]
    by (auto simp: prime_elem_iff_irreducible simp flip: prime_elem_nat_iff)
  hence "p mod 4  0"
    by (auto simp: mod_eq_0_iff_dvd)
  moreover have "p mod 4  2"
  proof
    assume "p mod 4 = 2"
    hence "p mod 4 mod 2 = 0"
      by (simp add: cong_def)
    thus False using prime p p > 2 prime_odd_nat[of p]
      by (auto simp: mod_mod_cancel)
  qed
  moreover have "p mod 4  {0,1,2,3}"
    by auto
  ultimately show ?thesis by (auto simp: cong_def)
qed auto

lemma of_nat_prod_mset: "of_nat (prod_mset A) = prod_mset (image_mset of_nat A)"
  by (induction A) auto

lemma multiplicity_0_left [simp]: "multiplicity 0 x = 0"
  by (cases "x = 0") (auto simp: not_dvd_imp_multiplicity_0)

lemma is_unit_power [intro]: "is_unit x  is_unit (x ^ n)"
  by (subst is_unit_power_iff) auto

lemma (in factorial_semiring) pow_divides_pow_iff:
  assumes "n > 0"
  shows   "a ^ n dvd b ^ n  a dvd b"
proof (cases "b = 0")
  case False
  show ?thesis
  proof
    assume dvd: "a ^ n dvd b ^ n"
    with b  0 have "a  0"
      using n > 0 by (auto simp: power_0_left)
    show "a dvd b"
    proof (rule multiplicity_le_imp_dvd)
      fix p :: 'a assume p: "prime p"
      from dvd b  0 have "multiplicity p (a ^ n)  multiplicity p (b ^ n)"
        by (intro dvd_imp_multiplicity_le) auto
      thus "multiplicity p a  multiplicity p b"
        using p a  0 b  0 n > 0 by (simp add: prime_elem_multiplicity_power_distrib)
    qed fact+
  qed (auto intro: dvd_power_same)
qed (use assms in auto simp: power_0_left)

lemma multiplicity_power_power:
  fixes p :: "'a :: {factorial_semiring, algebraic_semidom}"
  assumes "n > 0"
  shows   "multiplicity (p ^ n) (x ^ n) = multiplicity p x"
proof (cases "x = 0  p = 0  is_unit p")
  case True
  thus ?thesis using n > 0
    by (auto simp: power_0_left is_unit_power_iff multiplicity_unit_left)
next
  case False
  show ?thesis
  proof (intro antisym multiplicity_geI)
    have "(p ^ multiplicity p x) ^ n dvd x ^ n"
      by (intro dvd_power_same) (simp add: multiplicity_dvd)
    thus "(p ^ n) ^ multiplicity p x dvd x ^ n"
      by (simp add: mult_ac flip: power_mult)
  next
    have "(p ^ n) ^ multiplicity (p ^ n) (x ^ n) dvd x ^ n"
      by (simp add: multiplicity_dvd)
    hence "(p ^ multiplicity (p ^ n) (x ^ n)) ^ n dvd x ^ n"
      by (simp add: mult_ac flip: power_mult)
    thus "p ^ multiplicity (p ^ n) (x ^ n) dvd x"
      by (subst (asm) pow_divides_pow_iff) (use assms in auto)
  qed (use False n > 0 in auto simp: is_unit_power_iff)
qed

lemma even_square_cong_4_int:
  [x2 = 0] (mod 4) if even x for x :: int
proof -
  from that obtain y where x = 2 * y ..
  then show ?thesis by (simp add: cong_def)
qed

lemma even_square_cong_4_nat:
  [x2 = 0] (mod 4) if even x for x :: nat
  using that even_square_cong_4_int [of int x] by (simp flip: cong_int_iff)

lemma odd_square_cong_4_int:
  [x2 = 1] (mod 4) if odd x for x :: int
proof -
  from that obtain y where x = 2 * y + 1 ..
  then have x2 = 4 * (y2 + y) + 1
    by (simp add: power2_eq_square algebra_simps)
  also have  mod 4 = ((4 * (y2 + y)) mod 4 + 1 mod 4) mod 4
    by (simp only: mod_simps)
  also have  = 1 mod 4
    by simp
  finally show ?thesis
    by (simp only: cong_def)
qed

lemma odd_square_cong_4_nat:
  [x2 = 1] (mod 4) if odd x for x :: nat
  using that odd_square_cong_4_int [of int x] by (simp flip: cong_int_iff)


text ‹
  Gaussian integers will require a notion of an element being a power up to a unit,
  so we introduce this here. This should go in the library eventually.
›
definition is_nth_power_upto_unit where
  "is_nth_power_upto_unit n x  (u. is_unit u  is_nth_power n (u * x))"

lemma is_nth_power_upto_unit_base: "is_nth_power n x  is_nth_power_upto_unit n x"
  by (auto simp: is_nth_power_upto_unit_def intro: exI[of _ 1])

lemma is_nth_power_upto_unitI:
  assumes "normalize (x ^ n) = normalize y"
  shows   "is_nth_power_upto_unit n y"
proof -
  from associatedE1[OF assms] obtain u where "is_unit u" "u * y = x ^ n"
    by metis
  thus ?thesis
    by (auto simp: is_nth_power_upto_unit_def intro!: exI[of _ u])
qed

lemma is_nth_power_upto_unit_conv_multiplicity: 
  fixes x :: "'a :: factorial_semiring"
  assumes "n > 0"
  shows   "is_nth_power_upto_unit n x  (p. prime p  n dvd multiplicity p x)"
proof (cases "x = 0")
  case False
  show ?thesis
  proof safe
    fix p :: 'a assume p: "prime p"
    assume "is_nth_power_upto_unit n x"
    then obtain u y where uy: "is_unit u" "u * x = y ^ n"
      by (auto simp: is_nth_power_upto_unit_def elim!: is_nth_powerE)
    from p uy assms False have [simp]: "y  0" by (auto simp: power_0_left)
    have "multiplicity p (u * x) = multiplicity p (y ^ n)"
      by (subst uy(2) [symmetric]) simp
    also have "multiplicity p (u * x) = multiplicity p x"
      by (simp add: multiplicity_times_unit_right uy(1))
    finally show "n dvd multiplicity p x"
      using False and p and uy and assms
      by (auto simp: prime_elem_multiplicity_power_distrib)
  next
    assume *: "p. prime p  n dvd multiplicity p x"
    have "multiplicity p ((pprime_factors x. p ^ (multiplicity p x div n)) ^ n) = 
            multiplicity p x" if "prime p" for p
    proof -
      from that and * have "n dvd multiplicity p x" by blast
      have "multiplicity p x = 0" if "p  prime_factors x"
        using that and prime p by (simp add: prime_factors_multiplicity)
      with that and * and assms show ?thesis unfolding prod_power_distrib power_mult [symmetric]
        by (subst multiplicity_prod_prime_powers) (auto simp: in_prime_factors_imp_prime elim: dvdE)
    qed
    with assms False 
      have "normalize ((pprime_factors x. p ^ (multiplicity p x div n)) ^ n) = normalize x"
      by (intro multiplicity_eq_imp_eq) (auto simp: multiplicity_prod_prime_powers)
    thus "is_nth_power_upto_unit n x"
      by (auto intro: is_nth_power_upto_unitI)
  qed
qed (use assms in auto simp: is_nth_power_upto_unit_def)

lemma is_nth_power_upto_unit_0_left [simp, intro]: "is_nth_power_upto_unit 0 x  is_unit x"
proof
  assume "is_unit x"
  thus "is_nth_power_upto_unit 0 x"
    unfolding is_nth_power_upto_unit_def by (intro exI[of _ "1 div x"]) auto
next
  assume "is_nth_power_upto_unit 0 x"
  then obtain u where "is_unit u" "u * x = 1"
    by (auto simp: is_nth_power_upto_unit_def)
  thus "is_unit x"
    by (metis dvd_triv_right)
qed

lemma is_nth_power_upto_unit_unit [simp, intro]:
  assumes "is_unit x"
  shows   "is_nth_power_upto_unit n x"
  using assms by (auto simp: is_nth_power_upto_unit_def intro!: exI[of _ "1 div x"])

lemma is_nth_power_upto_unit_1_left [simp, intro]: "is_nth_power_upto_unit 1 x"
  by (auto simp: is_nth_power_upto_unit_def intro: exI[of _ 1])

lemma is_nth_power_upto_unit_mult_coprimeD1:
  fixes x y :: "'a :: factorial_semiring"
  assumes "coprime x y" "is_nth_power_upto_unit n (x * y)"
  shows   "is_nth_power_upto_unit n x"
proof -
  consider "n = 0" | "x = 0" "n > 0" | "x  0" "y = 0" "n > 0" | "n > 0" "x  0" "y  0"
    by force
  thus ?thesis
  proof cases
    assume [simp]: "n = 0"
    from assms have "is_unit (x * y)"
      by auto
    hence "is_unit x"
      using is_unit_mult_iff by blast
    thus ?thesis using assms by auto
  next
    assume "x = 0" "n > 0"
    thus ?thesis by (auto simp: is_nth_power_upto_unit_def)
  next
    assume *: "x  0" "y = 0" "n > 0"
    with assms show ?thesis by auto
  next
    assume *: "n > 0" and [simp]: "x  0" "y  0"
    show ?thesis
    proof (subst is_nth_power_upto_unit_conv_multiplicity[OF n > 0]; safe)
      fix p :: 'a assume p: "prime p"
      show "n dvd multiplicity p x"
      proof (cases "p dvd x")
        case False
        thus ?thesis
          by (simp add: not_dvd_imp_multiplicity_0)
      next
        case True
        have "n dvd multiplicity p (x * y)"
          using assms(2) n > 0 p by (auto simp: is_nth_power_upto_unit_conv_multiplicity)
        also have " = multiplicity p x + multiplicity p y"
          using p by (subst prime_elem_multiplicity_mult_distrib) auto
        also have "¬p dvd y"
          using coprime x y p dvd x p not_prime_unit coprime_common_divisor by blast
        hence "multiplicity p y = 0"
          by (rule not_dvd_imp_multiplicity_0)
        finally show ?thesis by simp
      qed
    qed
  qed
qed

lemma is_nth_power_upto_unit_mult_coprimeD2:
  fixes x y :: "'a :: factorial_semiring"
  assumes "coprime x y" "is_nth_power_upto_unit n (x * y)"
  shows   "is_nth_power_upto_unit n y"
  using assms is_nth_power_upto_unit_mult_coprimeD1[of y x]
  by (simp_all add: mult_ac coprime_commute)


subsection ‹Definition›

text ‹
  Gaussian integers are the ring $\mathbb{Z}[i]$ which is formed either by formally adjoining
  an element $i$ with $i^2 = -1$ to $\mathbb{Z}$ or by taking all the complex numbers with
  integer real and imaginary part.

  We define them simply by giving an appropriate ring structure to $\mathbb{Z}^2$, with the first
  component representing the real part and the second component the imaginary part:
›
codatatype gauss_int = Gauss_Int (ReZ: int) (ImZ: int)

text ‹
  The following is the imaginary unit $i$ in the Gaussian integers, which we will denote as
  𝗂:
›
primcorec gauss_i where
  "ReZ gauss_i = 0"
| "ImZ gauss_i = 1"

lemma gauss_int_eq_iff: "x = y  ReZ x = ReZ y  ImZ x = ImZ y"
  by (cases x; cases y) auto


(*<*)
bundle gauss_int_notation
begin

notation gauss_i ("𝗂")

end

bundle no_gauss_int_notation
begin

no_notation (output) gauss_i ("𝗂")

end

bundle gauss_int_output_notation
begin

notation (output) gauss_i ("𝔦")

end

unbundle gauss_int_notation
(*>*)


text ‹
  Next, we define the canonical injective homomorphism from the Gaussian integers into the
  complex numbers:
›
primcorec gauss2complex where
  "Re (gauss2complex z) = of_int (ReZ z)"
| "Im (gauss2complex z) = of_int (ImZ z)"

declare [[coercion gauss2complex]]

lemma gauss2complex_eq_iff [simp]: "gauss2complex z = gauss2complex u  z = u"
  by (simp add: complex_eq_iff gauss_int_eq_iff)

text ‹
  Gaussian integers also have conjugates, just like complex numbers:
›
primcorec gauss_cnj where
  "ReZ (gauss_cnj z) = ReZ z"
| "ImZ (gauss_cnj z) = -ImZ z"


text ‹
  In the remainder of this section, we prove that Gaussian integers are a commutative ring
  of characteristic 0 and several other trivial algebraic properties.
›

instantiation gauss_int :: comm_ring_1
begin

primcorec zero_gauss_int where
  "ReZ zero_gauss_int = 0"
| "ImZ zero_gauss_int = 0"

primcorec one_gauss_int where
  "ReZ one_gauss_int = 1"
| "ImZ one_gauss_int = 0"

primcorec uminus_gauss_int where
  "ReZ (uminus_gauss_int x) = -ReZ x"
| "ImZ (uminus_gauss_int x) = -ImZ x"

primcorec plus_gauss_int where
  "ReZ (plus_gauss_int x y) = ReZ x + ReZ y"
| "ImZ (plus_gauss_int x y) = ImZ x + ImZ y"

primcorec minus_gauss_int where
  "ReZ (minus_gauss_int x y) = ReZ x - ReZ y"
| "ImZ (minus_gauss_int x y) = ImZ x - ImZ y"

primcorec times_gauss_int where
  "ReZ (times_gauss_int x y) = ReZ x * ReZ y - ImZ x * ImZ y"
| "ImZ (times_gauss_int x y) = ReZ x * ImZ y + ImZ x * ReZ y"

instance
  by intro_classes (auto simp: gauss_int_eq_iff algebra_simps)

end

lemma gauss_i_times_i [simp]: "𝗂 * 𝗂 = (-1 :: gauss_int)"
  and gauss_cnj_i [simp]: "gauss_cnj 𝗂 = -𝗂"
  by (simp_all add: gauss_int_eq_iff)

lemma gauss_cnj_eq_0_iff [simp]: "gauss_cnj z = 0  z = 0"
  by (auto simp: gauss_int_eq_iff)

lemma gauss_cnj_eq_self: "Im z = 0  gauss_cnj z = z"
  and gauss_cnj_eq_minus_self: "Re z = 0  gauss_cnj z = -z"
  by (auto simp: gauss_int_eq_iff)

lemma ReZ_of_nat [simp]: "ReZ (of_nat n) = of_nat n"
  and ImZ_of_nat [simp]: "ImZ (of_nat n) = 0"
  by (induction n; simp)+

lemma ReZ_of_int [simp]: "ReZ (of_int n) = n"
  and ImZ_of_int [simp]: "ImZ (of_int n) = 0"
  by (induction n; simp)+

lemma ReZ_numeral [simp]: "ReZ (numeral n) = numeral n"
  and ImZ_numeral [simp]: "ImZ (numeral n) = 0"
  by (subst of_nat_numeral [symmetric], subst ReZ_of_nat ImZ_of_nat, simp)+

lemma gauss2complex_0 [simp]: "gauss2complex 0 = 0"
  and gauss2complex_1 [simp]: "gauss2complex 1 = 1"
  and gauss2complex_i [simp]: "gauss2complex 𝗂 = 𝗂"
  and gauss2complex_add [simp]: "gauss2complex (x + y) = gauss2complex x + gauss2complex y"
  and gauss2complex_diff [simp]: "gauss2complex (x - y) = gauss2complex x - gauss2complex y"
  and gauss2complex_mult [simp]: "gauss2complex (x * y) = gauss2complex x * gauss2complex y"
  and gauss2complex_uminus [simp]: "gauss2complex (-x) = -gauss2complex x"
  and gauss2complex_cnj [simp]: "gauss2complex (gauss_cnj x) = cnj (gauss2complex x)"
  by (simp_all add: complex_eq_iff)

lemma gauss2complex_of_nat [simp]: "gauss2complex (of_nat n) = of_nat n"
  by (simp add: complex_eq_iff)

lemma gauss2complex_eq_0_iff [simp]: "gauss2complex x = 0  x = 0"
  and gauss2complex_eq_1_iff [simp]: "gauss2complex x = 1  x = 1"
  and zero_eq_gauss2complex_iff [simp]: "0 = gauss2complex x  x = 0"
  and one_eq_gauss2complex_iff [simp]: "1 = gauss2complex x  x = 1"
  by (simp_all add: complex_eq_iff gauss_int_eq_iff)

lemma gauss_i_times_gauss_i_times [simp]: "𝗂 * (𝗂 * x) = (-x :: gauss_int)"
  by (subst mult.assoc [symmetric], subst gauss_i_times_i) auto

lemma gauss_i_neq_0 [simp]: "𝗂  0" "0  𝗂"
  and gauss_i_neq_1 [simp]: "𝗂  1" "1  𝗂"
  and gauss_i_neq_of_nat [simp]: "𝗂  of_nat n" "of_nat n  𝗂"
  and gauss_i_neq_of_int [simp]: "𝗂  of_int n" "of_int n  𝗂"
  and gauss_i_neq_numeral [simp]: "𝗂  numeral m" "numeral m  𝗂"
  by (auto simp: gauss_int_eq_iff)

lemma gauss_cnj_0 [simp]: "gauss_cnj 0 = 0"
  and gauss_cnj_1 [simp]: "gauss_cnj 1 = 1"
  and gauss_cnj_cnj [simp]: "gauss_cnj (gauss_cnj z) = z"
  and gauss_cnj_uminus [simp]: "gauss_cnj (-a) = -gauss_cnj a"
  and gauss_cnj_add [simp]: "gauss_cnj (a + b) = gauss_cnj a + gauss_cnj b"
  and gauss_cnj_diff [simp]: "gauss_cnj (a - b) = gauss_cnj a - gauss_cnj b"
  and gauss_cnj_mult [simp]: "gauss_cnj (a * b) = gauss_cnj a * gauss_cnj b"
  and gauss_cnj_of_nat [simp]: "gauss_cnj (of_nat n1) = of_nat n1"
  and gauss_cnj_of_int [simp]: "gauss_cnj (of_int n2) = of_int n2"
  and gauss_cnj_numeral [simp]: "gauss_cnj (numeral n3) = numeral n3"
  by (simp_all add: gauss_int_eq_iff)

lemma gauss_cnj_power [simp]: "gauss_cnj (a ^ n) = gauss_cnj a ^ n"
  by (induction n) auto

lemma gauss_cnj_sum [simp]: "gauss_cnj (sum f A) = (xA. gauss_cnj (f x))"
  by (induction A rule: infinite_finite_induct) auto

lemma gauss_cnj_prod [simp]: "gauss_cnj (prod f A) = (xA. gauss_cnj (f x))"
  by (induction A rule: infinite_finite_induct) auto

lemma of_nat_dvd_of_nat:
  assumes "a dvd b"
  shows   "of_nat a dvd (of_nat b :: 'a :: comm_semiring_1)"
  using assms by auto

lemma of_int_dvd_imp_dvd_gauss_cnj:
  fixes z :: gauss_int
  assumes "of_int n dvd z"
  shows   "of_int n dvd gauss_cnj z"
proof -
  from assms obtain u where "z = of_int n * u" by blast
  hence "gauss_cnj z = of_int n * gauss_cnj u"
    by simp
  thus ?thesis by auto
qed

lemma of_nat_dvd_imp_dvd_gauss_cnj:
  fixes z :: gauss_int
  assumes "of_nat n dvd z"
  shows   "of_nat n dvd gauss_cnj z"
  using of_int_dvd_imp_dvd_gauss_cnj[of "int n"] assms by simp

lemma of_int_dvd_of_int_gauss_int_iff:
  "(of_int m :: gauss_int) dvd of_int n  m dvd n"
proof
  assume "of_int m dvd (of_int n :: gauss_int)"
  then obtain a :: gauss_int where "of_int n = of_int m * a"
    by blast
  thus "m dvd n"
    by (auto simp: gauss_int_eq_iff)
qed auto

lemma of_nat_dvd_of_nat_gauss_int_iff:
  "(of_nat m :: gauss_int) dvd of_nat n  m dvd n"
  using of_int_dvd_of_int_gauss_int_iff[of "int m" "int n"] by simp

lemma gauss_cnj_dvd:
  assumes "a dvd b"
  shows   "gauss_cnj a dvd gauss_cnj b"
proof -
  from assms obtain c where "b = a * c"
    by blast
  hence "gauss_cnj b = gauss_cnj a * gauss_cnj c"
    by simp
  thus ?thesis by auto
qed

lemma gauss_cnj_dvd_iff: "gauss_cnj a dvd gauss_cnj b  a dvd b"
  using gauss_cnj_dvd[of a b] gauss_cnj_dvd[of "gauss_cnj a" "gauss_cnj b"] by auto

lemma gauss_cnj_dvd_left_iff: "gauss_cnj a dvd b  a dvd gauss_cnj b"
  by (subst gauss_cnj_dvd_iff [symmetric]) auto

lemma gauss_cnj_dvd_right_iff: "a dvd gauss_cnj b  gauss_cnj a dvd b"
  by (rule gauss_cnj_dvd_left_iff [symmetric])


instance gauss_int :: idom
proof
  fix z u :: gauss_int
  assume "z  0" "u  0"
  hence "gauss2complex z * gauss2complex u  0"
    by simp
  also have "gauss2complex z * gauss2complex u = gauss2complex (z * u)"
    by simp
  finally show "z * u  0"
    unfolding gauss2complex_eq_0_iff .
qed

instance gauss_int :: ring_char_0
  by intro_classes (auto intro!: injI simp: gauss_int_eq_iff)


subsection ‹Pretty-printing›

text ‹
  The following lemma collection provides better pretty-printing of Gaussian integers so that
  e.g.\ evaluation with the `value' command produces nicer results.
›
lemma gauss_int_code_post [code_post]:
  "Gauss_Int 0 0 = 0"
  "Gauss_Int 0 1 = 𝗂"
  "Gauss_Int 0 (-1) = -𝗂"
  "Gauss_Int 1 0 = 1"
  "Gauss_Int 1 1 = 1 + 𝗂"
  "Gauss_Int 1 (-1) = 1 - 𝗂"
  "Gauss_Int (-1) 0 = -1"
  "Gauss_Int (-1) 1 = -1 + 𝗂"
  "Gauss_Int (-1) (-1) = -1 - 𝗂"  
  "Gauss_Int (numeral b) 0 = numeral b"
  "Gauss_Int (-numeral b) 0 = -numeral b"
  "Gauss_Int (numeral b) 1 = numeral b + 𝗂"
  "Gauss_Int (-numeral b) 1 = -numeral b + 𝗂"
  "Gauss_Int (numeral b) (-1) = numeral b - 𝗂"
  "Gauss_Int (-numeral b) (-1) = -numeral b - 𝗂"
  "Gauss_Int 0 (numeral b) = numeral b * 𝗂"
  "Gauss_Int 0 (-numeral b) = -numeral b * 𝗂"
  "Gauss_Int 1 (numeral b) = 1 + numeral b * 𝗂"
  "Gauss_Int 1 (-numeral b) = 1 - numeral b * 𝗂"
  "Gauss_Int (-1) (numeral b) = -1 + numeral b * 𝗂"
  "Gauss_Int (-1) (-numeral b) = -1 - numeral b * 𝗂"
  "Gauss_Int (numeral a) (numeral b) = numeral a + numeral b * 𝗂"
  "Gauss_Int (numeral a) (-numeral b) = numeral a - numeral b * 𝗂"
  "Gauss_Int (-numeral a) (numeral b) = -numeral a + numeral b * 𝗂"
  "Gauss_Int (-numeral a) (-numeral b) = -numeral a - numeral b * 𝗂"
  by (simp_all add: gauss_int_eq_iff)

value "𝗂 ^ 3"
value "2 * (3 + 𝗂)"
value "(2 + 𝗂) * (2 - 𝗂)"


subsection ‹Norm›

text ‹
  The square of the complex norm (or complex modulus) on the Gaussian integers gives us a norm
  that always returns a natural number. We will later show that this is also a Euclidean norm
  (in the sense of a Euclidean ring).
›
definition gauss_int_norm :: "gauss_int  nat" where
  "gauss_int_norm z = nat (ReZ z ^ 2 + ImZ z ^ 2)"

lemma gauss_int_norm_0 [simp]: "gauss_int_norm 0 = 0"
  and gauss_int_norm_1 [simp]: "gauss_int_norm 1 = 1"
  and gauss_int_norm_i [simp]: "gauss_int_norm 𝗂 = 1"
  and gauss_int_norm_cnj [simp]: "gauss_int_norm (gauss_cnj z) = gauss_int_norm z"
  and gauss_int_norm_of_nat [simp]: "gauss_int_norm (of_nat n) = n ^ 2"
  and gauss_int_norm_of_int [simp]: "gauss_int_norm (of_int m) = nat (m ^ 2)"
  and gauss_int_norm_of_numeral [simp]: "gauss_int_norm (numeral n') = numeral (Num.sqr n')"
  by (simp_all add: gauss_int_norm_def nat_power_eq)

lemma gauss_int_norm_uminus [simp]: "gauss_int_norm (-z) = gauss_int_norm z"
  by (simp add: gauss_int_norm_def)

lemma gauss_int_norm_eq_0_iff [simp]: "gauss_int_norm z = 0  z = 0"
proof
  assume "gauss_int_norm z = 0"
  hence "ReZ z ^ 2 + ImZ z ^ 2  0"
    by (simp add: gauss_int_norm_def)
  moreover have "ReZ z ^ 2 + ImZ z ^ 2  0"
    by simp
  ultimately have "ReZ z ^ 2 + ImZ z ^ 2 = 0"
    by linarith
  thus "z = 0"
    by (auto simp: gauss_int_eq_iff)
qed auto

lemma gauss_int_norm_pos_iff [simp]: "gauss_int_norm z > 0  z  0"
  using gauss_int_norm_eq_0_iff[of z] by (auto intro: Nat.gr0I)

lemma real_gauss_int_norm: "real (gauss_int_norm z) = norm (gauss2complex z) ^ 2"
  by (auto simp: cmod_def gauss_int_norm_def)

lemma gauss_int_norm_mult: "gauss_int_norm (z * u) = gauss_int_norm z * gauss_int_norm u"
proof -
  have "real (gauss_int_norm (z * u)) = real (gauss_int_norm z * gauss_int_norm u)"
    unfolding of_nat_mult by (simp add: real_gauss_int_norm norm_power norm_mult power_mult_distrib)
  thus ?thesis by (subst (asm) of_nat_eq_iff)
qed

lemma self_mult_gauss_cnj: "z * gauss_cnj z = of_nat (gauss_int_norm z)"
  by (simp add: gauss_int_norm_def gauss_int_eq_iff algebra_simps power2_eq_square)

lemma gauss_cnj_mult_self: "gauss_cnj z * z = of_nat (gauss_int_norm z)"
  by (subst mult.commute, rule self_mult_gauss_cnj)

lemma self_plus_gauss_cnj: "z + gauss_cnj z = of_int (2 * ReZ z)"
  and self_minus_gauss_cnj: "z - gauss_cnj z = of_int (2 * ImZ z) * 𝗂"
  by (auto simp: gauss_int_eq_iff)

lemma gauss_int_norm_dvd_mono:
  assumes "a dvd b"
  shows "gauss_int_norm a dvd gauss_int_norm b"
proof -
  from assms obtain c where "b = a * c" by blast
  hence "gauss_int_norm b = gauss_int_norm (a * c)"
    by metis
  thus ?thesis by (simp add: gauss_int_norm_mult)
qed

text ‹
  A Gaussian integer is a unit iff its norm is 1, and this is the case precisely for the four
  elements ±1› and ±𝗂›:
›
lemma is_unit_gauss_int_iff: "x dvd 1  x  {1, -1, 𝗂, -𝗂 :: gauss_int}"
  and is_unit_gauss_int_iff': "x dvd 1  gauss_int_norm x = 1"
proof -
  have "x dvd 1" if "x  {1, -1, 𝗂, -𝗂}"
  proof -
    from that have *: "x * gauss_cnj x = 1"
      by (auto simp: gauss_int_norm_def)
    show "x dvd 1" by (subst * [symmetric]) simp
  qed
  moreover have "gauss_int_norm x = 1" if "x dvd 1"
    using gauss_int_norm_dvd_mono[OF that] by simp
  moreover have "x  {1, -1, 𝗂, -𝗂}" if "gauss_int_norm x = 1"
  proof -
    from that have *: "(ReZ x)2 + (ImZ x)2 = 1"
      by (auto simp: gauss_int_norm_def nat_eq_iff)
    hence "ReZ x ^ 2  1" and "ImZ x ^ 2  1"
      using zero_le_power2[of "ImZ x"] zero_le_power2[of "ReZ x"] by linarith+
    hence "¦ReZ x¦  1" and "¦ImZ x¦  1"
      by (auto simp: abs_square_le_1)
    hence "ReZ x  {-1, 0, 1}" and "ImZ x  {-1, 0, 1}"
      by auto
    thus "x  {1, -1, 𝗂, -𝗂 :: gauss_int}"
      using * by (auto simp: gauss_int_eq_iff)    
  qed
  ultimately show "x dvd 1  x  {1, -1, 𝗂, -𝗂 :: gauss_int}"
              and "x dvd 1  gauss_int_norm x = 1"
    by blast+
qed

lemma is_unit_gauss_i [simp, intro]: "(gauss_i :: gauss_int) dvd 1"
  by (simp add: is_unit_gauss_int_iff)

lemma gauss_int_norm_eq_Suc_0_iff: "gauss_int_norm x = Suc 0  x dvd 1"
  by (simp add: is_unit_gauss_int_iff')

lemma is_unit_gauss_cnj [intro]: "z dvd 1  gauss_cnj z dvd 1"
  by (simp add: is_unit_gauss_int_iff')

lemma is_unit_gauss_cnj_iff [simp]: "gauss_cnj z dvd 1  z dvd 1"
  by (simp add: is_unit_gauss_int_iff')


subsection ‹Division and normalisation›

text ‹
  We define a rounding operation that takes a complex number and returns a Gaussian integer
  by rounding the real and imaginary parts separately:
›
primcorec round_complex :: "complex  gauss_int" where
  "ReZ (round_complex z) = round (Re z)"
| "ImZ (round_complex z) = round (Im z)"

text ‹
  The distance between a rounded complex number and the original one is no more than
  $\frac{1}{2}\sqrt{2}$:
›
lemma norm_round_complex_le: "norm (z - gauss2complex (round_complex z)) ^ 2  1 / 2"
proof -
  have "(Re z - ReZ (round_complex z)) ^ 2  (1 / 2) ^ 2"
    using of_int_round_abs_le[of "Re z"]
    by (subst abs_le_square_iff [symmetric]) (auto simp: abs_minus_commute)
  moreover have "(Im z - ImZ (round_complex z)) ^ 2  (1 / 2) ^ 2"
    using of_int_round_abs_le[of "Im z"]
    by (subst abs_le_square_iff [symmetric]) (auto simp: abs_minus_commute)
  ultimately have "(Re z - ReZ (round_complex z)) ^ 2 + (Im z - ImZ (round_complex z)) ^ 2 
                     (1 / 2) ^ 2 + (1 / 2) ^ 2"
    by (rule add_mono)
  thus "norm (z - gauss2complex (round_complex z)) ^ 2  1 / 2"
    by (simp add: cmod_def power2_eq_square)
qed

lemma dist_round_complex_le: "dist z (gauss2complex (round_complex z))  sqrt 2 / 2"
proof -
  have "dist z (gauss2complex (round_complex z)) ^ 2 =
        norm (z - gauss2complex (round_complex z)) ^ 2"
    by (simp add: dist_norm)
  also have "  1 / 2"
    by (rule norm_round_complex_le)
  also have " = (sqrt 2 / 2) ^ 2"
    by (simp add: power2_eq_square)
  finally show ?thesis
    by (rule power2_le_imp_le) auto
qed


text ‹
  We can now define division on Gaussian integers simply by performing the division in the
  complex numbers and rounding the result. This also gives us a remainder operation defined
  accordingly for which the norm of the remainder is always smaller than the norm of the divisor.

  We can also define a normalisation operation that returns a canonical representative for each
  association class. Since the four units of the Gaussian integers are ±1› and ±𝗂›, each
  association class (other than 0›) has four representatives, one in each quadrant. We simply
  define the on in the upper-right quadrant (i.e.\ the one with non-negative imaginary part
  and positive real part) as the canonical one.

  Thus, the Gaussian integers form a Euclidean ring. This gives us many things, most importantly
  the existence of GCDs and LCMs and unique factorisation.
›
instantiation gauss_int :: algebraic_semidom
begin

definition divide_gauss_int :: "gauss_int  gauss_int  gauss_int" where
  "divide_gauss_int a b = round_complex (gauss2complex a / gauss2complex b)"

instance proof
  fix a :: gauss_int
  show "a div 0 = 0"
    by (auto simp: gauss_int_eq_iff divide_gauss_int_def)
next
  fix a b :: gauss_int assume "b  0"
  thus "a * b div b = a"
    by (auto simp: gauss_int_eq_iff divide_gauss_int_def)
qed

end

instantiation gauss_int :: semidom_divide_unit_factor
begin

definition unit_factor_gauss_int :: "gauss_int  gauss_int" where
  "unit_factor_gauss_int z =
     (if z = 0 then 0 else 
      if ImZ z  0  ReZ z > 0 then 1
      else if ReZ z  0  ImZ z > 0 then 𝗂
      else if ImZ z  0  ReZ z < 0 then -1
      else -𝗂)"

instance proof
  show "unit_factor (0 :: gauss_int) = 0"
    by (simp add: unit_factor_gauss_int_def)
next
  fix z :: gauss_int
  assume "is_unit z"
  thus "unit_factor z = z"
    by (subst (asm) is_unit_gauss_int_iff) (auto simp: unit_factor_gauss_int_def)
next
  fix z :: gauss_int
  assume z: "z  0"
  thus "is_unit (unit_factor z)"
    by (subst is_unit_gauss_int_iff) (auto simp: unit_factor_gauss_int_def)
next
  fix z u :: gauss_int
  assume "is_unit z"
  hence "z  {1, -1, 𝗂, -𝗂}"
    by (subst (asm) is_unit_gauss_int_iff)
  thus "unit_factor (z * u) = z * unit_factor u"
    by (safe; auto simp: unit_factor_gauss_int_def gauss_int_eq_iff[of u 0])
qed

end

instantiation gauss_int :: normalization_semidom
begin

definition normalize_gauss_int :: "gauss_int  gauss_int" where
  "normalize_gauss_int z =
     (if z = 0 then 0 else 
      if ImZ z  0  ReZ z > 0 then z
      else if ReZ z  0  ImZ z > 0 then -𝗂 * z
      else if ImZ z  0  ReZ z < 0 then -z
      else 𝗂 * z)"

instance proof
  show "normalize (0 :: gauss_int) = 0"
    by (simp add: normalize_gauss_int_def)
next
  fix z :: gauss_int
  show "unit_factor z * normalize z = z"
    by (auto simp: normalize_gauss_int_def unit_factor_gauss_int_def algebra_simps)
qed

end

lemma normalize_gauss_int_of_nat [simp]: "normalize (of_nat n :: gauss_int) = of_nat n"
  and normalize_gauss_int_of_int [simp]: "normalize (of_int m :: gauss_int) = of_int ¦m¦"
  and normalize_gauss_int_of_numeral [simp]: "normalize (numeral n' :: gauss_int) = numeral n'"
  by (auto simp: normalize_gauss_int_def)

lemma normalize_gauss_i [simp]: "normalize 𝗂 = 1"
  by (simp add: normalize_gauss_int_def)

lemma gauss_int_norm_normalize [simp]: "gauss_int_norm (normalize x) = gauss_int_norm x"
  by (simp add: normalize_gauss_int_def gauss_int_norm_mult)

lemma normalized_gauss_int:
  assumes "normalize z = z"
  shows   "ReZ z  0" "ImZ z  0"
  using assms
  by (cases "ReZ z" "0 :: int" rule: linorder_cases;
      cases "ImZ z" "0 :: int" rule: linorder_cases;
      simp add: normalize_gauss_int_def gauss_int_eq_iff)+

lemma normalized_gauss_int':
  assumes "normalize z = z" "z  0"
  shows   "ReZ z > 0" "ImZ z  0"
  using assms
  by (cases "ReZ z" "0 :: int" rule: linorder_cases;
      cases "ImZ z" "0 :: int" rule: linorder_cases;
      simp add: normalize_gauss_int_def gauss_int_eq_iff)+

lemma normalized_gauss_int_iff:
  "normalize z = z  z = 0  ReZ z > 0  ImZ z  0"
  by (cases "ReZ z" "0 :: int" rule: linorder_cases;
      cases "ImZ z" "0 :: int" rule: linorder_cases;
      simp add: normalize_gauss_int_def gauss_int_eq_iff)+

instantiation gauss_int :: idom_modulo
begin

definition modulo_gauss_int :: "gauss_int  gauss_int  gauss_int" where
  "modulo_gauss_int a b = a - a div b * b"

instance proof
  fix a b :: gauss_int
  show "a div b * b + a mod b = a"
    by (simp add: modulo_gauss_int_def)
qed

end

lemma gauss_int_norm_mod_less_aux:
  assumes [simp]: "b  0"
  shows   "2 * gauss_int_norm (a mod b)  gauss_int_norm b"
proof -
  define a' b' where "a' = gauss2complex a" and "b' = gauss2complex b"
  have [simp]: "b'  0" by (simp add: b'_def)
  have "gauss_int_norm (a mod b) = 
          norm (gauss2complex (a - round_complex (a' / b') * b)) ^ 2"
    unfolding modulo_gauss_int_def
    by (subst real_gauss_int_norm [symmetric]) (auto simp add: divide_gauss_int_def a'_def b'_def)
  also have "gauss2complex (a - round_complex (a' / b') * b) =
               a' - gauss2complex (round_complex (a' / b')) * b'"
    by (simp add: a'_def b'_def)
  also have " = (a' / b' - gauss2complex (round_complex (a' / b'))) * b'"
    by (simp add: field_simps)
  also have "norm  ^ 2 = norm (a' / b' - gauss2complex (round_complex (a' / b'))) ^ 2 * norm b' ^ 2"
    by (simp add: norm_mult power_mult_distrib)
  also have "  1 / 2 * norm b' ^ 2"
    by (intro mult_right_mono norm_round_complex_le) auto
  also have "norm b' ^ 2 = gauss_int_norm b"
    by (simp add: b'_def real_gauss_int_norm)
  finally show ?thesis by linarith
qed

lemma gauss_int_norm_mod_less:
  assumes [simp]: "b  0"
  shows   "gauss_int_norm (a mod b) < gauss_int_norm b"
proof -
  have "gauss_int_norm b > 0" by simp
  thus "gauss_int_norm (a mod b) < gauss_int_norm b"
    using gauss_int_norm_mod_less_aux[OF assms, of a] by presburger
qed

lemma gauss_int_norm_dvd_imp_le:
  assumes "b  0"
  shows   "gauss_int_norm a  gauss_int_norm (a * b)"
proof (cases "a = 0")
  case False
  thus ?thesis using assms by (intro dvd_imp_le gauss_int_norm_dvd_mono) auto
qed auto

instantiation gauss_int :: euclidean_ring
begin

definition euclidean_size_gauss_int :: "gauss_int  nat" where
  [simp]: "euclidean_size_gauss_int = gauss_int_norm"

instance proof
  show "euclidean_size (0 :: gauss_int) = 0"
    by simp
next
  fix a b :: gauss_int assume [simp]: "b  0"
  show "euclidean_size (a mod b) < euclidean_size b"
    using gauss_int_norm_mod_less[of b a] by simp
  show "euclidean_size a  euclidean_size (a * b)"
    by (simp add: gauss_int_norm_dvd_imp_le)
qed

end

instance gauss_int :: normalization_euclidean_semiring ..

instantiation gauss_int :: euclidean_ring_gcd
begin

definition gcd_gauss_int :: "gauss_int  gauss_int  gauss_int" where
  "gcd_gauss_int  normalization_euclidean_semiring_class.gcd"
definition lcm_gauss_int :: "gauss_int  gauss_int  gauss_int" where
  "lcm_gauss_int  normalization_euclidean_semiring_class.lcm"
definition Gcd_gauss_int :: "gauss_int set  gauss_int" where
  "Gcd_gauss_int  normalization_euclidean_semiring_class.Gcd"
definition Lcm_gauss_int :: "gauss_int set  gauss_int" where
  "Lcm_gauss_int  normalization_euclidean_semiring_class.Lcm"

instance 
  by intro_classes
     (simp_all add: gcd_gauss_int_def lcm_gauss_int_def Gcd_gauss_int_def Lcm_gauss_int_def)

end

lemma multiplicity_gauss_cnj: "multiplicity (gauss_cnj a) (gauss_cnj b) = multiplicity a b"
  unfolding multiplicity_def gauss_cnj_power [symmetric] gauss_cnj_dvd_iff ..
    
lemma multiplicity_gauss_int_of_nat:
  "multiplicity (of_nat a) (of_nat b :: gauss_int) = multiplicity a b"
  unfolding multiplicity_def of_nat_power [symmetric] of_nat_dvd_of_nat_gauss_int_iff ..

lemma gauss_int_dvd_same_norm_imp_associated:
  assumes "z1 dvd z2" "gauss_int_norm z1 = gauss_int_norm z2"
  shows   "normalize z1 = normalize z2"
proof (cases "z1 = 0")
  case [simp]: False
  from assms(1) obtain u where u: "z2 = z1 * u" by blast
  from assms have "gauss_int_norm u = 1"
    by (auto simp: gauss_int_norm_mult u)
  hence "is_unit u"
    by (simp add: is_unit_gauss_int_iff')
  with u show ?thesis by simp
qed (use assms in auto)

lemma gcd_of_int_gauss_int: "gcd (of_int a :: gauss_int) (of_int b) = of_int (gcd a b)"
proof (induction "nat ¦b¦" arbitrary: a b rule: less_induct)
  case (less b a)
  show ?case
  proof (cases "b = 0")
    case False
    have "of_int (gcd a b) = (of_int (gcd b (a mod b)) :: gauss_int)"
      by (subst gcd_red_int) auto
    also have " = gcd (of_int b) (of_int (a mod b))"
      using False by (intro less [symmetric]) (auto intro!: abs_mod_less)
    also have "a mod b = (a - a div b * b)"
      by (simp add: minus_div_mult_eq_mod)
    also have "of_int  = of_int (-(a div b)) * of_int b + (of_int a :: gauss_int)"
      by (simp add: algebra_simps)
    also have "gcd (of_int b)  = gcd (of_int b) (of_int a)"
      by (rule gcd_add_mult)
    finally show ?thesis by (simp add: gcd.commute)
  qed auto
qed

lemma coprime_of_int_gauss_int: "coprime (of_int a :: gauss_int) (of_int b) = coprime a b"
  unfolding coprime_iff_gcd_eq_1 gcd_of_int_gauss_int by auto

lemma gcd_of_nat_gauss_int: "gcd (of_nat a :: gauss_int) (of_nat b) = of_nat (gcd a b)"
  using gcd_of_int_gauss_int[of "int a" "int b"] by simp

lemma coprime_of_nat_gauss_int: "coprime (of_nat a :: gauss_int) (of_nat b) = coprime a b"
  unfolding coprime_iff_gcd_eq_1 gcd_of_nat_gauss_int by auto

lemma gauss_cnj_dvd_self_iff: "gauss_cnj z dvd z  ReZ z = 0  ImZ z = 0  ¦ReZ z¦ = ¦ImZ z¦"
proof
  assume "gauss_cnj z dvd z"
  hence "normalize (gauss_cnj z) = normalize z"
    by (rule gauss_int_dvd_same_norm_imp_associated) auto
  then obtain u :: gauss_int where "is_unit u"  and u: "gauss_cnj z = u * z"
    using associatedE1 by blast
  hence "u  {1, -1, 𝗂, -𝗂}"
    by (simp add: is_unit_gauss_int_iff)
  thus "ReZ z = 0  ImZ z = 0  ¦ReZ z¦ = ¦ImZ z¦"
  proof (elim insertE emptyE)
    assume [simp]: "u = 𝗂"
    have "ReZ z = ReZ (gauss_cnj z)"
      by simp
    also have "gauss_cnj z = 𝗂 * z"
      using u by simp
    also have "ReZ  = -ImZ z"
      by simp
    finally show "ReZ z = 0  ImZ z = 0  ¦ReZ z¦ = ¦ImZ z¦"
      by auto
  next
    assume [simp]: "u = -𝗂"
    have "ReZ z = ReZ (gauss_cnj z)"
      by simp
    also have "gauss_cnj z = -𝗂 * z"
      using u by simp
    also have "ReZ  = ImZ z"
      by simp
    finally show "ReZ z = 0  ImZ z = 0  ¦ReZ z¦ = ¦ImZ z¦"
      by auto
  next
    assume [simp]: "u = 1"
    have "ImZ z = -ImZ (gauss_cnj z)"
      by simp
    also have "gauss_cnj z = z"
      using u by simp
    finally show "ReZ z = 0  ImZ z = 0  ¦ReZ z¦ = ¦ImZ z¦"
      by auto
  next
    assume [simp]: "u = -1"
    have "ReZ z = ReZ (gauss_cnj z)"
      by simp
    also have "gauss_cnj z = -z"
      using u by simp
    also have "ReZ  = -ReZ z"
      by simp
    finally show "ReZ z = 0  ImZ z = 0  ¦ReZ z¦ = ¦ImZ z¦"
      by auto
  qed
next
  assume "ReZ z = 0  ImZ z = 0  ¦ReZ z¦ = ¦ImZ z¦"
  thus "gauss_cnj z dvd z"
  proof safe
    assume "¦ReZ z¦ = ¦ImZ z¦"
    then obtain u :: int where "is_unit u" and u: "ImZ z = u * ReZ z"
      using associatedE2[of "ReZ z" "ImZ z"] by auto
    from is_unit u have "u  {1, -1}"
      by auto
    hence "z = gauss_cnj z * (of_int u * 𝗂)"
      using u by (auto simp: gauss_int_eq_iff)
    thus ?thesis
      by (metis dvd_triv_left)
  qed (auto simp: gauss_cnj_eq_self gauss_cnj_eq_minus_self)
qed

lemma self_dvd_gauss_cnj_iff: "z dvd gauss_cnj z  ReZ z = 0  ImZ z = 0  ¦ReZ z¦ = ¦ImZ z¦"
  using gauss_cnj_dvd_self_iff[of z] by (subst (asm) gauss_cnj_dvd_left_iff) auto


subsection ‹Prime elements›

text ‹
  Next, we analyse what the prime elements of the Gaussian integers are. First, note that
  according to the conventions of Isabelle's computational algebra library, a prime element
  is called a prime iff it is also normalised, i.e.\ in our case it lies in the upper right
  quadrant.

  As a first fact, we can show that a Gaussian integer whose norm is ℤ›-prime must be
  $\mathbb{Z}[i]$-prime:
›

lemma prime_gauss_int_norm_imp_prime_elem:
  assumes "prime (gauss_int_norm q)"
  shows   "prime_elem q"
proof -
  have "irreducible q"
  proof (rule irreducibleI)
    fix a b assume "q = a * b"
    hence "gauss_int_norm q = gauss_int_norm a * gauss_int_norm b"
      by (simp_all add: gauss_int_norm_mult)
    thus "is_unit a  is_unit b"
      using assms by (auto dest!: prime_product simp: gauss_int_norm_eq_Suc_0_iff)
  qed (use assms in auto simp: is_unit_gauss_int_iff')
  thus "prime_elem q"
    using irreducible_imp_prime_elem_gcd by blast
qed

text ‹
  Also, a conjugate is a prime element iff the original element is a prime element:
›
lemma prime_elem_gauss_cnj [intro]: "prime_elem z  prime_elem (gauss_cnj z)"
  by (auto simp: prime_elem_def gauss_cnj_dvd_left_iff)

lemma prime_elem_gauss_cnj_iff [simp]: "prime_elem (gauss_cnj z)  prime_elem z"
  using prime_elem_gauss_cnj[of z] prime_elem_gauss_cnj[of "gauss_cnj z"] by auto


subsubsection ‹The factorisation of 2›

text ‹
  2 factors as $-i (1 + i)^2$ in the Gaussian integers, where $-i$ is a unit and
  $1 + i$ is prime.
›

lemma gauss_int_2_eq: "2 = -𝗂 * (1 + 𝗂) ^ 2"
  by (simp add: gauss_int_eq_iff power2_eq_square)

lemma prime_elem_one_plus_i_gauss_int: "prime_elem (1 + 𝗂)"
  by (rule prime_gauss_int_norm_imp_prime_elem) (auto simp: gauss_int_norm_def)

lemma prime_one_plus_i_gauss_int: "prime (1 + 𝗂)"
  by (simp add: prime_def prime_elem_one_plus_i_gauss_int
                gauss_int_eq_iff normalize_gauss_int_def)

lemma prime_factorization_2_gauss_int:
  "prime_factorization (2 :: gauss_int) = {#1 + 𝗂, 1 + 𝗂#}"
proof -
  have "prime_factorization (2 :: gauss_int) =
        (prime_factorization (prod_mset {#1 + gauss_i, 1 + gauss_i#}))"
    by (subst prime_factorization_unique) (auto simp: gauss_int_eq_iff normalize_gauss_int_def)
  also have "prime_factorization (prod_mset {#1 + gauss_i, 1 + gauss_i#}) =
               {#1 + gauss_i, 1 + gauss_i#}"
    using prime_one_plus_i_gauss_int by (subst prime_factorization_prod_mset_primes) auto
  finally show ?thesis .
qed


subsubsection ‹Inert primes›

text ‹
  Any ℤ›-prime congruent 3 modulo 4 is also a Gaussian prime. These primes are called
  ‹inert›, because they do not decompose when moving from ℤ› to $\mathbb{Z}[i]$.
›

lemma gauss_int_norm_not_3_mod_4: "[gauss_int_norm z  3] (mod 4)"
proof -
  have A: "ReZ z mod 4  {0..3}" "ImZ z mod 4  {0..3}" by auto
  have B: "{0..3} = {0, 1, 2, 3 :: int}" by auto

  have "[ReZ z ^ 2 + ImZ z ^ 2 = (ReZ z mod 4) ^ 2 + (ImZ z mod 4) ^ 2] (mod 4)"
    by (intro cong_add cong_pow) (auto simp: cong_def)
  moreover have "((ReZ z mod 4) ^ 2 + (ImZ z mod 4) ^ 2) mod 4  3 mod 4"
    using A unfolding B by auto
  ultimately have "[ReZ z ^ 2 + ImZ z ^ 2  3] (mod 4)"
    unfolding cong_def by metis
  hence "[int (nat (ReZ z ^ 2 + ImZ z ^ 2))  int 3] (mod (int 4))"
    by simp
  thus ?thesis unfolding gauss_int_norm_def
    by (subst (asm) cong_int_iff)
qed

lemma prime_elem_gauss_int_of_nat:
  fixes n :: nat
  assumes prime: "prime n" and "[n = 3] (mod 4)"
  shows   "prime_elem (of_nat n :: gauss_int)"
proof (intro irreducible_imp_prime_elem irreducibleI)
  from assms show "of_nat n  (0 :: gauss_int)"
    by (auto simp: gauss_int_eq_iff)
next
  show "¬is_unit (of_nat n :: gauss_int)"
    using assms by (subst is_unit_gauss_int_iff) (auto simp: gauss_int_eq_iff)
next
  fix a b :: gauss_int
  assume *: "of_nat n = a * b"
  hence "gauss_int_norm (a * b) = gauss_int_norm (of_nat n)"
    by metis
  hence *: "gauss_int_norm a * gauss_int_norm b = n ^ 2"
    by (simp add: gauss_int_norm_mult power2_eq_square flip: nat_mult_distrib)
  from prime_power_mult_nat[OF prime this] obtain i j :: nat
    where ij: "gauss_int_norm a = n ^ i" "gauss_int_norm b = n ^ j" by blast
  
  have "i + j = 2"
  proof -
    have "n ^ (i + j) = n ^ 2"
      using ij * by (simp add: power_add)
    from prime_power_inj[OF prime this] show ?thesis by simp
  qed
  hence "i = 0  j = 2  i = 1  j = 1  i = 2  j = 0"
    by auto
  thus "is_unit a  is_unit b"
  proof (elim disjE)
    assume "i = 1  j = 1"
    with ij have "gauss_int_norm a = n"
      by auto
    hence "[gauss_int_norm a = n] (mod 4)"
      by simp
    also have "[n = 3] (mod 4)" by fact
    finally have "[gauss_int_norm a = 3] (mod 4)" .
    moreover have "[gauss_int_norm a  3] (mod 4)"
      by (rule gauss_int_norm_not_3_mod_4)
    ultimately show ?thesis by contradiction
  qed (use ij in auto simp: is_unit_gauss_int_iff')
qed

theorem prime_gauss_int_of_nat:
  fixes n :: nat
  assumes prime: "prime n" and "[n = 3] (mod 4)"
  shows   "prime (of_nat n :: gauss_int)"
  using prime_elem_gauss_int_of_nat[OF assms]
  unfolding prime_def by simp


subsubsection ‹Non-inert primes›

text ‹
  Any ℤ›-prime congruent 1 modulo 4 factors into two conjugate Gaussian primes.
›

lemma minimal_QuadRes_neg1:
  assumes "QuadRes n (-1)" "n > 1" "odd n"
  obtains x :: nat where "x  (n - 1) div 2" and "[x ^ 2 + 1 = 0] (mod n)"
proof -
  from QuadRes n (-1) obtain x where "[x ^ 2 = (-1)] (mod (int n))"
    by (auto simp: QuadRes_def)
  hence "[x ^ 2 + 1 = -1 + 1] (mod (int n))"
    by (intro cong_add) auto
  also have "x ^ 2 + 1 = int (nat ¦x¦ ^ 2 + 1)"
    by simp
  finally have "[int (nat ¦x¦ ^ 2 + 1) = int 0] (mod (int n))"
    by simp
  hence "[nat ¦x¦ ^ 2 + 1 = 0] (mod n)"
    by (subst (asm) cong_int_iff)

  define x' where
    "x' = (if nat ¦x¦ mod n  (n - 1) div 2 then nat ¦x¦ mod n else n - (nat ¦x¦ mod n))"
  have x'_quadres: "[x' ^ 2 + 1 = 0] (mod n)"
  proof (cases "nat ¦x¦ mod n  (n - 1) div 2")
    case True
    hence "[x' ^ 2 + 1 = (nat ¦x¦ mod n) ^ 2 + 1] (mod n)"
      by (simp add: x'_def)
    also have "[(nat ¦x¦ mod n) ^ 2 + 1 = nat ¦x¦ ^ 2 + 1] (mod n)"
      by (intro cong_add cong_pow) (auto simp: cong_def)
    also have "[nat ¦x¦ ^ 2 + 1 = 0] (mod n)" by fact
    finally show ?thesis .
  next
    case False
    hence "[int (x' ^ 2 + 1) = (int n - int (nat ¦x¦ mod n)) ^ 2 + 1] (mod int n)"
      using n > 1 by (simp add: x'_def of_nat_diff add_ac)
    also have "[(int n - int (nat ¦x¦ mod n)) ^ 2 + 1 =
                (0 - int (nat ¦x¦ mod n)) ^ 2 + 1] (mod int n)"
      by (intro cong_add cong_pow) (auto simp: cong_def)
    also have "[(0 - int (nat ¦x¦ mod n)) ^ 2 + 1 = int ((nat ¦x¦ mod n) ^ 2 + 1)] (mod (int n))"
      by (simp add: add_ac)
    finally have "[x' ^ 2 + 1 = (nat ¦x¦ mod n)2 + 1] (mod n)"
      by (subst (asm) cong_int_iff)
    also have "[(nat ¦x¦ mod n)2 + 1 = nat ¦x¦ ^ 2 + 1] (mod n)"
      by (intro cong_add cong_pow) (auto simp: cong_def)
    also have "[nat ¦x¦ ^ 2 + 1 = 0] (mod n)" by fact
    finally show ?thesis .
  qed
  moreover have x'_le: "x'  (n - 1) div 2"
    using odd n by (auto elim!: oddE simp: x'_def)
  ultimately show ?thesis by (intro that[of x'])
qed

text ‹
  Let p› be some prime number that is congruent 1 modulo 4.
›
locale noninert_gauss_int_prime =
  fixes p :: nat
  assumes prime_p: "prime p" and cong_1_p: "[p = 1] (mod 4)"
begin

lemma p_gt_2: "p > 2" and odd_p: "odd p"
proof -
  from prime_p and cong_1_p have "p > 1" "p  2"
    by (auto simp: prime_gt_Suc_0_nat cong_def)
  thus "p > 2" by auto
  with prime_p show "odd p"
    using primes_dvd_imp_eq two_is_prime_nat by blast
qed

text ‹
  -1 is a quadratic residue modulo p›, so there exists some x› such that
  $x^2 + 1$ is divisible by p›. Moreover, we can choose x› such that it is positive and
  no greater than $\frac{1}{2}(p-1)$:
›
lemma minimal_QuadRes_neg1:
  obtains x where "x > 0" "x  (p - 1) div 2" "[x ^ 2 + 1 = 0] (mod p)"
proof -
  have "[Legendre (-1) (int p) = (- 1) ^ ((p - 1) div 2)] (mod (int p))"
    using prime_p p_gt_2 by (intro euler_criterion) auto
  also have "[p - 1 = 1 - 1] (mod 4)"
    using p_gt_2 by (intro cong_diff_nat cong_refl) (use cong_1_p in auto)
  hence "2 * 2 dvd p - 1"
    by (simp add: cong_0_iff)
  hence "even ((p - 1) div 2)"
    using dvd_mult_imp_div by blast
  hence "(-1) ^ ((p - 1) div 2) = (1 :: int)"
    by simp
  finally have "Legendre (-1) (int p) mod p = 1"
    using p_gt_2 by (auto simp: cong_def)
  hence "Legendre (-1) (int p) = 1"
    using p_gt_2 by (auto simp: Legendre_def cong_def zmod_minus1 split: if_splits)
  hence "QuadRes p (-1)"
    by (simp add: Legendre_def split: if_splits)
  from minimal_QuadRes_neg1[OF this] p_gt_2 odd_p
    obtain x where x: "x  (p - 1) div 2" "[x ^ 2 + 1 = 0] (mod p)" by auto
  have "x > 0"
    using x p_gt_2 by (auto intro!: Nat.gr0I simp: cong_def)
  from x and this show ?thesis by (intro that[of x]) auto
qed

text ‹
  We can show from this that p› is not prime as a Gaussian integer.
›
lemma not_prime: "¬prime_elem (of_nat p :: gauss_int)"
proof
  assume prime: "prime_elem (of_nat p :: gauss_int)"
  obtain x where x: "x > 0" "x  (p - 1) div 2" "[x ^ 2 + 1 = 0] (mod p)"
    using  minimal_QuadRes_neg1 .

  have "of_nat p dvd (of_nat (x ^ 2 + 1) :: gauss_int)"
    using x by (intro of_nat_dvd_of_nat) (auto simp: cong_0_iff)
  also have eq: "of_nat (x ^ 2 + 1) = ((of_nat x + 𝗂) * (of_nat x - 𝗂) :: gauss_int)"
    using x > 0 by (simp add: algebra_simps gauss_int_eq_iff power2_eq_square of_nat_diff)
  finally have "of_nat p dvd ((of_nat x + 𝗂) * (of_nat x - 𝗂) :: gauss_int)" .

  from prime and this
    have "of_nat p dvd (of_nat x + 𝗂 :: gauss_int)  of_nat p dvd (of_nat x - 𝗂 :: gauss_int)"
    by (rule prime_elem_dvd_multD)
  hence dvd: "of_nat p dvd (of_nat x + 𝗂 :: gauss_int)" "of_nat p dvd (of_nat x - 𝗂 :: gauss_int)"
    by (auto dest: of_nat_dvd_imp_dvd_gauss_cnj)

  have "of_nat (p ^ 2) = (of_nat p * of_nat p :: gauss_int)"
    by (simp add: power2_eq_square)
  also from dvd have " dvd ((of_nat x + 𝗂) * (of_nat x - 𝗂))"
    by (intro mult_dvd_mono)
  also have " = of_nat (x ^ 2 + 1)"
    by (rule eq [symmetric])
  finally have "p ^ 2 dvd (x ^ 2 + 1)"
    by (subst (asm) of_nat_dvd_of_nat_gauss_int_iff)
  hence "p ^ 2  x ^ 2 + 1"
    by (intro dvd_imp_le) auto
  moreover have "p ^ 2 > x ^ 2 + 1"
  proof -
    have "x ^ 2 + 1  ((p - 1) div 2) ^ 2 + 1"
      using x by (intro add_mono power_mono) auto
    also have "  (p - 1) ^ 2 + 1"
      by auto
    also have "(p - 1) * (p - 1) < (p - 1) * (p + 1)"
      using p_gt_2 by (intro mult_strict_left_mono) auto
    hence "(p - 1) ^ 2 + 1 < p ^ 2"
      by (simp add: algebra_simps power2_eq_square)
    finally show ?thesis .
  qed
  ultimately show False by linarith
qed

text ‹
  Any prime factor of p› in the Gaussian integers must have norm p›.
›
lemma norm_prime_divisor:
  fixes q :: gauss_int
  assumes q: "prime_elem q" "q dvd of_nat p"
  shows "gauss_int_norm q = p"
proof -
  from assms obtain r where r: "of_nat p = q * r"
    by auto
  have "p ^ 2 = gauss_int_norm (of_nat p)"
    by simp
  also have " = gauss_int_norm q * gauss_int_norm r"
    by (auto simp: r gauss_int_norm_mult)
  finally have *: "gauss_int_norm q * gauss_int_norm r = p ^ 2"
    by simp
  hence "i j. gauss_int_norm q = p ^ i  gauss_int_norm r = p ^ j"
    using prime_p by (intro prime_power_mult_nat)
  then obtain i j where ij: "gauss_int_norm q = p ^ i" "gauss_int_norm r = p ^ j"
    by blast
  have ij_eq_2: "i + j = 2"
  proof -
    from * have "p ^ (i