# 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)"
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"
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"
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"
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:
fixes x :: int
assumes "even x"
shows   "[x ^ 2 = 0] (mod 4)"
proof -
from assms have "even ¦x¦"
by simp
hence [simp]: "¦x¦ mod 2 = 0"
by presburger
have "(¦x¦ ^ 2) mod 4 = ((¦x¦ mod 4) ^ 2) mod 4"
also from assms have "¦x¦ mod 4 = 0 ∨ ¦x¦ mod 4 = 2"
using mod_double_modulus[of 2 "¦x¦"] by simp
hence "((¦x¦ mod 4) ^ 2) mod 4 = 0"
by auto
finally show ?thesis by (simp add: cong_def)
qed

lemma even_square_cong_4_nat: "even (x::nat) ⟹ [x ^ 2 = 0] (mod 4)"
using even_square_cong_4_int[of "int x"] by (auto simp flip: cong_int_iff)

lemma odd_square_cong_4_int:
fixes x :: int
assumes "odd x"
shows   "[x ^ 2 = 1] (mod 4)"
proof -
from assms have "odd ¦x¦"
by simp
hence [simp]: "¦x¦ mod 2 = 1"
by presburger
have "(¦x¦ ^ 2) mod 4 = ((¦x¦ mod 4) ^ 2) mod 4"
also from assms have "¦x¦ mod 4 = 1 ∨ ¦x¦ mod 4 = 3"
using mod_double_modulus[of 2 "¦x¦"] by simp
hence "((¦x¦ mod 4) ^ 2) mod 4 = 1"
by auto
finally show ?thesis by (simp add: cong_def)
qed

lemma odd_square_cong_4_nat: "odd (x::nat) ⟹ [x ^ 2 = 1] (mod 4)"
using odd_square_cong_4_int[of "int x"] by (auto 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"
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 ((∏p∈prime_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 ((∏p∈prime_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
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]

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"

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 𝗂⇩ℤ = -𝗂⇩ℤ"

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)"

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

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"

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"

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) = (∑x∈A. gauss_cnj (f x))"
by (induction A rule: infinite_finite_induct) auto

lemma gauss_cnj_prod [simp]: "gauss_cnj (prod f A) = (∏x∈A. 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 * 𝗂⇩ℤ"

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')"

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

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"
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"

lemma gauss_int_norm_eq_Suc_0_iff: "gauss_int_norm x = Suc 0 ⟷ x dvd 1"

lemma is_unit_gauss_cnj [intro]: "z dvd 1 ⟹ gauss_cnj z dvd 1"

lemma is_unit_gauss_cnj_iff [simp]: "gauss_cnj z dvd 1 ⟷ z dvd 1"

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"
thus "norm (z - gauss2complex (round_complex z)) ^ 2 ≤ 1 / 2"
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"
also have "… ≤ 1 / 2"
by (rule norm_round_complex_le)
also have "… = (sqrt 2 / 2) ^ 2"
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"
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"
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"

lemma gauss_int_norm_normalize [simp]: "gauss_int_norm (normalize x) = gauss_int_norm x"

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;

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;

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;

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"
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'"
also have "… = (a' / b' - gauss2complex (round_complex (a' / b'))) * b'"
also have "norm … ^ 2 = norm (a' / b' - gauss2complex (round_complex (a' / b'))) ^ 2 * norm b' ^ 2"
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"
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)"
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"
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)"
also have "of_int … = of_int (-(a div b)) * of_int b + (of_int a :: gauss_int)"
also have "gcd (of_int b) … = gcd (of_int b) (of_int a)"
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, 𝗂⇩ℤ, -𝗂⇩ℤ}"
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

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"
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"

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 + 𝗂⇩ℤ)"
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"
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)
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.
›

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))"
hence "[x ^ 2 + 1 = -1 + 1] (mod (int n))"
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)"
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)"
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))"
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)$:
›
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"
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)
by (simp add: Legendre_def split: if_splits)
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)"

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)"
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"
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: "`