Theory HOL-Examples.Gauss_Numbers

(*   Author:      Florian Haftmann, TU Muenchen; based on existing material on complex numbers›
*)

section ‹Gauss Numbers: integral gauss numbers›

theory Gauss_Numbers
  imports "HOL-Library.Centered_Division"
begin

codatatype gauss = Gauss (Re: int) (Im: int)

lemma gauss_eqI [intro?]:
  ‹x = y› if ‹Re x = Re y› ‹Im x = Im y›
  by (rule gauss.expand) (use that in simp)

lemma gauss_eq_iff:
  ‹x = y ⟷ Re x = Re y ∧ Im x = Im y›
  by (auto intro: gauss_eqI)


subsection ‹Basic arithmetic›

instantiation gauss :: comm_ring_1
begin

primcorec zero_gauss :: ‹gauss›
  where
    ‹Re 0 = 0›
  | ‹Im 0 = 0›

primcorec one_gauss :: ‹gauss›
  where
    ‹Re 1 = 1›
  | ‹Im 1 = 0›

primcorec plus_gauss :: ‹gauss ⇒ gauss ⇒ gauss›
  where
    ‹Re (x + y) = Re x + Re y›
  | ‹Im (x + y) = Im x + Im y›

primcorec uminus_gauss :: ‹gauss ⇒ gauss›
  where
    ‹Re (- x) = - Re x›
  | ‹Im (- x) = - Im x›

primcorec minus_gauss :: ‹gauss ⇒ gauss ⇒ gauss›
  where
    ‹Re (x - y) = Re x - Re y›
  | ‹Im (x - y) = Im x - Im y›

primcorec times_gauss :: ‹gauss ⇒ gauss ⇒ gauss›
  where
    ‹Re (x * y) = Re x * Re y - Im x * Im y›
  | ‹Im (x * y) = Re x * Im y + Im x * Re y›

instance
  by standard (simp_all add: gauss_eq_iff algebra_simps)

end

lemma of_nat_gauss:
  ‹of_nat n = Gauss (int n) 0›
  by (induction n) (simp_all add: gauss_eq_iff)

lemma numeral_gauss:
  ‹numeral n = Gauss (numeral n) 0›
proof -
  have ‹numeral n = (of_nat (numeral n) :: gauss)›
    by simp
  also have ‹… = Gauss (of_nat (numeral n)) 0›
    by (simp add: of_nat_gauss)
  finally show ?thesis
    by simp
qed

lemma of_int_gauss:
  ‹of_int k = Gauss k 0›
  by (simp add: gauss_eq_iff of_int_of_nat of_nat_gauss)

lemma conversion_simps [simp]:
  ‹Re (numeral m) = numeral m›
  ‹Im (numeral m) = 0›
  ‹Re (of_nat n) = int n›
  ‹Im (of_nat n) = 0›
  ‹Re (of_int k) = k›
  ‹Im (of_int k) = 0›
  by (simp_all add: numeral_gauss of_nat_gauss of_int_gauss)

lemma gauss_eq_0:
  ‹z = 0 ⟷ (Re z)⇧2 + (Im z)⇧2 = 0›
  by (simp add: gauss_eq_iff sum_power2_eq_zero_iff)

lemma gauss_neq_0:
  ‹z ≠ 0 ⟷ (Re z)⇧2 + (Im z)⇧2 > 0›
  by (simp add: gauss_eq_0 sum_power2_ge_zero less_le)

lemma Re_sum [simp]:
  ‹Re (sum f s) = (∑x∈s. Re (f x))›
  by (induct s rule: infinite_finite_induct) auto

lemma Im_sum [simp]:
  ‹Im (sum f s) = (∑x∈s. Im (f x))›
  by (induct s rule: infinite_finite_induct) auto

instance gauss :: idom
proof
  fix x y :: gauss
  assume ‹x ≠ 0› ‹y ≠ 0›
  then show ‹x * y ≠ 0›
    by (simp_all add: gauss_eq_iff)
      (smt (verit, best) mult_eq_0_iff mult_neg_neg mult_neg_pos mult_pos_neg mult_pos_pos)
qed



subsection ‹The Gauss Number $i$›

primcorec imaginary_unit :: gauss  (‹𝗂›)
  where
    ‹Re 𝗂 = 0›
  | ‹Im 𝗂 = 1›

lemma Gauss_eq:
  ‹Gauss a b = of_int a + 𝗂 * of_int b›
  by (simp add: gauss_eq_iff)

