Theory SG_Generalization

(* Author: Benoît Ballenghien, Université Paris-Saclay, CNRS, ENS Paris-Saclay, LMF *)

(*<*)
theory SG_Generalization
  imports SG_Theorem
begin 
  (*>*)


section ‹Sophie Germain's Theorem: generalized Version›

text ‹The proof we give here is from \<^cite>‹kiefer2012theoreme›.›

subsection ‹Auxiliary Primes›

abbreviation non_consecutivity_condition :: ‹nat ⇒ nat ⇒ bool› (‹NC›)
  where ‹NC p q ≡ ∄x y :: int. [x ≠ 0] (mod q) ∧ [y ≠ 0] (mod q) ∧ [x ^ p = 1 + y ^ p] (mod q)›

lemma non_consecutivity_condition_bis :
  ‹NC p q ⟷ (∄x y a b. [a :: int ≠ 0] (mod q) ∧ [a ^ p = x] (mod q) ∧
                         [b :: int ≠ 0] (mod q) ∧ [b ^ p = y] (mod q) ∧ [x = 1 + y] (mod q))›
  by (simp add: cong_def) (metis mod_add_right_eq)


abbreviation not_pth_power :: ‹nat ⇒ nat ⇒ bool› (‹PPP›)
  where ‹PPP p q ≡ ∄x :: int. [p = x ^ p] (mod q)›