lemma gauss_eq:
  ‹a = of_int (Re a) + 𝗂 * of_int (Im a)›
  by (simp add: gauss_eq_iff)

lemma gauss_i_not_zero [simp]:
  ‹𝗂 ≠ 0›
  by (simp add: gauss_eq_iff)

lemma gauss_i_not_one [simp]:
  ‹𝗂 ≠ 1›
  by (simp add: gauss_eq_iff)

lemma gauss_i_not_numeral [simp]:
  ‹𝗂 ≠ numeral n›
  by (simp add: gauss_eq_iff)

lemma gauss_i_not_neg_numeral [simp]:
  ‹𝗂 ≠ - numeral n›
  by (simp add: gauss_eq_iff)

lemma i_mult_i_eq [simp]:
  ‹𝗂 * 𝗂 = - 1›
  by (simp add: gauss_eq_iff)

lemma gauss_i_mult_minus [simp]:
  ‹𝗂 * (𝗂 * x) = - x›
  by (simp flip: mult.assoc)

lemma i_squared [simp]:
  ‹𝗂⇧2 = - 1›
  by (simp add: power2_eq_square)

lemma i_even_power [simp]:
  ‹𝗂 ^ (n * 2) = (- 1) ^ n›
  unfolding mult.commute [of n] power_mult by simp

lemma Re_i_times [simp]:
  ‹Re (𝗂 * z) = - Im z›
  by simp

lemma Im_i_times [simp]:
  ‹Im (𝗂 * z) = Re z›
  by simp

lemma i_times_eq_iff:
  ‹𝗂 * w = z ⟷ w = - (𝗂 * z)›
  by auto

lemma is_unit_i [simp]:
  ‹𝗂 dvd 1›
  by (rule dvdI [of _ _ ‹- 𝗂›]) simp

lemma gauss_numeral [code_post]:
  ‹Gauss 0 0 = 0›
  ‹Gauss 1 0 = 1›
  ‹Gauss (- 1) 0 = - 1›
  ‹Gauss (numeral n) 0 = numeral n›
  ‹Gauss (- numeral n) 0 = - numeral n›
  ‹Gauss 0 1 = 𝗂›
  ‹Gauss 0 (- 1) = - 𝗂›
  ‹Gauss 0 (numeral n) = numeral n * 𝗂›
  ‹Gauss 0 (- numeral n) = - numeral n * 𝗂›
  ‹Gauss 1 1 = 1 + 𝗂›
  ‹Gauss (- 1) 1 = - 1 + 𝗂›
  ‹Gauss (numeral n) 1 = numeral n + 𝗂›
  ‹Gauss (- numeral n) 1 = - numeral n + 𝗂›
  ‹Gauss 1 (- 1) = 1 - 𝗂›
  ‹Gauss 1 (numeral n) = 1 + numeral n * 𝗂›
  ‹Gauss 1 (- numeral n) = 1 - numeral n * 𝗂›
  ‹Gauss (- 1) (- 1) = - 1 - 𝗂›
  ‹Gauss (numeral n) (- 1) = numeral n - 𝗂›
  ‹Gauss (- numeral n) (- 1) = - numeral n - 𝗂›
  ‹Gauss (- 1) (numeral n) = - 1 + numeral n * 𝗂›
  ‹Gauss (- 1) (- numeral n) = - 1 - numeral n * 𝗂›
  ‹Gauss (numeral m) (numeral n) = numeral m + numeral n * 𝗂›
  ‹Gauss (- numeral m) (numeral n) = - numeral m + numeral n * 𝗂›
  ‹Gauss (numeral m) (- numeral n) = numeral m - numeral n * 𝗂›
  ‹Gauss (- numeral m) (- numeral n) = - numeral m - numeral n * 𝗂›
  by (simp_all add: gauss_eq_iff)


subsection ‹Gauss Conjugation›

primcorec cnj :: ‹gauss ⇒ gauss›
  where
    ‹Re (cnj z) = Re z›
  | ‹Im (cnj z) = - Im z›

lemma gauss_cnj_cancel_iff [simp]:
  ‹cnj x = cnj y ⟷ x = y›
  by (simp add: gauss_eq_iff)

lemma gauss_cnj_cnj [simp]:
  ‹cnj (cnj z) = z›
  by (simp add: gauss_eq_iff)

lemma gauss_cnj_zero [simp]:
  ‹cnj 0 = 0›
  by (simp add: gauss_eq_iff)

lemma gauss_cnj_zero_iff [iff]:
  ‹cnj z = 0 ⟷ z = 0›
  by (simp add: gauss_eq_iff)

lemma gauss_cnj_one_iff [simp]:
  ‹cnj z = 1 ⟷ z = 1›
  by (simp add: gauss_eq_iff)

lemma gauss_cnj_add [simp]:
  ‹cnj (x + y) = cnj x + cnj y›
  by (simp add: gauss_eq_iff)

lemma cnj_sum [simp]:
  ‹cnj (sum f s) = (∑x∈s. cnj (f x))›
  by (induct s rule: infinite_finite_induct) auto

lemma gauss_cnj_diff [simp]:
  ‹cnj (x - y) = cnj x - cnj y›
  by (simp add: gauss_eq_iff)

lemma gauss_cnj_minus [simp]:
  ‹cnj (- x) = - cnj x›
  by (simp add: gauss_eq_iff)

lemma gauss_cnj_one [simp]:
  ‹cnj 1 = 1›
  by (simp add: gauss_eq_iff)

lemma gauss_cnj_mult [simp]:
  ‹cnj (x * y) = cnj x * cnj y›
  by (simp add: gauss_eq_iff)

lemma cnj_prod [simp]:
  ‹cnj (prod f s) = (∏x∈s. cnj (f x))›
  by (induct s rule: infinite_finite_induct) auto

lemma gauss_cnj_power [simp]:
  ‹cnj (x ^ n) = cnj x ^ n›
  by (induct n) simp_all

lemma gauss_cnj_numeral [simp]:
  ‹cnj (numeral w) = numeral w›
  by (simp add: gauss_eq_iff)

lemma gauss_cnj_of_nat [simp]:
  ‹cnj (of_nat n) = of_nat n›
  by (simp add: gauss_eq_iff)

lemma gauss_cnj_of_int [simp]:
  ‹cnj (of_int z) = of_int z›
  by (simp add: gauss_eq_iff)

lemma gauss_cnj_i [simp]:
  ‹cnj 𝗂 = - 𝗂›
  by (simp add: gauss_eq_iff)

lemma gauss_add_cnj:
  ‹z + cnj z = of_int (2 * Re z)›
  by (simp add: gauss_eq_iff)

lemma gauss_diff_cnj:
  ‹z - cnj z = of_int (2 * Im z) * 𝗂›
  by (simp add: gauss_eq_iff)

lemma gauss_mult_cnj:
  ‹z * cnj z = of_int ((Re z)⇧2 + (Im z)⇧2)›
  by (simp add: gauss_eq_iff power2_eq_square)

lemma cnj_add_mult_eq_Re:
  ‹z * cnj w + cnj z * w = of_int (2 * Re (z * cnj w))›
  by (simp add: gauss_eq_iff)

lemma gauss_In_mult_cnj_zero [simp]:
  ‹Im (z * cnj z) = 0›
  by simp


subsection ‹Algebraic division›

instantiation gauss :: idom_modulo
begin

primcorec divide_gauss :: ‹gauss ⇒ gauss ⇒ gauss›
  where
    ‹Re (x div y) = (Re x * Re y + Im x * Im y) cdiv ((Re y)⇧2 + (Im y)⇧2)›
  | ‹Im (x div y) = (Im x * Re y - Re x * Im y) cdiv ((Re y)⇧2 + (Im y)⇧2)›

primcorec modulo_gauss :: ‹gauss ⇒ gauss ⇒ gauss›
  where
    ‹Re (x mod y) = Re x -
      ((Re x * Re y + Im x * Im y) cdiv ((Re y)⇧2 + (Im y)⇧2) * Re y -
       (Im x * Re y - Re x * Im y) cdiv ((Re y)⇧2 + (Im y)⇧2) * Im y)›
  | ‹Im (x mod y) = Im x -
      ((Re x * Re y + Im x * Im y) cdiv ((Re y)⇧2 + (Im y)⇧2) * Im y +
       (Im x * Re y - Re x * Im y) cdiv ((Re y)⇧2 + (Im y)⇧2) * Re y)›