definition auxiliary_prime :: ‹nat ⇒ nat ⇒ bool› (‹aux'_prime›)
  where ‹aux_prime p q ≡ prime p ∧ prime q ∧ NC p q ∧ PPP p q›

lemma auxiliary_primeI :
  ‹⟦prime p; prime q; NC p q; PPP p q⟧ ⟹ aux_prime p q›
  unfolding auxiliary_prime_def by auto

lemma auxiliary_primeD :
  ‹prime p› ‹prime q› ‹NC p q› ‹PPP p q› if ‹aux_prime p q›
  using that by (auto simp add: auxiliary_prime_def)


text ‹We do not give code equation yet, let us first work around these notions.›



lemma gen_mult_group_mod_prime_as_ord : ‹ord p g = p - 1›
  if ‹prime p› ‹{1 .. p - 1} = {g ^ k mod p |k. k ∈ UNIV}›
proof -
  from that(2) have ‹g mod p ∈ {1 .. p - 1}›
    by (simp add: set_eq_iff) (metis power_one_right)
  hence ‹[g ≠ 0] (mod p)› by (simp add: cong_def)

  with cong_0_iff prime_imp_coprime ‹prime p›
  have ‹ord p g = (LEAST d. 0 < d ∧ [g ^ d = 1] (mod p))›
    unfolding ord_def by auto
  also have ‹… = p - 1›
  proof (rule ccontr)
    assume ‹(LEAST d. 0 < d ∧ [g ^ d = 1] (mod p)) ≠ p - 1›
    with fermat_theorem ‹prime p› ‹[g ≠ 0] (mod p)›
    obtain k where ‹0 < k› ‹k < p - 1› ‹[g ^ k = 1] (mod p)›
      by (metis calculation cong_0_iff coprime_ord linorder_neqE_nat
          lucas_coprime_lemma ord_minimal prime_gt_1_nat zero_less_diff)
    { fix l m
      have ‹g ^ (m + (l * k)) mod p = (g ^ m mod p * ((g ^ k) ^ l mod p)) mod p›
        by (simp add: mod_mult_eq mult.commute power_add power_mult)
      also from ‹[g ^ k = 1] (mod p)› have ‹((g ^ k) ^ l mod p) = 1›
        by (metis cong_def cong_pow mod_if power_one prime_nat_iff ‹prime p›)
      finally have ‹g ^ (m + (l * k)) mod p = g ^ m mod p› by simp
    } note $ = this

    have ‹UNIV = (⋃l. {m + (l * k) |m. m ∈ {0..k - 1}})›
      by auto (metis Suc_pred ‹0 < k› add.commute div_mod_decomp mod_Suc_le_divisor)
    hence ‹{g ^ k mod p |k. k ∈ UNIV} =
             (⋃l. {g ^ (m + (l * k)) mod p |m. m ∈ {0..k - 1}})›
      by (simp add: set_eq_iff) metis
    also have ‹… = {g ^ m mod p |m. m ∈ {0..k - 1}}› by (simp add: "$") 
    finally have ‹card {g ^ k mod p |k. k ∈ UNIV} = card …› by simp
    also have ‹{g ^ m mod p |m. m ∈ {0..k - 1}} =
                 (λm. g ^ m mod p) ` {0..k - 1}› by auto
    also from card_image_le have ‹card … ≤ card {0..k - 1}› by blast
    also have ‹… = k› by (simp add: ‹0 < k›)
    finally show False
      by (metis that(2) ‹k < p - 1› card_atLeastAtMost diff_Suc_1 linorder_not_less)
  qed
  finally show ‹ord p g = p - 1› .
qed


lemma exists_nth_power_mod_prime_iff :
  fixes p n assumes ‹prime p›
  defines d_def : ‹d ≡ gcd n (p - 1)›
  shows ‹(∃x :: int. [a = x ^ n] (mod p)) ⟷
         (n ≠ 0 ∧ [a = 0] (mod p)) ∨ [a ^ ((p - 1) div d) = 1] (mod p)›
proof (cases ‹n = 0›)
  show ‹n = 0 ⟹ ?thesis›
    by (simp add: d_def)
      (metis ‹prime p› Suc_0_not_prime_nat Suc_pred div_self power_one_right prime_gt_0_nat)
next
  show ?thesis if ‹n ≠ 0›
  proof (cases ‹[a = 0] (mod p)›)
    show ‹[a = 0] (mod p) ⟹ ?thesis›
      by (auto simp add: cong_def power_0_left ‹n ≠ 0› intro!: exI[of _ 0])
  next
    have ‹0 < int p› by (simp add: prime_gt_0_nat ‹prime p›)
    from ‹prime p› residue_prime_mult_group_has_gen gen_mult_group_mod_prime_as_ord
    obtain g where * : ‹{1 .. p - 1} = {g ^ k mod p |k. k ∈ UNIV}› and ‹ord p g = p - 1› by blast
    have ‹[g ≠ 0] (mod p)›
      by (metis ‹ord p g = p - 1› ‹prime p› nat_less_le ord_0_right_nat
          ord_cong prime_nat_iff zero_less_diff)

    show ‹?thesis› if ‹[a ≠ 0] (mod p)›
    proof (rule iffI)
      assume ‹∃x. [a = x ^ n] (mod p)›
      then obtain x where ‹[a = x ^ n] (mod p)› ..
      from cong_less_unique_int[OF ‹0 < int p›, of x]
      obtain y :: nat where ‹0 ≤ y› ‹y < p› ‹[x = y] (mod p)›
        by (metis int_nat_eq less_eq_nat.simps(1) nat_less_iff)
      from ‹[a ≠ 0] (mod p)› ‹[a = x ^ n] (mod p)› have ‹[x ≠ 0] (mod p)›
        by (metis cong_pow cong_sym cong_trans power_0_left ‹n ≠ 0›)
      with ‹[x = y] (mod p)› have ‹y ≠ 0› by (metis of_nat_0)
      with ‹0 ≤ y› ‹y < p› have ‹y ∈ {1 .. p - 1}› by simp
      with "*" ‹[x = y] (mod p)› zmod_int obtain k where ‹[x = g ^ k] (mod p)› by auto

      with ‹[a = x ^ n] (mod p)› have ‹[a = g ^ (k * n)] (mod p)›
        by (metis (no_types, lifting) cong_pow cong_trans of_nat_power power_mult)
      hence ‹[a ^ ((p - 1) div d) = (g ^ (k * n)) ^ ((p - 1) div d)] (mod p)›
        by (simp add: cong_pow)
      moreover have ‹[(g ^ (k * n)) ^ ((p - 1) div d) = g ^ (k * n * (p - 1) div d)] (mod p)›
        by (metis (no_types) d_def cong_refl div_mult_swap gcd_dvd2 power_mult)
      moreover have ‹[g ^ (k * n * (p - 1) div d) = (g ^ (k * n div d)) ^ (p - 1)] (mod p)›
        by (metis (no_types) d_def cong_def dvd_div_mult dvd_mult gcd_dvd1 power_mult)
      moreover have ‹[(g ^ (k * n div d)) ^ (p - 1) = 1] (mod p)›
        by (rule fermat_theorem[OF ‹prime p›])
          (metis ‹[g ≠ 0] (mod p)› cong_0_iff prime_dvd_power_nat ‹prime p›)
      ultimately have ‹[a ^ ((p - 1) div d) = 1] (mod p)›
        by (metis (no_types, opaque_lifting) cong_def cong_int_iff of_nat_1)
      thus ‹ n ≠ 0 ∧ [a = 0] (mod p) ∨ [a ^ ((p - 1) div d) = 1] (mod p)› ..
    next
      assume ‹ n ≠ 0 ∧ [a = 0] (mod p) ∨ [a ^ ((p - 1) div d) = 1] (mod p)›
      with ‹[a ≠ 0] (mod p)› have ‹[a ^ ((p - 1) div d) = 1] (mod p)› by blast

      from cong_less_unique_int[OF ‹0 < int p›, of a]
      obtain b :: nat where ‹0 ≤ b› ‹b < p› ‹[a = b] (mod p)›
        by (metis int_nat_eq less_eq_nat.simps(1) nat_less_iff)
      from ‹[a ≠ 0] (mod p)› ‹[a = b] (mod p)› have ‹b ≠ 0› by (metis of_nat_0)
      with ‹0 ≤ b› ‹b < p› have ‹b ∈ {1 .. p - 1}› by simp
      with "*" have ‹b ∈ {g ^ k mod p |k. k ∈ UNIV}› by blast
      with ‹[a = b] (mod p)› zmod_int obtain k where ‹[a = g ^ k] (mod p)› by auto
      from this[THEN cong_pow, of ‹(p - 1) div d›] ‹[a ^ ((p - 1) div d) = 1] (mod p)›
      have ‹[(g ^ k) ^ ((p - 1) div d) = 1] (mod p)›
        by (simp flip: cong_int_iff) (metis (no_types) cong_def)
      hence ‹[g ^ (k * (p - 1) div d) = 1] (mod p)›
        by (metis (no_types) d_def div_mult_swap gcd_dvd2 power_mult)
      hence ‹p - 1 dvd k * (p - 1) div d›
        by (simp add: ord_divides' ‹ord p g = p - 1›)
      hence ‹d dvd k›
        by (metis ‹prime p› d_def div_mult_swap dvd_div_eq_0_iff dvd_mult_div_cancel
            dvd_times_right_cancel_iff gcd_dvd2 less_numeral_extra(3) prime_gt_1_nat zero_less_diff)
      then obtain l where ‹k = l * d› by (metis dvd_div_mult_self)
      moreover from bezout_nat[OF ‹n ≠ 0›]
      obtain u v where ‹u * n = v * (p - 1) + d› by (metis d_def mult.commute)
      ultimately have ‹l * u * n = l * v * (p - 1) + k›
        by (metis distrib_left mult.assoc)
      hence ‹(g ^ (l * u)) ^ n = (g ^ (l * v)) ^ (p - 1) * g ^ k›
        by (metis power_add power_mult)
      hence ‹[(g ^ (l * u)) ^ n = (g ^ (l * v)) ^ (p - 1) * g ^ k] (mod p)› by simp
      moreover have ‹[(g ^ (l * v)) ^ (p - 1) = 1] (mod p)›
        by (metis ‹ord p g = p - 1› dvd_triv_right ord_divides power_mult)
      ultimately have ‹[(g ^ (l * u)) ^ n = g ^ k] (mod p)›
        by (metis cong_scalar_right cong_trans mult_1)
      with ‹[a = g ^ k] (mod p)› have ‹[a = (g ^ (l * u)) ^ n] (mod p)›
        by (meson cong_int_iff cong_sym cong_trans)
      thus ‹∃x. [a = x ^ n] (mod p)› by auto
    qed
  qed
qed


corollary not_pth_power_iff :
  ‹PPP p q ⟷ [p ≠ 0] (mod q) ∧ [p ^ ((q - 1) div gcd p (q - 1)) ≠ 1] (mod q)›
  if ‹prime p› ‹prime q›
  by (subst exists_nth_power_mod_prime_iff[OF ‹prime q›, of p p])
    (metis cong_int_iff not_prime_0 of_nat_0 of_nat_1 of_nat_power ‹prime p›)

corollary not_pth_power_iff_mod :
  ‹PPP p q ⟷ ¬ q dvd p ∧ p ^ ((q - 1) div gcd p (q - 1)) mod q ≠ 1›
  if ‹prime p› and ‹prime q›
  by (subst not_pth_power_iff[OF ‹prime p› ‹prime q›])
    (simp add: cong_def mod_eq_0_iff_dvd prime_gt_Suc_0_nat)


lemma non_consecutivity_condition_iff_enum_mod :
  ― ‹This version is oriented towards code generation.›
  ‹NC p q ⟷
   (∀x ∈ {1..q - 1}. let x_p_mod = x ^ p mod q
                      in ∀y ∈ {1..q - 1}. x_p_mod ≠ (1 + y ^ p mod q) mod q)›
  if ‹q ≠ 0›
proof (unfold Let_def, intro iffI conjI ballI)
  fix x y assume ‹NC p q› ‹x ∈ {1..q - 1}› ‹ y ∈ {1..q - 1}›
  show ‹x ^ p mod q ≠ (1 + y ^ p mod q) mod q›
  proof (rule ccontr)
    assume ‹¬ x ^ p mod q ≠ (1 + y ^ p mod q) mod q›
    hence ‹[x ^ p = 1 + y ^ p] (mod q)›
      by (simp add: cong_def) presburger
    with ‹NC p q› have ‹[x = 0] (mod q) ∨ [y = 0] (mod q)›
      by (metis (mono_tags, opaque_lifting) cong_int_iff int_ops(1)
          of_nat_Suc of_nat_power_eq_of_nat_cancel_iff plus_1_eq_Suc)
    with cong_less_modulus_unique_nat
    have ‹x ∉ {1..q - 1} ∨ y ∉ {1..q - 1}› by force
    with ‹x ∈ {1..q - 1}› ‹ y ∈ {1..q - 1}› show False by blast
  qed
next
  show ‹NC p q› if * : ‹∀x∈{1..q - 1}. ∀y∈{1..q - 1}. x ^ p mod q ≠ (1 + y ^ p mod q) mod q›
  proof (rule ccontr)
    assume ‹¬ NC p q›
    then obtain x y :: int where ‹[x ≠ 0] (mod q)› ‹[y ≠ 0] (mod q)›
      ‹[x ^ p = 1 + y ^ p] (mod q)› by blast

    from ‹[x ≠ 0] (mod q)› ‹q ≠ 0› have ‹x mod q ∈ {1..q - 1}›
      by (simp add: cong_0_iff)
        (metis linorder_not_le mod_by_1 mod_eq_0_iff_dvd
          mod_pos_pos_trivial of_nat_0_less_iff pos_mod_sign)
    then obtain x' :: nat where ‹x' ∈ {1..q - 1}› ‹x' = x mod q›
      by (cases ‹x mod q›) simp_all

    from ‹[y ≠ 0] (mod q)› ‹q ≠ 0› have ‹y mod q ∈ {1..q - 1}›
      by (simp add: cong_0_iff)
        (metis linorder_not_le mod_by_1 mod_eq_0_iff_dvd
          mod_pos_pos_trivial of_nat_0_less_iff pos_mod_sign)
    then obtain y' :: nat where ‹y' ∈ {1..q - 1}› ‹y' = y mod q›
      by (cases ‹y mod q›) simp_all

    from ‹[x ^ p = 1 + y ^ p] (mod q)›
    have ‹(x mod q) ^ p mod q = (1 + (y mod q) ^ p mod q) mod q›
      by (simp add: cong_def mod_add_right_eq power_mod)
    hence ‹x' ^ p mod q = (1 + y' ^ p mod q) mod q›
      by (metis ‹x' = x mod int q› ‹y' = y mod int q› nat_mod_as_int
          of_nat_Suc of_nat_power plus_1_eq_Suc zmod_int)
    with "*" ‹x' ∈ {1..q - 1}› ‹y' ∈ {1..q - 1}› show False by blast
  qed
qed


lemma auxiliary_prime_iff_enum_mod [code] :
  ― ‹We will have a more optimized version later.›
  ‹aux_prime p q ⟷
   prime p ∧ prime q ∧
   ¬ q dvd p ∧ p ^ ((q - 1) div gcd p (q - 1)) mod q ≠ 1 ∧
   (∀x ∈ {1..q - 1}. let x_p_mod = x ^ p mod q
                      in ∀y ∈ {1..q - 1}. x_p_mod ≠ (1 + y ^ p mod q) mod q)›
proof (cases ‹prime p›; cases ‹prime q›)
  assume ‹prime p› and ‹prime q›
  from ‹prime q› have ‹q ≠ 0› by auto
  show ?thesis
    unfolding auxiliary_prime_def not_pth_power_iff_mod[OF ‹prime p› ‹prime q›]
      non_consecutivity_condition_iff_enum_mod[OF ‹q ≠ 0›] by blast
next
  show ‹¬ prime q ⟹ ?thesis›
    and ‹¬ prime q ⟹ ?thesis›
    and ‹¬ prime p ⟹ ?thesis›
    by (simp_all add: auxiliary_prime_def)
qed


text ‹We can for example compute pairs of auxiliary primes less than \<^term>‹110 :: nat›.›

value ‹[(p, q). p ← [1..110], q ← [1..110], aux_prime (nat p) (nat q)]›


lemma auxiliary_primeI' :
  ‹⟦prime p; prime q; ¬ q dvd p; p ^ ((q - 1) div gcd p (q - 1)) mod q ≠ 1;
    ⋀x y. x ∈ {1..q - 1} ⟹ y ∈ {1..q - 1} ⟹ [x ^ p ≠ 1 + y ^ p] (mod q)⟧
   ⟹ aux_prime p q›
  by (simp add: auxiliary_prime_iff_enum_mod cong_def mod_Suc_eq)



lemma two_is_not_auxiliary_prime : ‹¬ aux_prime p 2›
  by (simp add: auxiliary_prime_iff_enum_mod) presburger


lemma auxiliary_prime_of_2 : ‹aux_prime 2 q ⟷ q = 3 ∨ q = 5›
proof (rule iffI)
  show ‹q = 3 ∨ q = 5 ⟹ aux_prime 2 q›
  proof (elim disjE)
    show ‹q = 3 ⟹ aux_prime 2 q›
    proof (intro auxiliary_primeI')
      show ‹prime (2 :: nat)› and ‹q = 3 ⟹ prime q› by simp_all
    next
      fix x y assume ‹q = 3› ‹x ∈ {1..q - 1}› ‹y ∈ {1..q - 1}›
      hence ‹x = 1 ∧ y = 1 ∨ x = 1 ∧ y = 2 ∨ x = 2 ∧ y = 1 ∨ x = 2 ∧ y = 2›
        by simp (metis le_Suc_eq le_antisym numeral_2_eq_2)
      with ‹q = 3› show ‹[x⇧2 ≠ 1 + y⇧2] (mod q)›
        by (fastforce simp add: cong_def)
    next
      show ‹q = 3 ⟹ ¬ q dvd 2› by simp
    next
      show ‹q = 3 ⟹ 2 ^ ((q - 1) div gcd 2 (q - 1)) mod q ≠ 1› by simp
    qed
  next
    show ‹q = 5 ⟹ aux_prime 2 q›
    proof (intro auxiliary_primeI')
      show ‹prime (2 :: nat)› and ‹q = 5 ⟹ prime q› by simp_all
    next
      fix x y assume ‹q = 5› ‹x ∈ {1..q - 1}› ‹y ∈ {1..q - 1}›
      hence ‹(x = 1 ∨ x = 2 ∨ x = 3 ∨ x = 4) ∧ (y = 1 ∨ y = 2 ∨ y = 3 ∨ y = 4)›
        by (simp add: numeral_eq_Suc) linarith
      with ‹q = 5› show ‹[x⇧2 ≠ 1 + y⇧2] (mod q)›
        by (fastforce simp add: cong_def)
    next
      show ‹q = 5 ⟹ ¬ q dvd 2› by simp
    next
      show ‹q = 5 ⟹ 2 ^ ((q - 1) div gcd 2 (q - 1)) mod q ≠ 1› by simp
    qed
  qed
next
  assume aux_q : ‹aux_prime 2 q›
  with two_is_not_auxiliary_prime have ‹q ≠ 2› by blast
  show ‹q = 3 ∨ q = 5›
  proof (rule ccontr)
    assume ‹¬ (q = 3 ∨ q = 5)›
    with ‹q ≠ 2› auxiliary_primeD(1-2)[OF aux_q] prime_prime_factor
      le_neq_implies_less prime_ge_2_nat
    have ‹prime q› ‹odd q› ‹2 < q› by metis+

    hence ‹5 ≤ q›
      by (metis Suc_1 ‹¬ (q = 3 ∨ q = 5)› add.commute add_Suc_right eval_nat_numeral(3)
          even_numeral not_less_eq_eq numeral_Bit0 order_antisym_conv plus_1_eq_Suc prime_ge_2_nat)
    with ‹prime q› have ‹gcd 4 q = (1 :: int)›
      by (metis coprime_imp_gcd_eq_1 eval_nat_numeral(3) gcd.commute less_Suc_eq
          of_nat_1 order_less_le_trans prime_nat_iff'' zero_less_numeral)
    with cong_solve_dvd_int obtain inv_4 :: int
      where inv_4: ‹[4 * inv_4 = 1] (mod q)›
      by (metis dvd_refl gcd_int_int_eq of_nat_numeral)
    define x where ‹x ≡ 1 + inv_4›
    define y where ‹y ≡ 1 - inv_4›

    from inv_4 have ‹[x⇧2 = 1 + y⇧2] (mod q)›
      by (simp add: x_def y_def power2_eq_square cong_iff_dvd_diff ring_class.ring_distribs)
    moreover obtain x' y' :: nat where ‹[x' = x] (mod q)› ‹[y' = y] (mod q)›
      by (metis ‹prime q› cong_less_unique_int cong_sym int_eq_iff of_nat_0_less_iff prime_gt_0_nat)
    ultimately have ‹[x'⇧2 = 1 + y'⇧2] (mod q)›
      by (simp flip: cong_int_iff)
        (meson cong_add cong_pow cong_refl cong_sym_eq cong_trans)
    moreover have ‹[x' ≠ 0] (mod q)›
    proof (rule ccontr)
      assume ‹¬ [x' ≠ 0] (mod q)›
      with ‹[x' = x] (mod q)› have ‹[x = 0] (mod q)›
        by (metis cong_0_iff cong_dvd_iff int_dvd_int_iff)
      hence ‹[4 * x = 0] (mod q)›
        by (metis cong_scalar_left mult_zero_right)
      with cong_add[OF cong_refl[of 4] inv_4] have ‹q dvd 5›
        by (simp add: x_def) (metis cong_0_iff cong_dvd_iff int_ops(3) of_nat_dvd_iff)
      with ‹prime q› have ‹q = 5› by (auto intro: primes_dvd_imp_eq)
      with ‹¬ (q = 3 ∨ q = 5)› show False by blast
    qed
    moreover have ‹[y' ≠ 0] (mod q)›
    proof (rule ccontr)
      assume ‹¬ [y' ≠ 0] (mod q)›
      with ‹[y' = y] (mod q)› have ‹[y = 0] (mod q)›
        by (metis cong_0_iff cong_dvd_iff int_dvd_int_iff)
      hence ‹[4 * y = 0] (mod q)›
        by (metis cong_scalar_left mult_zero_right)
      with cong_diff[OF cong_refl[of 4] inv_4] have ‹q dvd 3›
        by (simp add: y_def) (metis cong_0_iff cong_dvd_iff int_ops(3) of_nat_dvd_iff)
      with ‹prime q› have ‹q = 3› by (auto intro: primes_dvd_imp_eq)
      with ‹¬ (q = 3 ∨ q = 5)› show False by blast
    qed
    ultimately have ‹[(int x')⇧2 = 1 + (int y')⇧2] (mod q) ∧
      [int x' ≠ 0] (mod q) ∧ [int y' ≠ 0] (mod q)›
      by (metis cong_int_iff of_nat_0 of_nat_Suc of_nat_power plus_1_eq_Suc)
    with auxiliary_primeD(3) aux_q show False by blast
  qed
qed




text ‹An auxiliary prime ‹q› of ‹p› is generally of the form \<^term>‹q = 2 * n * p + 1›.›

lemma auxiliary_prime_pattern_aux :
  ‹∃x y. [x ≠ 0] (mod q) ∧ [y ≠ 0] (mod q) ∧ [x ^ p = 1 + y ^ p] (mod q)›
  if ‹p ≠ 0› ‹prime q› ‹coprime p (q - 1)› ‹odd q›
proof -
  from bezout_nat ‹coprime p (q - 1)› ‹p ≠ 0›
  obtain u v where bez : ‹u * p = 1 + v * (q - 1)›
    by (metis add.commute coprime_imp_gcd_eq_1 mult.commute)
  have ‹[x ≠ 0] (mod q) ⟹ [(x ^ v) ^ (q - 1) = 1] (mod q)› for x
    by (meson cong_0_iff fermat_theorem prime_dvd_power ‹prime q›)
  hence * : ‹[(x ^ u) ^ p = x] (mod q)› for x
    by (fold power_mult, unfold bez, unfold power_add power_mult)
      (metis cong_0_iff cong_def cong_scalar_left ‹prime q›
        mult.right_neutral power_one_right prime_dvd_mult_iff)
  obtain x0 y0 where ‹[x0 ≠ 0] (mod q)› ‹[y0 ≠ 0] (mod q)› ‹[x0 = 1 + y0] (mod q)›
    by (metis ‹odd q› ‹prime q› cong_0_iff cong_refl not_prime_unit
        one_add_one prime_prime_factor two_is_prime_nat)
  from "*" this(3) have ‹[(x0 ^ u) ^ p = 1 + (y0 ^ u) ^ p] (mod q)›
    by (metis cong_add_lcancel_nat cong_def)
  moreover from ‹[x0 ≠ 0] (mod q)› ‹[y0 ≠ 0] (mod q)›
  have ‹[x0 ^ u ≠ 0] (mod q)› ‹[y0 ^ u ≠ 0] (mod q)›
    by (meson cong_0_iff prime_dvd_power ‹prime q›)+
  ultimately show ?thesis by blast
qed


theorem auxiliary_prime_pattern :
  ‹p = 2 ∧ (q = 3 ∨ q = 5) ∨ odd p ∧ (∃n≥1. q = 2 * n * p + 1)› if aux_p : ‹aux_prime p q›
proof -
  from auxiliary_prime_of_2 consider ‹p = 2› ‹q = 3 ∨ q = 5› | ‹odd p› ‹q ≠ 2›
    by (metis aux_p auxiliary_prime_def prime_prime_factor two_is_not_auxiliary_prime)
  thus ?thesis
  proof cases
    show ‹p = 2 ⟹ q = 3 ∨ q = 5 ⟹ ?thesis› by blast
  next
    assume ‹odd p› ‹q ≠ 2›
    have ‹2 < q› ‹odd q›
      by (use ‹q ≠ 2› auxiliary_prime_def le_neq_implies_less prime_ge_2_nat that in presburger)
        (metis ‹q ≠ 2› auxiliary_prime_def prime_prime_factor that two_is_prime_nat)
    from aux_p have ‹prime p› and ‹prime q› by (simp_all add: auxiliary_primeD(1-2))
    from euler_criterion[OF ‹prime q› ‹2 < q›]
    have * : ‹[Legendre p q = p ^ ((q - 1) div 2)] (mod q)› by simp
    have ‹¬ coprime p (q - 1)›
      by (metis auxiliary_prime_def cong_0_iff coprime_iff_gcd_eq_1
          div_by_1 fermat_theorem not_pth_power_iff aux_p)
    with ‹prime p› prime_imp_coprime have ‹p dvd q - 1› by blast
    then obtain k where ‹q = k * p + 1›
      by (metis add.commute ‹prime q› dvd_div_mult_self
          le_add_diff_inverse less_or_eq_imp_le prime_gt_1_nat)
    with ‹odd q› ‹odd p› have ‹even k› by simp
    then obtain n where ‹k = 2 * n› by fast
    with ‹q = k * p + 1› have ‹q = 2 * n * p + 1› by blast
    with ‹2 < q› have ‹1 ≤ n› 
      by (metis One_nat_def add.commute less_one linorder_not_less
          mult_is_0 one_le_numeral plus_1_eq_Suc)
    with ‹odd p› ‹q = 2 * n * p + 1› show ?thesis by blast
  qed
qed



lemma auxiliary_prime_imp_less : ‹aux_prime p q ⟹ p < q›
  by (auto dest: auxiliary_prime_pattern simp add: less_Suc_eq)

lemma auxiliary_primeE :
  assumes ‹aux_prime p q›
  obtains ‹p = 2› ‹q = 3› | ‹p = 2› ‹q = 5›
  | n where ‹odd p› ‹1 ≤ n› ‹q = 2 * n * p + 1›
    ‹NC p (2 * n * p + 1)› ‹PPP p (2 * n * p + 1)›
  by (metis assms auxiliary_prime_def auxiliary_prime_pattern)



text ‹With this, we can quickly eliminate numbers that cannot be auxiliary.›

declare auxiliary_prime_iff_enum_mod [code del]

lemma auxiliary_prime_iff_enum_mod_optimized [code] :
  ‹aux_prime p q ⟷
   p = 2 ∧ (q = 3 ∨ q = 5) ∨
   p < q ∧
   2 * p dvd q - 1 ∧
   prime p ∧ prime q ∧
   ¬ q dvd p ∧ p ^ ((q - 1) div gcd p (q - 1)) mod q ≠ 1 ∧
   (∀x ∈ {1..q - 1}. let x_p_mod = x ^ p mod q
                      in ∀y ∈ {1..q - 1}. x_p_mod ≠ (1 + y ^ p mod q) mod q)›
  by (fold auxiliary_prime_iff_enum_mod)
    (metis add_diff_cancel_right' auxiliary_prime_imp_less auxiliary_prime_of_2
      auxiliary_prime_pattern dvd_refl even_mult_iff mult_dvd_mono)

value ‹[(p, q). p ← [1..1000], q ← [1..110], aux_prime (nat p) (nat q)]›



subsection ‹Sophie Germain Primes are auxiliary›

text ‹When ‹p› is an \<^const>‹odd› \<^const>‹prime› and \<^term>‹2 * p + 1 :: nat› is also
      a \<^const>‹prime› (what we call a Sophie Germain \<^const>‹prime›),
      \<^term>‹2 * p + 1 :: nat› is automatically an \<^const>‹auxiliary_prime›.›

lemma SophGer_prime_imp_auxiliary_prime :
  fixes p assumes ‹SG p› defines q_def: ‹q ≡ 2 * p + 1›
  shows ‹aux_prime p q›
proof (rule auxiliary_primeI)
  from SophGer_primeD(2-3)[OF ‹SG p›]
  show ‹prime p› and ‹prime q› by (unfold q_def)
next
  from SophGer_primeD[OF ‹SG p›, folded q_def]
  have ‹odd p› ‹prime p› ‹prime q› ‹prime (int q)› ‹p ≠ 0› by simp_all
  show ‹NC p q›
  proof (rule ccontr)
    assume ‹¬ NC p q›
    then obtain x y :: int where ‹[x ≠ 0] (mod q)› ‹[y ≠ 0] (mod q)›
      ‹[x ^ p = 1 + y ^ p] (mod q)› by blast
    from SG_simps.p_th_power_mod_q ‹[x ≠ 0] (mod q)› ‹SG p› cong_0_iff q_def
    consider ‹[x ^ p = 1] (mod q)› | ‹[x ^ p = - 1] (mod q)› by blast
    thus False
    proof cases
      assume ‹[x ^ p = 1] (mod q)›
      with ‹[x ^ p = 1 + y ^ p] (mod q)› have ‹[y ^ p = 0] (mod q)›
        by (metis add.right_neutral cong_add_lcancel cong_sym cong_trans)
      with ‹[y ≠ 0] (mod q)› show False
        by (meson ‹prime (int q)› cong_0_iff prime_dvd_power_int)
    next
      assume ‹[x ^ p = - 1] (mod q)›
      with ‹[x ^ p = 1 + y ^ p] (mod q)›
      have ‹[- 1 = 1 + y ^ p] (mod q)› by (metis cong_def)
      moreover have ‹- (1::int) = 1 + - 2› by force
      ultimately have ‹[y ^ p = - 2] (mod q)›
        by (metis cong_add_lcancel cong_sym)
      with ‹odd p› cong_minus_minus_iff have ‹[(- y) ^ p = 2] (mod q)› by force
      with cong_sym have ‹∃x :: int. [2 = x ^ p] (mod q)› by blast
      with ‹p ≠ 0› exists_nth_power_mod_prime_iff[OF ‹prime q›]
      have ‹([2 = 0] (mod q) ∨ [4 = 1] (mod q))›
        by (simp add: q_def flip: cong_int_iff)
      hence ‹p ≤ 1›
      proof (elim disjE)
        from ‹p ≠ 0› show ‹[2 = 0] (mod q) ⟹ p ≤ 1›
          by (auto simp add: q_def cong_def)
      next
        from linorder_not_less show ‹[4 = 1] (mod q) ⟹ p ≤ 1›
          by (fastforce simp add: q_def cong_def)
      qed
      with ‹SG p› show False
        by (metis ‹prime p› linorder_not_less prime_nat_iff)
    qed
  qed
next
  from ‹SG p›[THEN SophGer_primeD(3), folded q_def]
  have ‹prime q› ‹prime (int q)› by simp_all
  from SG_simps.p_th_power_mod_q[OF ‹SG p›]
  have ‹¬ q dvd x ⟹ [x ^ p = 1] (mod q) ∨ [x ^ p = - 1] (mod q)› for x :: int
    unfolding q_def .
  moreover have ‹[p ≠ (1 :: int)] (mod q)› ‹[p ≠ - 1] (mod q)›
    using SG_simps.notcong_p(1, 3)[OF ‹SG p›] cong_sym unfolding q_def by blast+
  ultimately have ‹¬ q dvd x ⟹ [p ≠ x ^ p] (mod q)› for x :: int
    using cong_trans by blast
  moreover have ‹q dvd x ⟹ [p ≠ x ^ p] (mod q)› for x :: int
    by (metis SG_simps.pos Suc_eq_plus1 ‹SG p› cong_dvd_iff dvd_power dvd_trans
        int_dvd_int_iff less_add_Suc1 mult.commute mult_2_right nat_dvd_not_less q_def)
  ultimately show ‹PPP p q› by blast
qed



subsection ‹Main Theorems›

theorem Sophie_Germain_auxiliary_prime :
  ‹q dvd x ∨ q dvd y ∨ q dvd z›
  if ‹x ^ p + y ^ p = z ^ p› and ‹aux_prime p q› for x y z :: int
proof (rule ccontr)
  assume not_dvd : ‹¬ (q dvd x ∨ q dvd y ∨ q dvd z)›
  from auxiliary_primeD[OF ‹aux_prime p q›]
  have ‹prime p› ‹prime q› ‹NC p q› by simp_all

  have ‹coprime x q›
    by (meson coprime_commute not_dvd prime_imp_coprime prime_nat_int_transfer ‹prime q›)
  with bezout_int[of x q] obtain inv_x v :: int where ‹inv_x * x + v * q = 1› by auto
  hence inv_x : ‹[x * inv_x = 1] (mod q)› by (metis cong_iff_lin mult.commute)

  from ‹x ^ p + y ^ p = z ^ p› have ‹z ^ p = x ^ p + y ^ p› ..
  hence ‹(inv_x * z) ^ p = (inv_x * x) ^ p + (inv_x * y) ^ p›
    by (metis distrib_left power_mult_distrib)
  with inv_x have ‹[(inv_x * z) ^ p = 1 + (inv_x * y) ^ p] (mod q)›
    by (metis cong_add_rcancel cong_pow mult.commute power_one)
  moreover from inv_x ‹prime q› have ‹[(inv_x * z) ^ p ≠ 0] (mod q)›
    by (metis cong_dvd_iff dvd_0_right not_dvd not_prime_unit
        prime_dvd_mult_eq_int prime_dvd_power prime_nat_int_transfer)
  moreover from inv_x ‹prime q› have ‹[(inv_x * y) ^ p ≠ 0] (mod q)›
    by (metis cong_dvd_iff dvd_0_right not_dvd not_prime_unit
        prime_dvd_mult_eq_int prime_dvd_power prime_nat_int_transfer)
  moreover obtain y' z' :: int where ‹[y' = inv_x * y] (mod q)› ‹[z' = inv_x * z] (mod q)›
    by (metis ‹prime q› cong_less_unique_int cong_sym int_eq_iff of_nat_0_less_iff prime_gt_0_nat)
  ultimately show False
    by (metis ‹NC p q› ‹prime p› ‹prime q› cong_0_iff
        prime_dvd_power_iff prime_gt_0_nat prime_nat_int_transfer)
qed



theorem Sophie_Germain_generalization :
  ‹∄x y z :: int. x ^ p + y ^ p = z ^ p ∧
                  [x ≠ 0] (mod p⇧2) ∧ [y ≠ 0] (mod p⇧2) ∧ [z ≠ 0] (mod p⇧2)›
  if odd_p : ‹odd p› and aux_prime : ‹aux_prime p q›
proof (rule ccontr) ― ‹The proof is done by contradiction.›
  from ‹aux_prime p q› have prime_p : ‹prime p›
    by (metis auxiliary_primeD(1))
  hence not_p_0 : ‹p ≠ 0› and prime_int_p : ‹prime (int p)› by simp_all
  from ‹aux_prime p q› have prime_q : ‹prime q›
    by (metis auxiliary_primeD(2))
  hence prime_int_q : ‹prime (int q)› by simp
  from ‹odd p› ‹prime p› have p_ge_3 : ‹3 ≤ p›
    by (simp add: numeral_eq_Suc)
      (metis Suc_le_eq dvd_refl le_antisym not_less_eq_eq prime_gt_Suc_0_nat)

  assume ‹¬ (∄x y z. x ^ p + y ^ p = z ^ p ∧ [x ≠ 0] (mod int (p⇧2)) ∧
                     [y ≠ 0] (mod int (p⇧2)) ∧ [z ≠ 0] (mod (int p)⇧2))›
  then obtain x y z :: int
    where fermat : ‹x ^ p + y ^ p = z ^ p›
      and not_cong_0 : ‹[x ≠ 0] (mod p⇧2)› ‹[y ≠ 0] (mod p⇧2)› ‹[z ≠ 0] (mod p⇧2)› by auto

― ‹We first assume without loss of generality that
    \<^term>‹x›, \<^term>‹y› and \<^term>‹z› are setwise \<^const>‹coprime›.›
  let ?gcd = ‹Gcd {x, y, z}›
  wlog coprime : ‹?gcd = 1› goal False generalizing x y z keeping fermat not_cong_0
    using FLT_setwise_coprime_reduction_mod_version[OF fermat not_cong_0]
      hypothesis by blast

― ‹Then we can deduce that \<^term>‹x›, \<^term>‹y› and \<^term>‹z› are pairwise \<^const>‹coprime›.›
  have coprime_x_y : ‹coprime x y›
    by (fact FLT_setwise_coprime_imp_pairwise_coprime[OF ‹p ≠ 0› fermat coprime])
  have coprime_y_z : ‹coprime y z›
  proof (subst coprime_minus_right_iff[symmetric],
      rule FLT_setwise_coprime_imp_pairwise_coprime[OF ‹p ≠ 0›])
    from fermat ‹odd p› show ‹y ^ p + (- z) ^ p = (- x) ^ p› by simp
  next
    show ‹Gcd {y, - z, - x} = 1›
      by (metis Gcd_insert coprime gcd_neg1_int insert_commute)
  qed
  have coprime_x_z : ‹coprime x z›
  proof (subst coprime_minus_right_iff[symmetric],
      rule FLT_setwise_coprime_imp_pairwise_coprime[OF ‹p ≠ 0›])
    from fermat ‹odd p› show ‹x ^ p + (- z) ^ p = (- y) ^ p› by simp
  next
    show ‹Gcd {x, - z, - y} = 1› 
      by (metis Gcd_insert coprime gcd_neg1_int insert_commute)
  qed

― ‹From @{thm [show_question_marks = false] Sophie_Germain_auxiliary_prime}
    we have that among \<^term>‹x›, \<^term>‹y› and \<^term>‹z›,
    one (and only one, see below) is a multiple of \<^term>‹q›.›
  from Sophie_Germain_auxiliary_prime[OF fermat aux_prime]
  have q_dvd_xyz : ‹q dvd x ∨ q dvd y ∨ q dvd z› .

― ‹Without loss of generality, we can assume that \<^term>‹x› is a multiple of \<^term>‹q›.›
  wlog q_dvd_z : ‹q dvd z› goal False generalizing x y z
    keeping fermat not_cong_0 coprime_x_y coprime_y_z coprime_x_z
  proof -
    from negation q_dvd_xyz have ‹q dvd x ∨ q dvd y› by simp
    thus False
    proof (elim disjE)
      show ‹q dvd x ⟹ False›
        by (erule hypothesis[of x ‹- y› z])
          (use fermat not_cong_0 ‹odd p› in
            ‹simp_all add: cong_0_iff coprime_commute coprime_x_y coprime_x_z coprime_y_z›)
    next
      show ‹q dvd y ⟹ False›
        by (erule hypothesis[of y ‹- x› z])
          (use fermat not_cong_0 ‹odd p› in
            ‹simp_all add: cong_0_iff coprime_commute coprime_x_y coprime_x_z coprime_y_z›)
    qed
  qed

  define S where ‹S ≡ λy z :: int. ∑k = 0..p - 1. (- z) ^ (p - 1 - k) * y ^ k›

― ‹Now we prove that \<^term>‹x›, \<^term>‹y› or \<^term>‹x› is dividable by \<^term>‹p›.›
  have p_dvd_xyz : ‹p dvd x ∨ p dvd y ∨ p dvd z›
  proof (rule ccontr)
    assume ‹¬ (p dvd x ∨ p dvd y ∨ p dvd z)›
    hence ‹[x ≠ 0] (mod p)› ‹[y ≠ 0] (mod p)› ‹[z ≠ 0] (mod p)› 
      by (simp_all add: cong_0_iff)
    from Sophie_Germain_lemma[OF ‹odd p› ‹prime p›, of ‹- z› x y]
      coprime_x_y fermat ‹[z ≠ 0] (mod p)›
    obtain a α where ‹x + y = a ^ p› ‹S x y = α ^ p›
      by (simp add: S_def ‹odd p› coprime_commute) (meson cong_0_iff dvd_minus_iff)
    from Sophie_Germain_lemma[OF ‹odd p› ‹prime p›, of ‹- x› z ‹- y›]
      coprime_y_z fermat ‹[x ≠ 0] (mod int p)›
    obtain b β where ‹z - y = b ^ p› ‹S z (- y) = β ^ p›
      by (simp add: S_def ‹odd p› coprime_commute) (meson cong_0_iff dvd_minus_iff)
    from Sophie_Germain_lemma[OF ‹odd p› ‹prime p›, of ‹- y› z ‹- x›]
      coprime_x_z fermat ‹[y ≠ 0] (mod p)›
    obtain c γ where ‹z - x = c ^ p› ‹S z (- x) = γ ^ p›
      by (simp add: S_def ‹odd p› coprime_commute) (meson cong_0_iff dvd_minus_iff)

    have ‹a ^ p + b ^ p + c ^ p = x + y + (z - y) + (z - x)›
      by (simp add: ‹x + y = a ^ p› ‹z - y = b ^ p› ‹z - x = c ^ p›)
    also have ‹… = 2 * z› by simp
    also from ‹q dvd z› have ‹[… = 0] (mod q)› by (simp add: cong_0_iff)
    finally have ‹[a ^ p + b ^ p + c ^ p = 0] (mod q)› .

    have ‹[a = 0] (mod q) ∨ [b = 0] (mod q) ∨ [c = 0] (mod q)›
    proof (rule ccontr)
      assume ‹¬ ([a = 0] (mod q) ∨ [b = 0] (mod q) ∨ [c = 0] (mod q))›
      hence ‹[a ≠ 0] (mod q)› ‹[b ≠ 0] (mod q)› ‹[c ≠ 0] (mod q)› by simp_all
      from ‹[c ≠ 0] (mod q)› have ‹gcd c q = 1›
        by (meson aux_prime auxiliary_prime_def cong_0_iff coprime_iff_gcd_eq_1
            residues_prime.p_coprime_right_int residues_prime_def)
      from bezout_int[of c q, unfolded this]
      obtain u v where ‹u * c + v * int q = 1› by blast
      with ‹[a ≠ 0] (mod q)› have ‹[u ≠ 0] (mod q)›
        by (metis cong_0_iff cong_dvd_iff cong_iff_lin dvd_mult mult.commute unit_imp_dvd)

      from ‹[a ^ p + b ^ p + c ^ p = 0] (mod q)›
      have ‹[(u * a) ^ p + (u * b) ^ p + (u * c) ^ p = 0] (mod q)›
        by (simp add: power_mult_distrib)
          (metis cong_scalar_left distrib_left mult.commute mult_zero_left)
      also from ‹u * c + v * int q = 1› have ‹u * c = 1 - v * int q› by simp
      finally have ‹[(u * a) ^ p + (u * b) ^ p + (1 - v * int q) ^ p = 0] (mod q)› .
      moreover have ‹[(1 - v * int q) ^ p = 1] (mod q)›
        by (metis add_uminus_conv_diff cong_0_iff cong_add_lcancel_0
            cong_pow dvd_minus_iff dvd_triv_right power_one)
      ultimately have ‹[(u * a) ^ p + (u * b) ^ p + 1 = 0] (mod q)›
        by (meson cong_add_lcancel cong_sym cong_trans)
      hence ‹[1 + (u * b) ^ p = (- (u * a)) ^ p] (mod q)›
        by (simp add: ‹odd p› cong_iff_dvd_diff) presburger
      hence ‹[(- (u * a)) ^ p = 1 + (u * b) ^ p] (mod q)› by (fact cong_sym)
      moreover from ‹[a ≠ 0] (mod q)› ‹[u ≠ 0] (mod q)›
        cong_prime_prod_zero_int[OF _ ‹prime (int q)›, of u a] cong_minus_minus_iff
      have ‹[- u * a ≠ 0] (mod q)› by force
      moreover from ‹[b ≠ 0] (mod q)› ‹[u ≠ 0] (mod q)›
        cong_prime_prod_zero_int[OF _ ‹prime (int q)›, of u b]
      have ‹[u * b ≠ 0] (mod q)› by blast
      ultimately show False
        using aux_prime auxiliary_primeD(3) by auto
    qed
    hence ‹q dvd a›
    proof (elim disjE)
      show ‹[a = 0] (mod q) ⟹ q dvd a› by (simp add: cong_0_iff)
    next
      assume ‹[b = 0] (mod q)›
      with ‹z - y = b ^ p› ‹q dvd z› ‹prime p› have ‹q dvd y›
        by (metis cong_0_iff cong_dvd_iff cong_iff_dvd_diff
            dvd_power dvd_trans prime_gt_0_nat)
      with ‹prime (int q)› ‹q dvd z› ‹coprime y z› have False
        by (metis coprime_def not_prime_unit)
      thus ‹q dvd a› ..
    next
      assume ‹[c = 0] (mod q)›
      with ‹z - x = c ^ p› ‹q dvd z› ‹prime p› have ‹q dvd x›
        by (metis cong_0_iff cong_dvd_iff cong_iff_dvd_diff
            dvd_power dvd_trans prime_gt_0_nat)
      with ‹prime (int q)› ‹q dvd z› ‹coprime x z› have False
        by (metis coprime_def not_prime_unit)
      thus ‹q dvd a› ..
    qed
    with ‹x + y = a ^ p› ‹p ≠ 0› ‹prime (int q)› have ‹[y = - x] (mod q)›
      by (metis add.commute add.inverse_inverse add_uminus_conv_diff
          bot_nat_0.not_eq_extremum cong_iff_dvd_diff prime_dvd_power_int_iff)
    hence ‹[S x y = p * x ^ (p - 1)] (mod q)›
      unfolding S_def by (fact Sophie_Germain_lemma_computation_cong_simp[OF ‹p ≠ 0›])
    moreover from ‹z - x = c ^ p› ‹q dvd z› have ‹[x = (- c) ^ p] (mod q)›
      by (simp add: ‹odd p› cong_iff_dvd_diff)
        (metis add_diff_cancel_left' diff_diff_eq2)
    ultimately have ‹[S x y = p * ((- c) ^ p) ^ (p - 1)] (mod q)›
      by (meson cong_pow cong_scalar_left cong_trans)
    also have ‹S x y = α ^ p› by (fact ‹S x y = α ^ p›)
    also have ‹((- c) ^ p) ^ (p - 1) = ((- c) ^ (p - 1)) ^ p›
      by (metis mult.commute power_mult)
    finally have ‹[α ^ p = p * ((- c) ^ (p - 1)) ^ p] (mod q)› .

    have ‹gcd q ((- c) ^ (p - 1)) = 1›
      by (metis ‹[x = (- c) ^ p] (mod int q)› ‹int q dvd z› cong_dvd_iff
          coprime_def coprime_imp_gcd_eq_1 coprime_x_z dvd_mult not_prime_unit
          power_eq_if prime_imp_coprime_int ‹p ≠ 0› ‹prime (int q)›)
    with bezout_int obtain u v
      where ‹u * int q + v * (- c) ^ (p - 1) = 1› by metis
    hence ‹v * (- c) ^ (p - 1) = 1 - u * int q› by simp
    have ‹[1 - u * int q = 1] (mod q)› by (simp add: cong_iff_lin)
    from ‹[α ^ p = p * ((- c) ^ (p - 1)) ^ p] (mod q)›
    have ‹[(v * α) ^ p = p * (v * (- c) ^ (p - 1)) ^ p] (mod q)›
      by (simp add: power_mult_distrib) (metis cong_scalar_left mult.left_commute)
    with ‹[1 - u * int q = 1] (mod int q)› have ‹[(v * α) ^ p = p] (mod q)›
      by (unfold ‹v * (- c) ^ (p - 1) = 1 - u * int q›)
        (metis cong_pow cong_scalar_left cong_trans mult.comm_neutral power_one)
    with ‹aux_prime p q›[THEN auxiliary_primeD(4)] cong_sym show False by blast
  qed

― ‹Without loss of generality, we can assume that it is \<^term>‹z›.›
  wlog p_dvd_z : ‹p dvd z› goal False generalizing x y z S
    keeping fermat not_cong_0 coprime_x_y coprime_y_z coprime_x_z S_def
  proof -
    from negation p_dvd_xyz have ‹p dvd x ∨ p dvd y› by simp
    thus False
    proof (elim disjE)
      show ‹p dvd x ⟹ False›
        by (erule hypothesis[of x ‹- y› z])
          (use fermat not_cong_0 ‹odd p› in
            ‹simp_all add: cong_0_iff coprime_commute coprime_x_y coprime_x_z coprime_y_z›)
    next
      show ‹p dvd y ⟹ False›
        by (erule hypothesis[of y ‹- x› z])
          (use fermat not_cong_0 ‹odd p› in
            ‹simp_all add: cong_0_iff coprime_commute coprime_x_y coprime_x_z coprime_y_z›)
    qed
  qed

― ‹The rest of the proof consists in deducing that actually \<^term>‹p⇧2›
    divides \<^term>‹z›, which contradicts @{thm ‹[z ≠ 0] (mod p⇧2)›}.›
  from ‹p dvd z› coprime_x_z coprime_y_z
  have ‹[x ≠ 0] (mod p)› ‹[y ≠ 0] (mod p)›
    by (simp_all add: cong_0_iff)
      (meson not_coprimeI not_prime_unit ‹prime (int p)›)+

  from Sophie_Germain_lemma[OF ‹odd p› ‹prime p›, of ‹- x› z ‹- y›]
    coprime_y_z fermat ‹[x ≠ 0] (mod int p)›
  obtain b β where ‹z - y = b ^ p› ‹S z (- y) = β ^ p›
    by (simp add: S_def ‹odd p› coprime_commute) (meson cong_0_iff dvd_minus_iff)
  from Sophie_Germain_lemma[OF ‹odd p› ‹prime p›, of ‹- y› z ‹- x›]
    coprime_x_z fermat ‹[y ≠ 0] (mod p)›
  obtain c γ where ‹z - x = c ^ p› ‹S z (- x) = γ ^ p›
    by (simp add: S_def ‹odd p› coprime_commute) (meson cong_0_iff dvd_minus_iff)

  from fermat have ‹(- z) ^ p + x ^ p + y ^ p = 0› by (simp add: ‹odd p›)
  from Sophie_Germain_lemma_computation[OF ‹odd p›] fermat
  have ‹(x + y) * S x y = z ^ p› by (simp add: S_def)
  with ‹[z ≠ 0] (mod p⇧2)› have ‹x + y ≠ 0› ‹S x y ≠ 0› by auto

  define z' where ‹z' ≡ z div p›
  from ‹p dvd z› ‹[z ≠ 0] (mod p⇧2)› ‹p ≠ 0›
  have ‹z = z' * p› ‹[z' ≠ 0] (mod p)›
    by (simp_all add: cong_0_iff z'_def dvd_div_iff_mult power2_eq_square)

  from Sophie_Germain_lemma_only_possible_prime_common_divisor[OF ‹prime p› _ ‹coprime x y›]
  have ‹∃q∈#prime_factorization r. q ≠ p ⟹ ¬ r dvd x + y ∨ ¬ r dvd S x y› for r :: nat
    unfolding S_def
    by (meson dvd_trans in_prime_factors_iff int_dvd_int_iff
        of_nat_eq_iff prime_nat_int_transfer)
  from this[OF Ex_other_prime_factor[OF _ _ ‹prime p›]]
  have ‹r dvd x + y ⟹ r dvd S x y ⟹ r = 0 ∨ (∃k. r = p ^ k)› for r :: nat by auto
  moreover have ‹¬ (p ^ k dvd x + y ∧ p ^ k dvd S x y)› if ‹1 < k› for k
  proof (rule ccontr)
    assume ‹¬ (¬ (p ^ k dvd x + y ∧ p ^ k dvd S x y))›
    moreover from ‹1 < k› have ‹p⇧2 dvd p ^ k›
      by (simp add: le_imp_power_dvd)
    ultimately have ‹p⇧2 dvd x + y› ‹p⇧2 dvd S x y›
      by (meson dvd_trans of_nat_dvd_iff)+
    from ‹p⇧2 dvd x + y› have ‹[y = - x] (mod p⇧2)›
      by (simp add: add.commute cong_iff_dvd_diff)
    hence ‹[S x y = p * x ^ (p - 1)] (mod p⇧2)›
      unfolding S_def by (fact Sophie_Germain_lemma_computation_cong_simp[OF ‹p ≠ 0›])
    moreover from ‹[x ≠ 0] (mod p)› ‹z = z' * p› ‹[z ≠ 0] (mod p⇧2)› ‹prime (int p)›
    have ‹¬ p⇧2 dvd p * x ^ (p - 1)›
      by (metis cong_0_iff dvd_mult_cancel_left mult_zero_right
          of_nat_power power2_eq_square prime_dvd_power_int)
    ultimately show False using ‹p⇧2 dvd S x y› cong_dvd_iff by blast
  qed
  ultimately have p_only_nontrivial_div :
    ‹r dvd x + y ⟹ r dvd S x y ⟹ r = 1 ∨ r = p› for r :: nat
    by (metis ‹S x y ≠ 0› dvd_0_left_iff less_one
        linorder_neqE_nat of_nat_eq_0_iff power_0 power_one_right)
      ― ‹\<^term>‹p› is therefore the only possible nontrivial common divisor.›

  define mul_x_plus_y where ‹mul_x_plus_y = multiplicity p (x + y)›
  define mul_S_x_y where ‹mul_S_x_y = multiplicity p (S x y)›
  from ‹(x + y) * S x y = z ^ p›
  have ‹(x + y) * S x y = z' ^ p * p ^ p›
    by (simp add: ‹z = z' * p› power_mult_distrib)


  have ‹mul_x_plus_y + mul_S_x_y = multiplicity p (z ^ p)›
    unfolding mul_x_plus_y_def mul_S_x_y_def
    by (metis ‹(x + y) * S x y = z ^ p› ‹S x y ≠ 0› ‹x + y ≠ 0› prime_def
        prime_elem_multiplicity_mult_distrib ‹prime (int p)›)
  also have ‹z ^ p = z' ^ p * p ^ p›
    by (simp add: ‹z = z' * p› power_mult_distrib)
  also have ‹multiplicity p … = p›
    by (metis ‹[z' ≠ 0] (mod int p)› aux_prime auxiliary_prime_def cong_0_iff
        mult_of_nat_commute multiplicity_decomposeI of_nat_eq_0_iff of_nat_power
        prime_dvd_power_iff prime_gt_0_nat ‹p ≠ 0› ‹prime (int p)›)
  finally have ‹mul_x_plus_y + mul_S_x_y = p› .

  moreover have ‹(2 ≤ mul_x_plus_y ⟶ mul_S_x_y    ≤ 1) ∧
                 (2 ≤ mul_S_x_y    ⟶ mul_x_plus_y ≤ 1)›
  proof (rule ccontr)
    assume ‹¬ ?thesis›
    hence ‹2 ≤ mul_x_plus_y ∧ 2 ≤ mul_S_x_y› by linarith
    hence ‹p⇧2 dvd (x + y) ∧ p⇧2 dvd (S x y)›
      by (simp add: mul_x_plus_y_def mul_S_x_y_def multiplicity_dvd')
    thus False
      by (metis cong_0_iff less_numeral_extra(3) one_eq_prime_power_iff
          p_dvd_z p_only_nontrivial_div pos2 ‹[z ≠ 0] (mod p⇧2)› ‹prime p›)
  qed
  ultimately consider ‹mul_x_plus_y = p› ‹mul_S_x_y = 0›
    | ‹mul_x_plus_y = 0› ‹mul_S_x_y = p›
    | ‹mul_x_plus_y = 1› ‹mul_S_x_y = p - 1›
    | ‹mul_x_plus_y = p - 1› ‹mul_S_x_y = 1›
    by (metis Nat.add_diff_assoc add_cancel_right_right add_diff_cancel_left'
        diff_is_0_eq not_less_eq_eq one_add_one plus_1_eq_Suc)
  thus False
  proof cases

    assume ‹mul_x_plus_y = p› ‹mul_S_x_y = 0›
    from ‹mul_x_plus_y = p› have ‹p dvd (x + y)›
      by (metis mul_x_plus_y_def not_dvd_imp_multiplicity_0 ‹p ≠ 0›)
    hence ‹[y = - x] (mod p)›
      by (simp add: add.commute cong_iff_dvd_diff)
    hence ‹[S x y = S x (- x)] (mod p)›
      unfolding S_def
      by (meson cong_minus_minus_iff cong_pow cong_scalar_right cong_sum)
    also have ‹S x (- x) = (∑k = 0..p - 1. x ^ (p - 1))›
      unfolding S_def
      by (rule sum.cong[OF refl], simp)
        (metis One_nat_def diff_Suc_eq_diff_pred le_add_diff_inverse2 power_add)
    also from ‹p ≠ 0› have ‹… = p * x ^ (p - 1)› by simp
    finally have ‹[S x y = p * x ^ (p - 1)] (mod p)› .
    with ‹[x ≠ 0] (mod p)› have ‹p dvd S x y› by (simp add: cong_dvd_iff)
    with ‹mul_S_x_y = 0› show False
      by (metis ‹S x y ≠ 0› mul_S_x_y_def not_one_le_zero not_prime_unit
          power_dvd_iff_le_multiplicity power_one_right ‹prime (int p)›)
  next

    assume ‹mul_x_plus_y = 0› ‹mul_S_x_y = p›
    from ‹mul_S_x_y = p› have ‹p dvd S x y›
      by (metis mul_S_x_y_def not_dvd_imp_multiplicity_0 ‹p ≠ 0›)
    from Sophie_Germain_lemma_computation[OF ‹odd p›]
    have ‹(x + y) * S x y = x ^ p + y ^ p› unfolding S_def .
    moreover from ‹p dvd S x y› have ‹[(x + y) * S x y = 0] (mod p)›
      by (simp add: cong_0_iff)
    moreover have ‹[x ^ p + y ^ p = x + y] (mod p)›
    proof (rule cong_add)
      have ‹x ^ p = x * x ^ (p - 1)›
        by (simp add: power_eq_if ‹p ≠ 0›)
      moreover from ‹[x ≠ 0] (mod p)› have ‹[x ^ (p - 1)  = 1] (mod p)›
        using cong_0_iff fermat_theorem_int ‹prime p› by blast
      ultimately show ‹[x ^ p = x] (mod p)›
        by (metis cong_scalar_left mult.right_neutral)
    next
      have ‹y ^ p = y * y ^ (p - 1)›
        by (simp add: power_eq_if ‹p ≠ 0›)
      moreover from ‹[y ≠ 0] (mod p)› have ‹[y ^ (p - 1)  = 1] (mod p)›
        using cong_0_iff fermat_theorem_int ‹prime p› by blast
      ultimately show ‹[y ^ p = y] (mod p)›
        by (metis cong_scalar_left mult.right_neutral)
    qed
    ultimately have ‹p dvd x + y›
      by (simp add: cong_0_iff cong_dvd_iff)
    with ‹mul_x_plus_y = 0› show False
      by (metis ‹x + y ≠ 0› mul_x_plus_y_def multiplicity_eq_zero_iff
          not_prime_unit ‹prime (int p)›)
  next

    define x_plus_y' where ‹x_plus_y' ≡ (x + y) div p›
    define S_x_y' where ‹S_x_y' ≡ (S x y) div p ^ (p - 1)›
    define s where ‹s ≡ x + y›
    let ?f  = ‹λk. (p choose k) * (- x) ^ k * s ^ (p - k)›
    let ?f' = ‹λk. (p choose k) * (- x) ^ k * s ^ (p - 1 - k)›

    assume ‹mul_x_plus_y = 1› ‹mul_S_x_y = p - 1›
    hence ‹s = p * x_plus_y'› ‹S x y = p ^ (p - 1) * S_x_y'›
      by (simp_all add: s_def x_plus_y'_def S_x_y'_def)
        (metis dvd_mult_div_cancel mul_x_plus_y_def multiplicity_dvd power_Suc0_right,
          metis dvd_mult_div_cancel mul_S_x_y_def multiplicity_dvd)

    from fermat have ‹s * S x y = (s - x) ^ p + x ^ p›
      by (simp add: s_def ‹(x + y) * S x y = z ^ p›)
    also have ‹s - x = (- x + s)› by simp
    also have ‹(- x + s) ^ p = (∑k≤p. (p choose k) * (- x) ^ k * s ^ (p - k))›
      by (fact binomial_ring)
    also have ‹… = (∑k∈{0..p - 1}. ?f k) + (∑k∈{p}. ?f k)›
      by (rule sum_Un_eq[symmetric])
        (auto simp add: linorder_not_le prime_gt_0_nat ‹prime p›)
    also have ‹(∑k∈{0..p - 1}. ?f k) = (∑k∈{0}. ?f k) + (∑k∈{1..p - 1}. ?f k)›
      by (rule sum_Un_eq[symmetric]) auto
    also have ‹(∑k∈{1..p - 1}. ?f k) = (∑k∈{1..p - 2}. ?f k) + (∑k∈{p - 1}. ?f k)›
      by (rule sum_Un_eq[symmetric]) (use ‹3 ≤ p› in auto)
    also have ‹(∑k∈{0}. ?f k) = s * s ^ (p - 1)› by (simp add: power_eq_if ‹p ≠ 0›)
    also have ‹(∑k∈{1..p - 2}. ?f k) = (∑k∈{1..p - 2}. s * ?f' k)›
      by (rule sum.cong, simp_all)
        (metis Suc_diff_Suc diff_less less_2_cases_iff less_zeroE
          linorder_neqE order.strict_iff_not power_Suc ‹p ≠ 0›)
    also have ‹… = s * (∑k∈{1..p - 2}. ?f' k)›
      by (simp add: mult.assoc sum_distrib_left)
    also have ‹(∑k∈{p - 1}. ?f k) = s * p * x ^ (p - 1)›
      by (simp del: One_nat_def, subst binomial_symmetric[symmetric])
        (use ‹3 ≤ p› in ‹auto simp add: ‹odd p››)
    finally have ‹s * S x y =
                  s * (s ^ (p - 1) + (∑k∈{1..p - 2}. ?f' k) + p * x ^ (p - 1))›
      by (simp add: ‹odd p› distrib_left int_distrib(4))
    hence S_x_y_developed : ‹S x y = s ^ (p - 1) + (∑k∈{1..p - 2}. ?f' k) + p * x ^ (p - 1)›
      using ‹x + y ≠ 0› mult_cancel_left unfolding s_def by blast
    have ‹[p * x ^ (p - 1) = 0] (mod p⇧2)›
    proof (rule cong_trans[OF _ cong_sym])
      show ‹[p * x ^ (p - 1) = 0 + 0 + p * x ^ (p - 1)] (mod p⇧2)› by simp
    next
      show ‹[0 = 0 + 0 + p * x ^ (p - 1)] (mod p⇧2)›
      proof (rule cong_trans)
        have ‹p ^ (p - 1) dvd S x y›
          by (simp add: ‹S x y = p ^ (p - 1) * S_x_y'›)
        moreover from ‹3 ≤ p› have ‹p ^ 2 dvd p ^ (p - 1)›
          by (auto intro: le_imp_power_dvd)
        ultimately show ‹[0 = S x y] (mod p⇧2)›
          by (metis cong_0_iff cong_sym dvd_trans of_nat_dvd_iff)
      next
        show ‹[S x y = 0 + 0 + p * x ^ (p - 1)] (mod p⇧2)›
        proof (unfold S_x_y_developed, rule cong_add[OF _ cong_refl])
          show ‹[s ^ (p - 1) + (∑k = 1..p - 2. ?f' k) = 0 + 0] (mod p⇧2)›
          proof (rule cong_add)
            have ‹p dvd s› by (simp add: ‹s = p * x_plus_y'›)
            hence ‹p ^ (p - 1) dvd s ^ (p - 1)› by (simp add: dvd_power_same)
            moreover from ‹3 ≤ p› have ‹p ^ 2 dvd p ^ (p - 1)›
              by (auto intro: le_imp_power_dvd)
            ultimately show ‹[s ^ (p - 1) = 0] (mod p⇧2)›
              by (metis cong_0_iff dvd_trans of_nat_dvd_iff)
          next
            show ‹[∑k = 1..p - 2. ?f' k = 0] (mod p⇧2)›
            proof (rule cong_trans)
              show ‹[∑k = 1..p - 2. ?f' k = ∑k∈{1..p - 2}. 0] (mod p⇧2)›
              proof (rule cong_sum)
                fix k assume ‹k ∈ {1..p - 2}›
                from ‹k ∈ {1..p - 2}› have ‹p dvd (p choose k)›
                  by (auto intro: dvd_choose_prime simp add: ‹prime p›)
                moreover from ‹k ∈ {1..p - 2}› have ‹p dvd s ^ (p - 1 - k)›
                  by (auto simp add: ‹prime p› ‹s = p * x_plus_y'› prime_dvd_power_int_iff)
                ultimately show ‹[?f' k = 0] (mod p⇧2)›
                  by (simp add: cong_0_iff mult_dvd_mono power2_eq_square)
              qed
            next
              show ‹[∑k = 1..p - 2. 0 = 0] (mod int (p⇧2))› by simp
            qed
          qed
        qed
      qed
    qed
    hence ‹p dvd x ^ (p - 1)› by (simp add: cong_iff_dvd_diff power2_eq_square ‹p ≠ 0›)
    with prime_dvd_power ‹prime (int p)› have ‹p dvd x› by blast
    with ‹[x ≠ 0] (mod p)› show False by (simp add: cong_0_iff)
  next

    define x_plus_y' where ‹x_plus_y' ≡ (x + y) div p ^ (p - 1)›
    define S_x_y' where ‹S_x_y' ≡ (S x y) div p›
    assume ‹mul_x_plus_y = p - 1› ‹mul_S_x_y = 1›
    with ‹(x + y) * S x y = z ^ p› ‹mul_x_plus_y + mul_S_x_y = p›
    have ‹x_plus_y' * S_x_y' = z' ^ p›
      unfolding x_plus_y'_def S_x_y'_def z'_def
      by (metis (no_types, opaque_lifting)
          div_mult_div_if_dvd div_power mul_S_x_y_def mul_x_plus_y_def
          multiplicity_dvd of_nat_power p_dvd_z power_add power_one_right)
    have ‹coprime x_plus_y' S_x_y'›
    proof (rule ccontr)
      assume ‹¬ coprime x_plus_y' S_x_y'›
      then obtain r where ‹prime r› ‹r dvd x_plus_y'› ‹r dvd S_x_y'›
        by (metis coprime_absorb_left not_coprime_nonzeroE
            prime_factor_int unit_imp_dvd zdvd1_eq)
      from ‹r dvd x_plus_y'› have ‹r dvd x + y›
        by (metis ‹mul_x_plus_y = p - 1› dvd_div_iff_mult dvd_mult_left
            mul_x_plus_y_def multiplicity_dvd of_nat_eq_0_iff of_nat_power
            power_not_zero ‹p ≠ 0› x_plus_y'_def)
      moreover from ‹r dvd S_x_y'› have ‹r dvd S x y›
        by (metis S_x_y'_def ‹mul_S_x_y = 1› dvd_mult_div_cancel dvd_trans
            dvd_triv_right mul_S_x_y_def multiplicity_dvd power_one_right)
      ultimately have ‹r = p›
        by (metis ‹prime r› not_prime_1 p_only_nontrivial_div
            pos_int_cases prime_gt_0_int prime_nat_int_transfer)
      with ‹[z' ≠ 0] (mod int p)› ‹r dvd S_x_y'› ‹prime p› ‹prime (int p)›
        ‹x_plus_y' * S_x_y' = z' ^ p› show False
        by (metis  cong_0_iff dvd_trans dvd_triv_right prime_dvd_power_int_iff prime_gt_0_nat)
    qed
    from prod_is_some_powerE[OF ‹coprime x_plus_y' S_x_y'› ‹x_plus_y' * S_x_y' = z' ^ p›]
    obtain a where ‹normalize x_plus_y' = a ^ p› by blast
    moreover from prod_is_some_powerE
      [OF coprime_commute[THEN iffD1, OF ‹coprime x_plus_y' S_x_y'›]]
    obtain α where ‹normalize S_x_y' = α ^ p›
      by (metis ‹x_plus_y' * S_x_y' = z' ^ p› mult.commute)
    ultimately have ‹x_plus_y' = (if 0 ≤ x_plus_y' then a ^ p else (- a) ^ p)›
      ‹S_x_y' = (if 0 ≤ S_x_y' then α ^ p else (- α) ^ p)›
      by (metis ‹odd p› abs_of_nonneg abs_of_nonpos add.inverse_inverse
          linorder_linear normalize_int_def power_minus_odd)+
    then obtain α a where ‹x_plus_y' = a ^ p› ‹S_x_y' = α ^ p› by meson

    from Sophie_Germain_lemma[OF ‹odd p› ‹prime p›, of ‹- x› z ‹- y›]
      coprime_y_z fermat ‹[x ≠ 0] (mod int p)›
    obtain b β where ‹z - y = b ^ p› ‹S z (- y) = β ^ p›
      by (simp add: S_def ‹odd p› coprime_commute) (meson cong_0_iff dvd_minus_iff)
    from Sophie_Germain_lemma[OF ‹odd p› ‹prime p›, of ‹- y› z ‹- x›]
      coprime_x_z fermat ‹[y ≠ 0] (mod p)›
    obtain c γ where ‹z - x = c ^ p› ‹S z (- x) = γ ^ p›
      by (simp add: S_def ‹odd p› coprime_commute) (meson cong_0_iff dvd_minus_iff)

    have ‹¬ p dvd c›
      by (metis ‹[x ≠ 0] (mod int p)› ‹z - x = c ^ p› cong_0_iff cong_dvd_iff
          cong_iff_dvd_diff dvd_def dvd_trans p_dvd_z power_eq_if ‹p ≠ 0›)
    have ‹¬ p dvd b›
      by (metis ‹[y ≠ 0] (mod int p)› ‹z - y = b ^ p› cong_0_iff cong_dvd_iff
          cong_iff_dvd_diff dvd_def dvd_trans p_dvd_z power_eq_if ‹p ≠ 0›)

    have ‹p dvd 2 * z - x - y›
      by (metis ‹mul_S_x_y = 1› ‹mul_x_plus_y + mul_S_x_y = p› comm_monoid_add_class.add_0 diff_diff_eq dvd_diff int_ops(2)
          mul_x_plus_y_def not_dvd_imp_multiplicity_0 one_dvd p_dvd_z prime_dvd_mult_iff wlog_keep.prime_int_p)
    hence ‹[2 * z - x - y = 0] (mod p)› by (simp add: cong_0_iff)
    also from ‹z - x = c ^ p› ‹z - y = b ^ p›
    have ‹2 * z - x - y = c ^ p + b ^ p› by presburger
    also have ‹… = c ^ (p - 1) * c + b ^ (p - 1) * b›
      by (simp add: power_eq_if ‹p ≠ 0›)
    finally have ‹[c ^ (p - 1) * c + b ^ (p - 1) * b = 0] (mod p)› .
    moreover have ‹[c ^ (p - 1) = 1] (mod p)›
      by (fact fermat_theorem_int[OF ‹prime p› ‹¬ p dvd c›])
    moreover have ‹[b ^ (p - 1) = 1] (mod p)›
      by (fact fermat_theorem_int[OF ‹prime p› ‹¬ p dvd b›])
    ultimately have ‹[1 * c + 1 * b = 0] (mod p)›
      by (meson cong_add cong_scalar_right cong_sym_eq cong_trans)
    hence ‹[c + b = 0] (mod p)› by simp
    hence ‹[b = - c] (mod p)› by (simp add: add.commute cong_iff_dvd_diff)
    hence ‹[S c b = p * c ^ (p - 1)] (mod p)›
      unfolding S_def
      by (fact Sophie_Germain_lemma_computation_cong_simp[OF ‹p ≠ 0›])
    hence ‹p dvd S c b› by (simp add: cong_dvd_iff)
    moreover from ‹[c + b = 0] (mod p)› have ‹p dvd c + b› by (simp add: cong_dvd_iff)
    moreover from Sophie_Germain_lemma_computation[OF ‹odd p›]
    have ‹c ^ p + b ^ p = (c + b) * S c b› unfolding S_def ..
    ultimately have ‹p⇧2 dvd c ^ p + b ^ p›
      by (simp add: mult_dvd_mono power2_eq_square)
    moreover have ‹p⇧2 dvd x + y›
      by (metis ‹mul_S_x_y = 1› ‹mul_x_plus_y + mul_S_x_y = p› add_0
          bot_nat_0.not_eq_extremum dvd_trans linorder_not_less
          mul_x_plus_y_def multiplicity_dvd' nat_dvd_not_less odd_even_add
          odd_p of_nat_power one_dvd prime_prime_factor ‹prime p›)
    ultimately have ‹p⇧2 dvd 2 * z›
      by (metis ‹2 * z - x - y = c ^ p + b ^ p› diff_diff_eq zdvd_zdiffD)
    moreover have ‹coprime (p⇧2) 2› by (simp add: ‹odd p›)
    ultimately have ‹p⇧2 dvd z›
      by (simp add: coprime_dvd_mult_left_iff mult.commute)
    with ‹[z ≠ 0] (mod p⇧2)› show False by (simp add: cong_0_iff)
  qed
qed



text ‹Since @{thm [show_question_marks = false] SophGer_prime_imp_auxiliary_prime},
      this result is a generalization of
      @{thm [show_question_marks = false] Sophie_Germain_theorem}.›



(*<*)
end
  (*>*)