instance proof
  fix x y :: gauss
  show ‹x div 0 = 0›
    by (simp add: gauss_eq_iff)
  show ‹x * y div y = x› if ‹y ≠ 0›
  proof -
    define Y where ‹Y = (Re y)⇧2 + (Im y)⇧2›
    moreover have ‹Y > 0›
      using that by (simp add: gauss_eq_0 less_le Y_def)
    have *: ‹Im y * (Im y * Re x) + Re x * (Re y * Re y) = Re x * Y›
      ‹Im x * (Im y * Im y) + Im x * (Re y * Re y) = Im x * Y›
      ‹(Im y)⇧2 + (Re y)⇧2 = Y›
      by (simp_all add: power2_eq_square algebra_simps Y_def)
    from ‹Y > 0› show ?thesis
      by (simp add: gauss_eq_iff algebra_simps) (simp add: * nonzero_mult_cdiv_cancel_right)
  qed
  show ‹x div y * y + x mod y = x›
    by (simp add: gauss_eq_iff)
qed

end

instantiation gauss :: euclidean_ring
begin

definition euclidean_size_gauss :: ‹gauss ⇒ nat›
  where ‹euclidean_size x = nat ((Re x)⇧2 + (Im x)⇧2)›

instance proof
  show ‹euclidean_size (0::gauss) = 0›
    by (simp add: euclidean_size_gauss_def)
  show ‹euclidean_size (x mod y) < euclidean_size y› if ‹y ≠ 0› for x y :: gauss
  proof-
    define X and Y and R and I
      where ‹X = (Re x)⇧2 + (Im x)⇧2› and ‹Y = (Re y)⇧2 + (Im y)⇧2›
        and ‹R = Re x * Re y + Im x * Im y› and ‹I = Im x * Re y - Re x * Im y›
    with that have ‹0 < Y› and rhs: ‹int (euclidean_size y) = Y›
      by (simp_all add: gauss_neq_0 euclidean_size_gauss_def)
    have ‹X * Y = R⇧2 + I⇧2›
      by (simp add: R_def I_def X_def Y_def power2_eq_square algebra_simps)
    let ?lhs = ‹X - I * (I cdiv Y) - R * (R cdiv Y)
        - I cdiv Y * (I cmod Y) - R cdiv Y * (R cmod Y)›
    have ‹?lhs = X + Y * (R cdiv Y) * (R cdiv Y) + Y * (I cdiv Y) * (I cdiv Y)
        - 2 * (R cdiv Y * R + I cdiv Y * I)›
      by (simp flip: minus_cmod_eq_mult_cdiv add: algebra_simps)
    also have ‹… = (Re (x mod y))⇧2 + (Im (x mod y))⇧2›
      by (simp add: X_def Y_def R_def I_def algebra_simps power2_eq_square)
    finally have lhs: ‹int (euclidean_size (x mod y)) = ?lhs›
      by (simp add: euclidean_size_gauss_def)
    have ‹?lhs * Y = (I cmod Y)⇧2 + (R cmod Y)⇧2›
      apply (simp add: algebra_simps power2_eq_square ‹X * Y = R⇧2 + I⇧2›)
      apply (simp flip: mult.assoc add.assoc minus_cmod_eq_mult_cdiv)
      apply (simp add: algebra_simps)
      done
    also have ‹… ≤ (Y div 2)⇧2 + (Y div 2)⇧2›
      by (rule add_mono) (use ‹Y > 0› abs_cmod_less_equal [of Y] in ‹simp_all add: power2_le_iff_abs_le›)
    also have ‹… < Y⇧2›
      using ‹Y > 0› by (cases ‹Y = 1›) (simp_all add: power2_eq_square mult_le_less_imp_less flip: mult.assoc)
    finally have ‹?lhs * Y < Y⇧2› .
    with ‹Y > 0› have ‹?lhs < Y›
      by (simp add: power2_eq_square)
    then have ‹int (euclidean_size (x mod y)) < int (euclidean_size y)›
      by (simp only: lhs rhs)
    then show ?thesis
      by simp
  qed
  show ‹euclidean_size x ≤ euclidean_size (x * y)› if ‹y ≠ 0› for x y :: gauss
  proof -
    from that have ‹euclidean_size y > 0›
      by (simp add: euclidean_size_gauss_def gauss_neq_0)
    then have ‹euclidean_size x ≤ euclidean_size x * euclidean_size y›
      by simp
    also have ‹… = nat (((Re x)⇧2 + (Im x)⇧2) * ((Re y)⇧2 + (Im y)⇧2))›
      by (simp add: euclidean_size_gauss_def nat_mult_distrib)
    also have ‹… = euclidean_size (x * y)›
      by (simp add: euclidean_size_gauss_def eq_nat_nat_iff) (simp add: algebra_simps power2_eq_square)
    finally show ?thesis .
  qed
qed

end

end