Session Probabilistic_Prime_Tests

Theory Legendre_Symbol

(*
  File:     Legendre_Symbol.thy
  Authors:  Daniel Stüwe

  Some more facts about the Legendre symbol that are missing from HOL-Number_Theory
*)
section ‹Additional Facts about the Legendre Symbol›
theory Legendre_Symbol
imports 
  "HOL-Number_Theory.Number_Theory"
begin

lemma basic_cong[simp]: 
  fixes p :: int
  assumes "2 < p"
  shows   "[-1   1] (mod p)"
          "[ 1  -1] (mod p)"
          "[ 0   1] (mod p)"
          "[ 1   0] (mod p)"
          "[ 0  -1] (mod p)"
          "[-1   0] (mod p)"
  using assms by (simp_all add: cong_iff_dvd_diff zdvd_not_zless)

lemma [simp]: "0 < n  (a mod 2) ^ n = a mod 2" for n :: nat and a :: int
 by (metis not_mod_2_eq_0_eq_1 power_one zero_power)

lemma Legendre_in_cong_eq: 
  fixes p :: int
  assumes "p > 2" and "b  {-1,0,1}"
  shows   "[Legendre a m = b] (mod p)  Legendre a m = b"
  using assms unfolding Legendre_def by auto

lemma Legendre_p_eq_2[simp]: "Legendre a 2 = a mod 2"  
  by (clarsimp simp: Legendre_def QuadRes_def cong_iff_dvd_diff) presburger

lemma Legendre_p_eq_1[simp]: "Legendre a 1 = 0" by (simp add: Legendre_def)

lemma euler_criterion_int:
  assumes "prime p" and "2 < p" 
  shows "[Legendre a p = a^((nat p-1) div 2)] (mod p)"
  using euler_criterion assms prime_int_nat_transfer
  by (metis int_nat_eq nat_numeral prime_gt_0_int zless_nat_conj) 

lemma QuadRes_neg[simp]: "QuadRes (-p) a = QuadRes p a" unfolding QuadRes_def by auto

lemma Legendre_neg[simp]: "Legendre a (-p) = Legendre a p" unfolding Legendre_def by auto

lemma Legendre_mult[simp]:
  assumes "prime p"
  shows "Legendre (a*b) p = Legendre a p * Legendre b p"
proof -
  consider "p = 2" | "p > 2"
    using assms order_le_less prime_ge_2_int by auto

  thus ?thesis proof (cases)
  case 1
    then show ?thesis
      by (metis Legendre_p_eq_2 mod_mult_eq mod_self mult_cancel_right2
                mult_eq_0_iff not_mod_2_eq_1_eq_0 one_mod_two_eq_one)
  next
  case 2
    hence "[Legendre (a*b) p = (a*b)^((nat p-1) div 2)] (mod p)"
      using euler_criterion_int assms by blast
  
    also have "[(a*b)^((nat p-1) div 2) = a^((nat p-1) div 2) * b^((nat p-1) div 2)] (mod p)"
      by (simp add: field_simps)
      
    also have "[a^((nat p-1) div 2) * b^((nat p-1) div 2) = Legendre a p * Legendre b p] (mod p)"
      using cong_sym[OF euler_criterion_int] assms 2 cong_mult by blast
  
    finally show ?thesis using Legendre_in_cong_eq[OF 2] by (simp add: Legendre_def)
  qed
qed

lemma QuadRes_mod[simp]: "p dvd n  QuadRes p (a mod n) = QuadRes p a"
  by (simp add: mod_mod_cancel QuadRes_def cong_def)

lemma Legendre_mod[simp]: "p dvd n  Legendre (a mod n) p = Legendre a p"
  by (simp add: mod_mod_cancel Legendre_def cong_def)

lemma two_cong_0_iff: "[2 = 0] (mod p)  p = 1  p = 2" for p :: nat
  unfolding cong_altdef_nat[of 0 2 p, simplified]
  using dvd_refl prime_nat_iff two_is_prime_nat by blast

lemma two_cong_0_iff_nat: "[2 = 0] (mod int p)  p = 1  p = 2"
  unfolding cong_iff_dvd_diff
  using two_is_prime_nat prime_nat_iff int_dvd_int_iff[of p 2]
  by auto

lemma two_cong_0_iff_int: "p > 0  [2 = 0] (mod p)  p = 1  p = 2" for p :: int
  by (metis of_nat_numeral pos_int_cases semiring_char_0_class.of_nat_eq_1_iff two_cong_0_iff_nat)

lemma QuadRes_2_2 [simp, intro]: "QuadRes 2 2"
  unfolding QuadRes_def
  unfolding cong_def
  by presburger

lemma Suc_mod_eq[simp]: "[Suc a = Suc b] (mod 2) = [a = b] (mod 2)"
  using Suc_eq_plus1_left cong_add_lcancel_nat by presburger

lemma div_cancel_aux: "c dvd a  (d + a * b) div c = (d div c) + a div c * b" for a b c :: nat
  by (metis div_plus_div_distrib_dvd_right dvd_div_mult dvd_trans dvd_triv_left)

corollary div_cancel_Suc: "c dvd a  1 < c  Suc (a * b) div c = a div c * b"
  using div_cancel_aux[where d = 1] by fastforce

lemma cong_aux_eq_1: "odd p  [(p - 1) div 2 - p div 4 = (p^2 - 1) div 8] (mod 2)" for p :: nat
proof (induction p rule: nat_less_induct)
  case (1 n)

  consider "n = 1" | "n > 1" using odd_pos[OF ‹odd n] by linarith

  then show ?case proof (cases) 
    assume "n > 1"

    then obtain m where m: "m = n - 2" and m': "odd m" "m < n" using ‹odd n by simp
    then obtain b where b: "m = 2 * b + 1" using oddE by blast
  
    have IH: "[(m - 1) div 2 - m div 4 = (m^2 - 1) div 8] (mod 2)" using "1.IH" m' by simp

    have [simp]: "n = 2 * b + 1 + 2" using m n > 1 b by auto

    have *: "(n2 - 1) div 8 = ((n - 2)2 - 1) div 8 + (n - 1) div 2"
      unfolding  power2_sum power2_eq_square by simp

    have "[(n - 1) div 2 - n div 4 = (n - 2 - 1) div 2 - (n - 2) div 4 + (n - 1) div 2] (mod 2)"
      by (rule cong_sym, cases "even b") (auto simp: cong_altdef_nat div_cancel_Suc elim: oddE)
    also have "[(n - 2 - 1) div 2 - (n - 2) div 4 + (n - 1) div 2 = (n2 - 1) div 8] (mod 2)"
        using IH cong_add_rcancel_nat unfolding * m by presburger
    finally show ?thesis . 
     
  qed simp
qed

lemma cong_2_pow[intro]: "(-1 :: int)^a = (-1)^b" if "[a = b] (mod 2)" for a b :: nat
proof -
  have "even a = even b"
    by (simp add: cong_dvd_iff that) 
  then show ?thesis by auto
qed

lemma card_Int: "card (A  B) = card A - card (A - B)" if "finite A"
  by (metis Diff_Diff_Int Diff_subset card_Diff_subset finite_Diff that)

text ‹Proofs are inspired by \cite{Quadratic_Reciprocity}.›

theorem supplement2_Legendre:
  fixes p :: int
  assumes "p > 2" "prime p"
  shows "Legendre 2 p = (-1) ^ (((nat p)^2 - 1) div 8)"
proof -
  interpret GAUSS "nat p" 2
    using assms
    unfolding GAUSS_def prime_int_nat_transfer
    by (simp add: two_cong_0_iff_int)

  have "card E = card ((λx. x * 2 mod p) `
          {0<..(p - 1) div 2}  {(p - 1) div 2<..})" (is "_ = card ?A")
    unfolding E_def C_def B_def A_def image_image using assms by simp
  also have "(λx. x * 2 mod p) ` {0<..(p - 1) div 2} = ((*) 2) ` {0<..(p - 1) div 2}"
    by (intro image_cong) auto
  also have "card (  {(p - 1) div 2<..}) =
               nat ((p - 1) div 2) - card ((*) 2 ` {0<..(p - 1) div 2} - {(p - 1) div 2<..})"
    using assms by (subst card_Int) (auto simp: card_image inj_onI)
  also have "card (((*) 2) ` {0<..(p - 1) div 2} - {(p - 1) div 2<..}) = card {0 <.. (p div 4)}"
    by (rule sym, intro bij_betw_same_card[of "(*) 2"] bij_betw_imageI inj_onI)
       (insert assms prime_odd_int[of p], auto simp: image_def elim!: oddE)
  also have " = nat p div 4" using assms by simp
  also have "nat ((p - 1) div 2) - nat p div 4 =  nat ((p - 1) div 2 - p div 4)"
    using assms by (simp add: nat_diff_distrib nat_div_distrib)
  finally have "card E = " .

  then have "Legendre 2 p = (-1) ^ nat ((p - 1) div 2 - p div 4)"
    using gauss_lemma assms by simp
  also have "nat ((p - 1) div 2 - p div 4) = (nat p - 1) div 2 - nat p div 4"
    using assms by (simp add: nat_div_distrib nat_diff_distrib)
  also have "(-1) ^  = ((-1) ^ (((nat p)^2 - 1) div 8) :: int)"
    using cong_aux_eq_1[of "nat p"] odd_p by blast 
  finally show ?thesis .
qed

theorem supplement1_Legendre:
  "prime p  2 < p  Legendre (-1) p = (-1)^((p-1) div 2)"
  using euler_criterion[of p "-1"] Legendre_in_cong_eq[symmetric, of p]
  by (simp add: minus_one_power_iff)

lemma QuadRes_1_right [intro, simp]: "QuadRes p 1"
  by (metis QuadRes_def cong_def power_one)

lemma Legendre_1_left [simp]: "prime p  Legendre 1 p = 1"
  by (auto simp add: Legendre_def cong_iff_dvd_diff not_prime_unit)

lemma cong_eq_0_not_coprime: "prime p  [a = 0] (mod p)  ¬coprime a p" for a p :: int
  unfolding cong_iff_dvd_diff prime_int_iff 
  by auto

lemma not_coprime_cong_eq_0: "prime p  ¬coprime a p  [a = 0] (mod p)" for a p :: int
  unfolding cong_iff_dvd_diff
  using prime_imp_coprime[of p a]
  by (auto simp: coprime_commute)

lemma prime_cong_eq_0_iff: "prime p  [a = 0] (mod p)  ¬coprime a p" for a p :: int
  using not_coprime_cong_eq_0[of p a] cong_eq_0_not_coprime[of p a]
  by auto

lemma Legendre_eq_0_iff [simp]: "prime p  Legendre a p = 0  ¬coprime a p"
  unfolding Legendre_def by (auto simp: prime_cong_eq_0_iff)

lemma Legendre_prod_mset [simp]: "prime p  Legendre (prod_mset M) p = (q∈#M. Legendre q p)"
  by (induction M) simp_all

lemma Legendre_0_eq_0[simp]: "Legendre 0 p = 0" unfolding Legendre_def by auto

lemma Legendre_values: "Legendre p q  {1, - 1, 0}"
  unfolding Legendre_def by auto

end

Theory Algebraic_Auxiliaries

(*
  File:     Algebraic_Auxiliaries.thy
  Authors:  Daniel Stüwe

  Miscellaneous facts about algebra and number theory
*)
section ‹Auxiliary Material›
theory Algebraic_Auxiliaries
  imports 
    "HOL-Algebra.Algebra"
    "HOL-Computational_Algebra.Squarefree"
    "HOL-Number_Theory.Number_Theory"
begin
          
hide_const (open) Divisibility.prime

lemma sum_of_bool_eq_card:
  assumes "finite S"
  shows "(a  S. of_bool (P a)) = real (card {a  S . P a })"
proof -
  have "(a  S. of_bool (P a) :: real) = (a  {xS. P x}. 1)"
    using assms by (intro sum.mono_neutral_cong_right) auto
  thus ?thesis by simp
qed

lemma mod_natE:
  fixes a n b :: nat
  assumes "a mod n = b"
  shows " l. a = n * l + b"
  using assms mod_mult_div_eq[of a n] by (metis add.commute)

lemma (in group) r_coset_is_image: "H #> a = (λ x. x  a) ` H"
  unfolding r_coset_def image_def
  by blast

lemma (in group) FactGroup_order:
  assumes "subgroup H G" "finite H"
  shows "order G = order (G Mod H) * card H"
using lagrange assms unfolding FactGroup_def order_def by simp

corollary (in group) FactGroup_order_div:
  assumes "subgroup H G" "finite H"
  shows "order (G Mod H) = order G div card H" 
using assms FactGroup_order subgroupE(2)[OF ‹subgroup H G] by (auto simp: order_def)

lemma group_hom_imp_group_hom_image:
  assumes "group_hom G G h"
  shows "group_hom G (Gcarrier := h  ` carrier G) h"
  using group_hom.axioms[OF assms] group_hom.img_is_subgroup[OF assms] group.subgroup_imp_group
  by(auto intro!: group_hom.intro simp: group_hom_axioms_def hom_def)

theorem homomorphism_thm:
  assumes "group_hom G G h"
  shows "G Mod kernel G (Gcarrier := h ` carrier G) h  G carrier := h ` carrier G"
  by (intro group_hom.FactGroup_iso group_hom_imp_group_hom_image assms) simp

lemma is_iso_imp_same_card:
  assumes "H  G "
  shows "order H = order G"
proof -
  from assms obtain h where "bij_betw h (carrier H) (carrier G)"
    unfolding is_iso_def iso_def
    by blast

  then show ?thesis
    unfolding order_def 
    by (rule bij_betw_same_card)
qed

corollary homomorphism_thm_order:
  assumes "group_hom G G h" 
  shows "order (Gcarrier := h ` carrier G) * card (kernel G (Gcarrier := h ` carrier G) h) = order G "
proof -
  have "order (Gcarrier := h ` carrier G) = order (G Mod (kernel G (Gcarrier := h ` carrier G) h))"
    using is_iso_imp_same_card[OF homomorphism_thm] ‹group_hom G G h
    by fastforce

  moreover have "group G" using ‹group_hom G G h group_hom.axioms by blast

  ultimately show ?thesis
    using ‹group_hom G G h and group_hom_imp_group_hom_image[OF ‹group_hom G G h] 
    unfolding FactGroup_def
    by (simp add: group.lagrange group_hom.subgroup_kernel order_def)
qed

lemma (in group_hom) kernel_subset: "kernel G H h  carrier G"
  using subgroup_kernel G.subgroupE(1) by blast

lemma (in group) proper_subgroup_imp_bound_on_card:
  assumes "H  carrier G" "subgroup H G" "finite (carrier G)"
  shows "card H  order G div 2"
proof -
  from ‹finite (carrier G) have "finite (rcosets H)"
    by (simp add: RCOSETS_def)

  note subgroup.subgroup_in_rcosets[OF ‹subgroup H G is_group]
  then obtain J where "J  H" "J  rcosets H"
    using rcosets_part_G[OF ‹subgroup H G] and H  carrier G
    by (metis Sup_le_iff inf.absorb_iff2 inf.idem inf.strict_order_iff)

  then have "2  card (rcosets H)"
    using H  rcosets H card_mono[OF ‹finite (rcosets H), of "{H, J}"]
    by simp

  then show ?thesis
    using mult_le_mono[of 2 "card (rcosets H)" "card H" "card H"]  
    unfolding lagrange[OF ‹subgroup H G]
    by force
qed

lemma cong_exp_trans[trans]: 
  "[a ^ b = c] (mod n)  [a = d] (mod n)  [d ^ b = c] (mod n)"
  "[c = a ^ b] (mod n)  [a = d] (mod n)  [c = d ^ b] (mod n)"
  using cong_pow cong_sym cong_trans by blast+

lemma cong_exp_mod[simp]: 
  "[(a mod n) ^ b = c] (mod n)  [a ^ b = c] (mod n)"
  "[c = (a mod n) ^ b] (mod n)  [c = a ^ b] (mod n)"
  by (auto simp add: cong_def mod_simps)

lemma cong_mult_mod[simp]:
  "[(a mod n) * b = c] (mod n)  [a * b = c] (mod n)"
  "[a * (b mod n) = c] (mod n)  [a * b = c] (mod n)"
  by (auto simp add: cong_def mod_simps)

lemma cong_add_mod[simp]:
  "[(a mod n) + b = c] (mod n)  [a + b = c] (mod n)"
  "[a + (b mod n) = c] (mod n)  [a + b = c] (mod n)"
  "[iA. f i mod n = c] (mod n)  [iA. f i = c] (mod n)"
  by (auto simp add: cong_def mod_simps)

lemma cong_add_trans[trans]:
  "[a = b + x] (mod n)  [x = y] (mod n)  [a = b + y] (mod n)"
  "[a = x + b] (mod n)  [x = y] (mod n)  [a = y + b] (mod n)"
  "[b + x = a] (mod n)  [x = y] (mod n)  [b + y = a] (mod n)"
  "[x + b = a] (mod n)  [x = y] (mod n)  [y + b = a] (mod n)"
  unfolding cong_def
  using mod_simps(1, 2)
  by metis+

lemma cong_mult_trans[trans]:
  "[a = b * x] (mod n)  [x = y] (mod n)  [a = b * y] (mod n)"
  "[a = x * b] (mod n)  [x = y] (mod n)  [a = y * b] (mod n)"
  "[b * x = a] (mod n)  [x = y] (mod n)  [b * y = a] (mod n)"
  "[x * b = a] (mod n)  [x = y] (mod n)  [y * b = a] (mod n)"
  unfolding cong_def
  using mod_simps(4, 5)
  by metis+

lemma cong_diff_trans[trans]:
  "[a = b - x] (mod n)  [x = y] (mod n)  [a = b - y] (mod n)" 
  "[a = x - b] (mod n)  [x = y] (mod n)  [a = y - b] (mod n)" 
  "[b - x = a] (mod n)  [x = y] (mod n)  [b - y = a] (mod n)"
  "[x - b = a] (mod n)  [x = y] (mod n)  [y - b = a] (mod n)"
  for a :: "'a :: {unique_euclidean_semiring, euclidean_ring_cancel}"
  unfolding cong_def
  by (metis mod_diff_eq)+

lemma eq_imp_eq_mod_int: "a = b  [a = b] (mod m)" for a b :: int by simp
lemma eq_imp_eq_mod_nat: "a = b  [a = b] (mod m)" for a b :: nat by simp

lemma cong_pow_I: "a = b  [x^a = x^b](mod n)" by simp

lemma gre1I: "(n = 0  False)  (1 :: nat)  n"
  by presburger

lemma gre1I_nat: "(n = 0  False)  (Suc 0 :: nat)  n"
  by presburger

lemma totient_less_not_prime:
  assumes "¬ prime n" "1 < n"
  shows "totient n < n - 1"
  using totient_imp_prime totient_less assms
  by (metis One_nat_def Suc_pred le_less_trans less_SucE zero_le_one)

lemma power2_diff_nat: "x  y  (x - y)2 = x2 + y2 - 2 * x * y" for x y :: nat
  by (simp add: algebra_simps power2_eq_square mult_2_right)
     (meson Nat.diff_diff_right le_add2 le_trans mult_le_mono order_refl)

lemma square_inequality: "1 < n  (n + n)  (n * n)" for n :: nat
  by (metis Suc_eq_plus1_left Suc_leI mult_le_mono1 semiring_normalization_rules(4))

lemma square_one_cong_one:
  assumes "[x = 1](mod n)"
  shows "[x^2 = 1](mod n)"
  using assms cong_pow by fastforce 

lemma cong_square_alt_int:
  "prime p  [a * a = 1] (mod p)  [a = 1] (mod p)  [a = p - 1] (mod p)"
  for a p :: "'a :: {normalization_semidom, linordered_idom, unique_euclidean_ring}"
  using dvd_add_triv_right_iff[of p "a - (p - 1)"]
  by (auto simp add: cong_iff_dvd_diff square_diff_one_factored dest!: prime_dvd_multD)

lemma cong_square_alt:
  "prime p  [a * a = 1] (mod p)  [a = 1] (mod p)  [a = p - 1] (mod p)"
  for a p :: nat
  using cong_square_alt_int[of "int p" "int a"] prime_nat_int_transfer[of p] prime_gt_1_nat[of p]
  by (simp flip: cong_int_iff add: of_nat_diff)  

lemma square_minus_one_cong_one:
  fixes n x :: nat
  assumes "1 < n" "[x = n - 1](mod n)"
  shows "[x^2 = 1](mod n)"
proof -
  have "[x^2 = (n - 1) * (n - 1)] (mod n)"
    using cong_mult[OF assms(2) assms(2)]
    by (simp add: algebra_simps power2_eq_square)

  also have "[(n - 1) * (n - 1) = Suc (n * n) - (n + n)] (mod n)" 
    using power2_diff_nat[of 1 n] 1 < n
    by (simp add: algebra_simps power2_eq_square)

  also have "[Suc (n * n) - (n + n) = Suc (n * n)] (mod n)"
  proof -
    have "n * n + 0 * n = n * n" by linarith
    moreover have "n * n - (n + n) + (n + n) = n * n"
      using square_inequality[OF 1 < n] le_add_diff_inverse2 by blast
    moreover have "(Suc 0 + 1) * n = n + n"
      by simp
    ultimately show ?thesis
      using square_inequality[OF 1 < n] 
      by (metis (no_types) Suc_diff_le add_Suc cong_iff_lin_nat)
  qed

  also have "[Suc (n * n) = 1] (mod n)"
    using cong_to_1'_nat by auto

  finally show ?thesis .
qed

lemma odd_prime_gt_2_int:
 "2 < p" if "odd p" "prime p" for p :: int
  using prime_ge_2_int[OF ‹prime p] ‹odd p
  by (cases "p = 2") auto

lemma odd_prime_gt_2_nat:
 "2 < p" if "odd p" "prime p" for p :: nat
  using prime_ge_2_nat[OF ‹prime p] ‹odd p
  by (cases "p = 2") auto

lemma gt_one_imp_gt_one_power_if_coprime:
  "1  x  1 < n  coprime x n  1  x ^ (n - 1) mod n"
  by (rule gre1I) (auto simp: coprime_commute dest: coprime_absorb_left)

lemma residue_one_dvd: "a mod n = 1  n dvd a - 1" for a n :: nat
  by (fastforce intro!: cong_to_1_nat simp: cong_def)

lemma coprimeI_power_mod:
  fixes x r n :: nat
  assumes "x ^ r mod n = 1" "r  0" "n  0"
  shows "coprime x n"
proof -
  have "coprime (x ^ r mod n) n"
    using coprime_1_right x ^ r mod n = 1
    by (simp add: coprime_commute)
  thus ?thesis using r  0 n  0 by simp
qed


(* MOVE - EXTRA *)
lemma prime_dvd_choose:
  assumes "0 < k" "k < p" "prime p" 
  shows "p dvd (p choose k)"
proof -
  have "k  p" using k < p by auto 

  have "p dvd fact p" using ‹prime p by (simp add: prime_dvd_fact_iff)   
  
  moreover have "¬ p dvd fact k * fact (p - k)"
    unfolding prime_dvd_mult_iff[OF ‹prime p] prime_dvd_fact_iff[OF ‹prime p] 
    using assms by simp
  
  ultimately show ?thesis
    unfolding binomial_fact_lemma[OF k  p, symmetric]
    using assms prime_dvd_multD by blast
qed

lemma cong_eq_0_I: "(iA. [f i mod n = 0] (mod n))  [iA. f i = 0] (mod n)"
  using cong_sum by fastforce

lemma power_mult_cong:
  assumes "[x^n = a](mod m)" "[y^n = b](mod m)"
  shows "[(x*y)^n = a*b](mod m)"
  using assms cong_mult[of "x^n" a m "y^n" b] power_mult_distrib
  by metis

lemma
  fixes n :: nat
  assumes "n > 1"
  shows odd_pow_cong: "odd m  [(n - 1) ^ m = n - 1] (mod n)"
  and even_pow_cong: "even m  [(n - 1) ^ m = 1] (mod n)"
proof (induction m)
  case (Suc m)
  case 1
  with Suc have IH: "[(n - 1) ^ m = 1] (mod n)" by auto
  show ?case using 1 < n cong_mult[OF cong_refl IH] by simp
next
  case (Suc m)
  case 2
  with Suc have IH: "[(n - 1) ^ m = n - 1] (mod n)" by auto
  show ?case 
    using cong_mult[OF cong_refl IH, of "(n - 1)"] and square_minus_one_cong_one[OF 1 < n, of "n - 1"]
    by (auto simp: power2_eq_square intro: cong_trans)
qed simp_all

lemma cong_mult_uneq':
  fixes a :: "'a::{unique_euclidean_ring, ring_gcd}"
  assumes "coprime d a"
  shows "[b  c] (mod a)  [d = e] (mod a)  [b * d  c * e] (mod a)"
  using cong_mult_rcancel[OF assms]
  using cong_trans[of "b*d" "c*e" a "c*d"]
  using cong_scalar_left cong_sym by blast

lemma p_coprime_right_nat: "prime p  coprime a p = (¬ p dvd a)" for p a :: nat
  by (meson coprime_absorb_left coprime_commute not_prime_unit prime_imp_coprime_nat)


lemma squarefree_mult_imp_coprime:
  assumes "squarefree (a * b :: 'a :: semiring_gcd)"
  shows   "coprime a b"
proof (rule coprimeI)
  fix l assume "l dvd a" "l dvd b"
  then obtain a' b' where "a = l * a'" "b = l * b'"
    by (auto elim!: dvdE)
  with assms have "squarefree (l2 * (a' * b'))"
    by (simp add: power2_eq_square mult_ac)
  thus "l dvd 1" by (rule squarefreeD) auto
qed

lemma prime_divisor_exists_strong:
  fixes m :: int
  assumes "m > 1" "¬prime m"
  shows   "n k. m = n * k  1 < n  n < m  1 < k  k < m"
proof -
  from assms obtain n k where nk: "n * k > 1" "n  0" "m = n * k" "n  1" "n  0" "k  1"
    using assms unfolding prime_int_iff dvd_def by auto
  from nk have "n > 1" by linarith

  from nk assms have "n * k > 0" by simp
  with n  0 have "k > 0"
    using zero_less_mult_pos by force
  with k  1 have "k > 1" by linarith
  from nk have "n > 1" by linarith

  from k > 1 nk have "n < m" "k < m" by simp_all
  with nk k > 1 n > 1 show ?thesis by blast
qed

lemma prime_divisor_exists_strong_nat:
  fixes m :: nat
  assumes "1 < m" "¬prime m"
  shows   "p k. m = p * k  1 < p  p < m  1 < k  k < m  prime p"
proof -
  obtain p where p_def: "prime p" "p dvd m" "p  m" "1 < p" 
    using assms prime_prime_factor and prime_gt_1_nat
    by blast

  moreover define k where "k = m div p"
  with p dvd m have "m = p * k" by simp

  moreover have "p < m" 
    using p  m dvd_imp_le[OF p dvd m] and m > 1
    by simp

  moreover have "1 < k" "k < m"
    using 1 < m 1 < p and p  m
    unfolding m = p * k
    by (force intro: Suc_lessI Nat.gr0I)+

  ultimately show ?thesis using 1 < m by blast
qed

(* TODO Remove *)
lemma prime_factorization_eqI:
  assumes "p. p ∈# P  prime p" "prod_mset P = n"
  shows   "prime_factorization n = P"
  using prime_factorization_prod_mset_primes[of P] assms by simp

lemma prime_factorization_prime_elem:
  assumes "prime_elem p"
  shows   "prime_factorization p = {#normalize p#}"
proof -
  have "prime_factorization p = prime_factorization (normalize p)"
    by (metis normalize_idem prime_factorization_cong)
  also have " = {#normalize p#}"
    by (rule prime_factorization_prime) (use assms in auto)
  finally show ?thesis .
qed

lemma size_prime_factorization_eq_Suc_0_iff [simp]:
  fixes n :: "'a :: factorial_semiring_multiplicative"
  shows "size (prime_factorization n) = Suc 0  prime_elem n"
proof
  assume size: "size (prime_factorization n) = Suc 0"
  hence [simp]: "n  0" by auto
  from size obtain p where *: "prime_factorization n = {#p#}"
    by (auto elim!: size_mset_SucE)
  hence p: "p  prime_factors n" by auto

  have "prime_elem (normalize p)"
    using p by (auto simp: in_prime_factors_iff)
  also have "p = prod_mset (prime_factorization n)"
    using * by simp
  also have "normalize  = normalize n"
    by (rule prod_mset_prime_factorization_weak) auto
  finally show "prime_elem n" by simp
qed (auto simp: prime_factorization_prime_elem)
(* END TODO *)

(* TODO Move *)
lemma squarefree_prime_elem [simp, intro]:
  fixes p :: "'a :: algebraic_semidom"
  assumes "prime_elem p"
  shows   "squarefree p"
proof (rule squarefreeI)
  fix x assume "x2 dvd p"
  show "is_unit x"
  proof (rule ccontr)
    assume "¬is_unit x"
    hence "¬is_unit (x2)"
      by (simp add: is_unit_power_iff)
    from assms and this and x2 dvd p have "prime_elem (x2)"
      by (rule prime_elem_mono)
    thus False by (simp add: prime_elem_power_iff)
  qed
qed

lemma squarefree_prime [simp, intro]: "prime p  squarefree p"
  by auto

lemma not_squarefree_primepow:
  assumes "primepow n"
  shows   "squarefree n  prime n"
  using assms by (auto simp: primepow_def squarefree_power_iff prime_power_iff)

lemma prime_factorization_normalize [simp]:
  "prime_factorization (normalize n) = prime_factorization n"
  by (rule prime_factorization_cong) auto

lemma one_prime_factor_iff_primepow:
  fixes n :: "'a :: factorial_semiring_multiplicative"
  shows "card (prime_factors n) = Suc 0  primepow (normalize n)"
proof
  assume "primepow (normalize n)"
  then obtain p k where pk: "prime p" "normalize n = p ^ k" "k > 0"
    by (auto simp: primepow_def)
  hence "card (prime_factors (normalize n)) = Suc 0"
    by (subst pk) (simp add: prime_factors_power prime_factorization_prime)
  thus "card (prime_factors n) = Suc 0"
    by simp
next
  assume *: "card (prime_factors n) = Suc 0"
  from * have "(pprime_factors n. p ^ multiplicity p n) = normalize n"
    by (intro prod_prime_factors) auto
  also from * have "card (prime_factors n) = 1" by simp
  then obtain p where p: "prime_factors n = {p}"
    by (elim card_1_singletonE)
  finally have "normalize n = p ^ multiplicity p n"
    by simp
  moreover from p have "prime p" "multiplicity p n > 0"
    by (auto simp: prime_factors_multiplicity)
  ultimately show "primepow (normalize n)"
    unfolding primepow_def by blast
qed

lemma squarefree_imp_prod_prime_factors_eq:
  fixes x :: "'a :: factorial_semiring_multiplicative"
  assumes "squarefree x"
  shows   "(prime_factors x) = normalize x"
proof -
  from assms have [simp]: "x  0" by auto
  have "(pprime_factors x. p ^ multiplicity p x) = normalize x"
    by (intro prod_prime_factors) auto
  also have "(pprime_factors x. p ^ multiplicity p x) = (pprime_factors x. p)"
    using assms by (intro prod.cong refl) (auto simp: squarefree_factorial_semiring')
  finally show ?thesis by simp
qed
(* END TODO *)

end

Theory Jacobi_Symbol

(*
  File:     Jacobi_Symbol.thy
  Authors:  Daniel Stüwe, Manuel Eberl

  The Jacobi symbol, a generalisation of the Legendre symbol.
  This is used in the Solovay--Strassen test.
*)
section ‹The Jacobi Symbol›
theory Jacobi_Symbol
imports 
  Legendre_Symbol
  Algebraic_Auxiliaries
begin

text ‹
  The Jacobi symbol is a generalisation of the Legendre symbol to non-primes \cite{Legendre_Symbol, Jacobi_Symbol}.
  It is defined as
  \[\left(\frac{a}{n}\right) =
      \left(\frac{a}{p_1}\right)^{k_1} \ldots \left(\frac{a}{p_l}\right)^{k_l}\]
  where $(\frac{a}{p})$ denotes the Legendre symbol, a› is an integer, n› is an odd natural
  number and $p_1^{k_1}\ldots p_l^{k_l}$ is its prime factorisation.

  There is, however, a fairly natural generalisation to all non-zero integers for n›.
  It is less clear what a good choice for n = 0› is; Mathematica and Maxima adopt
  the convention that $(\frac{\pm 1}{0}) = 1$ and $(\frac{a}{0}) = 0$ otherwise. However,
  we chose the slightly different convention $(\frac{a}{0}) = 0$ for ‹all› a› because then
  the Jacobi symbol is completely multiplicative in both arguments without any restrictions.
›
definition Jacobi :: "int  int  int" where
  "Jacobi a n = (if n = 0 then 0 else
                  (p∈#prime_factorization n. Legendre a p))"

lemma Jacobi_0_right [simp]: "Jacobi a 0 = 0"
  by (simp add: Jacobi_def)

lemma Jacobi_mult_left [simp]: "Jacobi (a * b) n = Jacobi a n * Jacobi b n"
proof (cases "n = 0")
  case False
  have *: "{# Legendre (a * b) p          . p ∈# prime_factorization n #} =
           {# Legendre a p * Legendre b p . p ∈# prime_factorization n #}"
    by (meson Legendre_mult in_prime_factors_imp_prime image_mset_cong)

  show ?thesis using False unfolding Jacobi_def * prod_mset.distrib by auto
qed auto

lemma Jacobi_mult_right [simp]: "Jacobi a (n * m) = Jacobi a n * Jacobi a m"
  by (cases "m = 0"; cases "n = 0")
     (auto simp: Jacobi_def prime_factorization_mult)

lemma prime_p_Jacobi_eq_Legendre[intro!]: "prime p  Jacobi a p = Legendre a p"
  unfolding Jacobi_def prime_factorization_prime by simp

lemma Jacobi_mod [simp]: "Jacobi (a mod m) n = Jacobi a n" if "n dvd m"
proof -
  have *: "{# Legendre (a mod m) p . p ∈# prime_factorization n #} =
           {# Legendre a p . p ∈# prime_factorization n #}" using that
    by (intro image_mset_cong, subst Legendre_mod)
       (auto intro: dvd_trans[OF in_prime_factors_imp_dvd])
  thus ?thesis by (simp add: Jacobi_def)
qed

lemma Jacobi_mod_cong: "[a = b] (mod n)  Jacobi a n = Jacobi b n"
  by (metis Jacobi_mod cong_def dvd_refl)

lemma Jacobi_1_eq_1 [simp]: "p  0  Jacobi 1 p = 1"
  by (simp add: Jacobi_def in_prime_factors_imp_prime cong: image_mset_cong)

lemma Jacobi_eq_0_not_coprime:
  assumes "n  0" "¬coprime a n"
  shows   "Jacobi a n = 0"
proof -
  from assms have "p. p dvd gcd a n  prime p"
    by (intro prime_divisor_exists) auto
  then obtain p where p: "p dvd a" "p dvd n" "prime p"
    by auto
  hence "Legendre a p = 0" using assms
    by (auto simp: prime_int_iff)
  thus ?thesis using p assms
    unfolding Jacobi_def 
    by (auto simp: image_iff prime_factors_dvd)
qed

lemma Jacobi_p_eq_2'[simp]: "n > 0  Jacobi a (2^n) = a mod 2"
  by (auto simp add: Jacobi_def prime_factorization_prime_power)

lemma Jacobi_prod_mset[simp]: "n  0  Jacobi (prod_mset M) n = (q∈#M. Jacobi q n)"
  by (induction M) simp_all

lemma non_trivial_coprime_neq:
  "1 < a  1 < b  coprime a b  a  b" for a b :: int by auto


lemma odd_odd_even: 
  fixes a b :: int 
  assumes "odd a" "odd b"
  shows "even ((a*b-1) div 2) = even ((a-1) div 2 + (b-1) div 2)"
  using assms by (auto elim!: oddE simp: algebra_simps)

lemma prime_nonprime_wlog [case_names primes nonprime sym]:
  assumes "p q. prime p  prime q  P p q"
  assumes "p q. ¬prime p  P p q"
  assumes "p q. P p q  P q p"
  shows   "P p q"
  by (cases "prime p"; cases "prime q") (auto intro: assms)

lemma Quadratic_Reciprocity_Jacobi:
  fixes p q :: int
  assumes "coprime p q"
      and "2 < p" "2 < q"
      and "odd p" "odd q"
    shows "Jacobi p q * Jacobi q p =
           (- 1) ^ (nat ((p - 1) div 2 * ((q - 1) div 2)))"
  using assms
proof (induction "nat p" "nat q" arbitrary: p q 
         rule: measure_induct_rule[where f = "λ(a, b). a + b", split_format(complete), simplified])
  case (1 p q)
  thus ?case
  proof (induction p q rule: prime_nonprime_wlog)
    case (sym p q)
    thus ?case by (simp only: add_ac coprime_commute mult_ac) blast
  next
    case (primes p q)
    from ‹prime p ‹prime q have "prime (nat p)" "prime (nat q)" "p  q"
      using prime_int_nat_transfer primes(4) non_trivial_coprime_neq prime_gt_1_int
      by blast+

    with Quadratic_Reciprocity_int and prime_p_Jacobi_eq_Legendre
    show ?case
      using ‹prime p ‹prime q primes(5-) 
      by presburger
  next
    case (nonprime p q)
    from ¬prime p obtain a b where *: "p = a * b" "1 < b" "1 < a"
      using 2 < p prime_divisor_exists_strong[of p] by auto

    hence odd_ab: "odd a" "odd b" using ‹odd p by simp_all

    moreover have "2 < b" and "2 < a" 
      using odd_ab and * by presburger+

    moreover have "coprime a q" and "coprime b q" using ‹coprime p q 
      unfolding * by simp_all

    ultimately have IH: "Jacobi a q * Jacobi q a = (- 1) ^ nat ((a - 1) div 2 * ((q - 1) div 2))"
                        "Jacobi b q * Jacobi q b = (- 1) ^ nat ((b - 1) div 2 * ((q - 1) div 2))"
      by (auto simp: * nonprime)

    have pos: "0 < q" "0 < p" "0 < a" "0 < b" 
      using * 2 < q by simp_all

    have "Jacobi p q * Jacobi q p = (Jacobi a q * Jacobi q a) * (Jacobi b q * Jacobi q b)"
      using * by simp

    also have "... = (- 1) ^ nat ((a - 1) div 2 * ((q - 1) div 2)) *
                     (- 1) ^ nat ((b - 1) div 2 * ((q - 1) div 2))"
      using IH by presburger

    also from odd_odd_even[OF odd_ab]
    have "... = (- 1) ^ nat ((p - 1) div 2 * ((q - 1) div 2))"
      unfolding * minus_one_power_iff using 2 < q *
      by (auto simp add: even_nat_iff pos_imp_zdiv_nonneg_iff)

    finally show ?case .
  qed
qed

lemma Jacobi_values: "Jacobi p q  {1, -1, 0}"
proof (cases "q = 0")
  case False
  hence "¦Legendre p x¦ = 1" if "x ∈# prime_factorization q" "Jacobi p q  0" for x
    using that prod_mset_zero_iff Legendre_values[of p x]
    unfolding Jacobi_def is_unit_prod_mset_iff set_image_mset
    by fastforce

  then have "is_unit (prod_mset (image_mset (Legendre p) (prime_factorization q)))"
    if "Jacobi p q  0"
    using that False
    unfolding Jacobi_def is_unit_prod_mset_iff 
    by auto

  thus ?thesis by (auto simp: Jacobi_def)
qed auto

lemma Quadratic_Reciprocity_Jacobi':
  fixes p q :: int
  assumes "coprime p q"
      and "2 < p" "2 < q"
      and "odd p" "odd q"
    shows "Jacobi q p = (if p mod 4 = 3  q mod 4 = 3 then -1 else 1) * Jacobi p q"
proof -
  have aux: "a  {1, -1, 0}  c  0  a*b = c  b = c * a" for b c a :: int by auto

  from Quadratic_Reciprocity_Jacobi[OF assms] 
  have "Jacobi q p = (-1) ^ nat ((p - 1) div 2 * ((q - 1) div 2)) * Jacobi p q"
    using Jacobi_values by (fastforce intro!: aux)

  also have "(-1 :: int) ^ nat ((p - 1) div 2 * ((q - 1) div 2)) = (if even ((p - 1) div 2)  even ((q - 1) div 2) then 1 else - 1)"
    unfolding minus_one_power_iff using 2 < p 2 < q
    by (auto simp: even_nat_iff)

  also have "... = (if p mod 4 = 3  q mod 4 = 3 then -1 else 1)"
    using ‹odd p ‹odd q by presburger

  finally show ?thesis .

qed

lemma dvd_odd_square: "8 dvd a2 - 1" if "odd a" for a :: int 
proof -
  obtain x where "a = 2*x + 1" using ‹odd a by (auto elim: oddE)
  thus ?thesis
    by(cases "odd x") 
      (auto elim: oddE simp: power2_eq_square algebra_simps)
qed 

lemma odd_odd_even': 
  fixes a b :: int 
  assumes "odd a" "odd b"
  shows "even (((a * b)2 - 1) div 8)  even (((a2 - 1) div 8) + ((b2 - 1) div 8))"
proof -
  obtain x where [simp]: "a = 2*x + 1" using ‹odd a by (auto elim: oddE)
  obtain y where [simp]: "b = 2*y + 1" using ‹odd b by (auto elim: oddE)
  show ?thesis
    by (cases "even x"; cases "even y"; elim oddE evenE)
       (auto simp: power2_eq_square algebra_simps)
qed

lemma odd_odd_even_nat': 
  fixes a b :: nat 
  assumes "odd a" "odd b"
  shows "even (((a * b)2 - 1) div 8)  even (((a2 - 1) div 8) + ((b2 - 1) div 8))"
proof -
  obtain x where [simp]: "a = 2*x + 1" using ‹odd a by (auto elim: oddE)
  obtain y where [simp]: "b = 2*y + 1" using ‹odd b by (auto elim: oddE)
  show ?thesis
    by (cases "even x"; cases "even y"; elim oddE evenE)
       (auto simp: power2_eq_square algebra_simps)
qed

lemma supplement2_Jacobi: "odd p  p > 1  Jacobi 2 p = (- 1) ^ (((nat p)2 - 1) div 8)"
proof (induction p rule: prime_divisors_induct)
  case (factor p x)

  then have "odd x" by force

  have "2 < p" 
    using ‹odd (p * x) prime_gt_1_int[OF ‹prime p] 
    by (cases "p = 2") auto

  have "odd p" using prime_odd_int[OF ‹prime p 2 < p] .

  have "0 < x"
    using 1 < (p * x) prime_gt_0_int[OF ‹prime p]
    and less_trans less_numeral_extra(1) zero_less_mult_pos by blast

  have base_case : "Jacobi 2 p = (- 1) ^ (((nat p)2 - 1) div 8)" 
    using 2 < p ‹prime p supplement2_Legendre and prime_p_Jacobi_eq_Legendre
    by presburger

  show ?case proof (cases "x = 1")
    case True
    thus ?thesis using base_case by force
  next
    case False
    have "Jacobi 2 (p * x) = Jacobi 2 p * Jacobi 2 x"
      using 2 < p 0 < x by simp

    also have "Jacobi 2 x = (- 1) ^ (((nat x)2 - 1) div 8)"
      using ‹odd x 0 < x x  1 by (intro factor.IH) auto

    also note base_case

    also have "(-1) ^ (((nat p)2 - 1) div 8) * (-1) ^ (((nat x)2 - 1) div 8)
             = (-1 :: int) ^ (((nat (p * x))2 - 1) div 8)" 
      unfolding minus_one_power_iff
      using 2 < p 0 < x ‹odd x ‹odd p and odd_odd_even_nat'
      using [[linarith_split_limit = 0]]
      by (force simp add: nat_mult_distrib even_nat_iff)

    finally show ?thesis .
  qed
qed simp_all

lemma mod_nat_wlog [consumes 1, case_names modulo]:
  fixes P :: "nat  bool"
  assumes "b > 0"
  assumes "k. k  {0..<b}  n mod b = k  P n"
  shows   "P n"
  using assms and mod_less_divisor 
  by fastforce

lemma mod_int_wlog [consumes 1, case_names modulo]:
  fixes P :: "int  bool"
  assumes "b > 0"
  assumes "k. 0  k  k < b  n mod b = k  P n"
  shows   "P n"
  using assms and pos_mod_conj by blast 

lemma supplement2_Jacobi':
  assumes "odd p" and "p > 1"
  shows "Jacobi 2 p = (if p mod 8 = 1  p mod 8 = 7 then 1 else -1)"
proof -
  have "0 < (4 :: nat)" by simp
  then have *: "even ((p2 - 1) div 8) = (p mod 8 = 1  p mod 8 = 7)" if "odd p" for p :: nat
  proof(induction p rule: mod_nat_wlog)
    case (modulo k)
    then consider "p mod 4 = 1" | "p mod 4 = 3"
      using ‹odd p
      by (metis dvd_0_right even_even_mod_4_iff even_numeral mod_exhaust_less_4)

    then show ?case proof (cases)
      case 1
      then obtain l where l: "p = 4 * l + 1" using mod_natE by blast
      have "even l = ((4 * l + 1) mod 8 = 1  (4 * l + 1) mod 8 = 7)" by presburger
      thus ?thesis by (simp add: l power2_eq_square algebra_simps)
    next
      case 2
      then obtain l where l: "p = 4 * l + 3" using mod_natE by blast
      have "odd l = ((3 + l * 4) mod 8 = Suc 0  (3 + l * 4) mod 8 = 7)" by presburger
      thus ?thesis by (simp add: l power2_eq_square algebra_simps)
    qed
  qed

  have [simp]: "nat p mod 8 = nat (p mod 8)"
    using p > 1 using nat_mod_distrib[of p 8] by simp
  from assms have "odd (nat p)" by (simp add: even_nat_iff)
  show ?thesis
    unfolding supplement2_Jacobi[OF assms]
              minus_one_power_iff *[OF ‹odd (nat p)]
    by (simp add: nat_eq_iff)
qed

theorem supplement1_Jacobi:
  "odd p  1 < p  Jacobi (-1) p = (-1) ^ (nat ((p - 1) div 2))"
proof (induction p rule: prime_divisors_induct)
  case (factor p x)
  then have "odd x" by force

  have "2 < p" 
    using ‹odd (p * x) prime_gt_1_int[OF ‹prime p]
    by (cases "p = 2") auto

  have "prime (nat p)"
    using ‹prime p prime_int_nat_transfer
    by blast

  have "Jacobi (-1) p = Legendre (-1) p"
    using prime_p_Jacobi_eq_Legendre[OF ‹prime p] .

  also have "... = (-1) ^ ((nat p - 1) div 2)"
    using ‹prime p 2 < p and supplement1_Legendre[of "nat p"]
    by (metis int_nat_eq nat_mono_iff nat_numeral_as_int prime_gt_0_int prime_int_nat_transfer) 

  also have "((nat p - 1) div 2) = nat ((p - 1) div 2)" by force

  finally have base_case: "Jacobi (-1) p = (-1) ^ nat ((p - 1) div 2)" .

  show ?case proof (cases "x = 1")
    case True
    then show ?thesis using base_case by simp
  next
    case False
    have "0 < x" 
      using 1 < (p * x) prime_gt_0_int[OF ‹prime p]
      by (meson int_one_le_iff_zero_less not_less not_less_iff_gr_or_eq zero_less_mult_iff)
  
    have "odd p" using ‹prime p 2 < p by (simp add: prime_odd_int) 

    have "Jacobi (-1) (p * x) = Jacobi (-1) p * Jacobi (-1) x"
      using 2 < p 0 < x by simp

    also note base_case

    also have "Jacobi (-1) x = (-1) ^ nat ((x - 1) div 2)"
      using 0 < x False ‹odd x factor.IH 
      by fastforce

    also have "(- 1) ^ nat ((p - 1) div 2) * (- 1) ^ nat ((x - 1) div 2) =
               (- 1 :: int) ^ nat ((p*x - 1) div 2)"
      unfolding minus_one_power_iff
      using 2 < p 0 < x and ‹odd x ‹odd p
      by (fastforce elim!: oddE simp: even_nat_iff algebra_simps)

    finally show ?thesis .
  qed
qed simp_all

theorem supplement1_Jacobi':
  "odd n  1 < n  Jacobi (-1) n = (if n mod 4 = 1 then 1 else -1)"
  by (simp add: even_nat_iff minus_one_power_iff supplement1_Jacobi)
     presburger?

lemma Jacobi_0_eq_0: "¬is_unit n  Jacobi 0 n = 0"
  by (cases "prime_factorization n = {#}")
     (auto simp: Jacobi_def prime_factorization_empty_iff image_iff intro: Nat.gr0I)

lemma is_unit_Jacobi_aux: "is_unit x  Jacobi a x = 1"
  unfolding Jacobi_def using prime_factorization_empty_iff[of x] by auto

lemma is_unit_Jacobi[simp]: "Jacobi a 1 = 1" "Jacobi a (-1) = 1"
  using is_unit_Jacobi_aux by simp_all

lemma Jacobi_neg_right [simp]:
  "Jacobi a (-n) = Jacobi a n"
proof -
  have * : "-n = (-1) * n" by simp
  show ?thesis unfolding *
    by (subst Jacobi_mult_right) auto
qed

lemma Jacobi_neg_left:
  assumes "odd n" "1 < n" 
  shows   "Jacobi (-a) n = (if n mod 4 = 1 then 1 else -1) * Jacobi a n"
proof -
  have * : "-a = (-1) * a" by simp
  show ?thesis unfolding * Jacobi_mult_left supplement1_Jacobi'[OF assms] ..
qed

function jacobi_code :: "int  int  int" where
"jacobi_code a n = ( 
        if n = 0 then 0
   else if n = 1 then 1
   else if a = 1 then 1
   else if n < 0 then jacobi_code a (-n)
   else if even n then if even a then 0 else jacobi_code a (n div 2)
   else if a < 0 then (if n mod 4 = 1 then 1 else -1) * jacobi_code (-a) n
   else if a = 0 then 0
   else if a  n then jacobi_code (a mod n) n
   else if even a      then (if n mod 8  {1, 7} then 1 else -1) * jacobi_code (a div 2) n
   else if coprime a n then (if n mod 4 = 3  a mod 4 = 3 then -1 else 1) * jacobi_code n a
   else 0)"
  by auto
termination
proof (relation "measure (λ(a, n). nat(abs(a) + abs(n)*2) + 
                   (if n < 0 then 1 else 0) + (if a < 0 then 1 else 0))", goal_cases)
  case (5 a n)
  thus ?case by (fastforce intro!: less_le_trans[OF pos_mod_bound])
qed auto

lemmas [simp del] = jacobi_code.simps

lemma Jacobi_code [code]: "Jacobi a n = jacobi_code a n"
proof (induction a n rule: jacobi_code.induct)
  case (1 a n)
  show ?case
  proof (cases "n = 0")
    case 2: False
    then show ?thesis proof (cases "n = 1")
      case 3: False
      then show ?thesis proof (cases "a = 1")
        case 4: False
          then show ?thesis proof (cases "n < 0")
            case True
            then show ?thesis using 2 3 4 1(1) by (subst jacobi_code.simps) simp
            next
            case 5: False
            then show ?thesis proof (cases "even n")
              case True
              then show ?thesis using 2 3 4 5 1(2)
                by (elim evenE, subst jacobi_code.simps) (auto simp: prime_p_Jacobi_eq_Legendre)
            next
              case 6: False
              then show ?thesis  proof (cases "a < 0")
                case True
                then show ?thesis using 2 3 4 5 6
                  by(subst jacobi_code.simps, subst 1(3)[symmetric]) (simp_all add: Jacobi_neg_left)
              next
                case 7: False
                then show ?thesis proof (cases "a = 0")
                  case True
                  have *: "¬ is_unit n" using 3 5 by simp
                  then show ?thesis
                    using Jacobi_0_eq_0[OF *] 2 3 4 5 7 True
                    by (subst jacobi_code.simps) simp
                next
                  case 8: False
                  then show ?thesis proof (cases "a  n")
                    case True
                    then show ?thesis using 2 3 4 5 6 7 8 1(4)
                      by (subst jacobi_code.simps) simp
                  next
                    case 9: False
                    then show ?thesis proof (cases "even a")
                      case True
                      hence "a = 2 * (a div 2)" by simp
                      also have "Jacobi  n = Jacobi 2 n * Jacobi (a div 2) n"
                        by simp
                      also have "Jacobi (a div 2) n = jacobi_code (a div 2) n"
                        using 2 3 4 5 6 7 8 9 True by (intro 1(5))
                      also have "Jacobi 2 n = (if n mod 8  {1, 7} then 1 else - 1)"
                        using 2 3 5 supplement2_Jacobi'[OF 6] by simp
                      also have " * jacobi_code (a div 2) n = jacobi_code a n"
                        using 2 3 4 5 6 7 8 9 True
                        by (subst (2) jacobi_code.simps) (simp only: if_False if_True HOL.simp_thms)
                      finally show ?thesis .
                    next
                      case 10: False
                      note foo = 1 2 3
                      then show ?thesis proof (cases "coprime a n")
                        case True
                        note this_case = 2 3 4 5 6 7 8 9 10 True
                        have "2 < a" using 10 4 7 by presburger
                        moreover have "2 < n" using 3 5 6 by presburger
                        ultimately have "jacobi_code a n = (if n mod 4 = 3  a mod 4 = 3 then - 1 else 1)
                                                        * jacobi_code n a"
                          using this_case by (subst jacobi_code.simps) simp
                        also have "jacobi_code n a = Jacobi n a"
                          using this_case by (intro 1(6) [symmetric]) auto
                        also have "(if n mod 4 = 3  a mod 4 = 3 then -1 else 1) *  = Jacobi a n"
                          using this_case and 2 < a
                          by (intro Quadratic_Reciprocity_Jacobi' [symmetric])
                             (auto simp: coprime_commute)
                        finally show ?thesis ..
                      next
                        case False
                        have *: "0 < a" "0 < n" using 5 7 8 9 by linarith+ 
                        show ?thesis
                          using 1 2 3 4 5 6 7 8 9 10 False *
                          by (subst jacobi_code.simps) (auto simp: Jacobi_eq_0_not_coprime)
                      qed
                    qed
                  qed
                qed
              qed
            qed
        qed
      qed (subst jacobi_code.simps, simp)
    qed (subst jacobi_code.simps, simp)
  qed (subst jacobi_code.simps, simp)
qed

lemma Jacobi_eq_0_imp_not_coprime:
  assumes "p  0" "p  1"
  shows   "Jacobi n p = 0  ¬coprime n p"
  using assms Jacobi_mod_cong coprime_iff_invertible_int by force

lemma Jacobi_eq_0_iff_not_coprime:
  assumes "p  0" "p  1"
  shows "Jacobi n p = 0  ¬coprime n p"
proof -
  from assms and Jacobi_eq_0_imp_not_coprime 
  show ?thesis using Jacobi_eq_0_not_coprime by auto
qed

end

Theory Residues_Nat

(*
  File:     Residues_Nat.thy
  Authors:  Daniel Stüwe, Manuel Eberl

  The multiplicative group of the ring of residues modulo n.
*)
section ‹Residue Rings of Natural Numbers›
theory Residues_Nat
  imports Algebraic_Auxiliaries
begin            

subsection ‹The multiplicative group of residues modulo n›

definition Residues_Mult :: "'a :: {linordered_semidom, euclidean_semiring}  'a monoid" where
  "Residues_Mult p =
     carrier = {x  {1..p} . coprime x p}, monoid.mult = λx y. x * y mod p, one = 1"

locale residues_mult_nat =
  fixes n :: nat and G
  assumes n_gt_1: "n > 1"
  defines "G  Residues_Mult n"
begin

lemma carrier_eq [simp]: "carrier G = totatives n"
  and mult_eq [simp]:    "(x G y) = (x * y) mod n"
  and one_eq [simp]:     "𝟭G = 1"
  by (auto simp: G_def Residues_Mult_def totatives_def)

lemma mult_eq': "(⊗G) = (λx y. (x * y) mod n)"
  by (intro ext; simp)+

sublocale group G
proof(rule groupI, goal_cases)
  case (1 x y)
  from 1 show ?case using n_gt_1
    by (auto intro!: Nat.gr0I simp: coprime_commute coprime_dvd_mult_left_iff
                                    coprime_absorb_left nat_dvd_not_less totatives_def)
next
  case (5 x)
  hence "(y. y  0  y < n  [x * y = Suc 0] (mod n))"
    using coprime_iff_invertible'_nat[of n x] n_gt_1
    by (auto simp: totatives_def)
  then obtain y where y: "y  0" "y < n" "[x * y = Suc 0] (mod n)" by blast

  from [x * y = Suc 0] (mod n) have "gcd (x * y) n = 1"
    by (simp add: cong_gcd_eq)
  hence "coprime y n" by fastforce

  with y n_gt_1 show "ycarrier G. y G x = 𝟭G"
    by (intro bexI[of _ y]) (auto simp: totatives_def cong_def mult_ac intro!: Nat.gr0I)
qed (use n_gt_1 in auto simp: mod_simps algebra_simps totatives_less›)

sublocale comm_group
  by unfold_locales (auto simp: mult_ac)

lemma nat_pow_eq [simp]: "x [^]G (k :: nat) = (x ^ k) mod n"
  using n_gt_1 by (induction k) (simp_all add: mod_mult_left_eq mod_mult_right_eq mult_ac)

lemma nat_pow_eq': "([^]G) = (λx k. (x ^ k) mod n)"
  by (intro ext) simp

lemma order_eq: "order G = totient n"
  by (simp add: order_def totient_def)

lemma order_less: "¬prime n  order G < n - 1"
  using totient_less_not_prime[of n] n_gt_1
  by (auto simp: order_eq)

lemma ord_residue_mult_group:
  assumes "a  totatives n"
  shows   "local.ord a = Pocklington.ord n a"
proof (rule dvd_antisym)
  have "[a ^ local.ord a = 1] (mod n)"
    using pow_ord_eq_1[of a] assms by (auto simp: cong_def)
  thus "Pocklington.ord n a dvd local.ord a"
    by (subst (asm) ord_divides)
next
  show "local.ord a dvd Pocklington.ord n a"
    using assms Pocklington.ord[of a n] n_gt_1 pow_eq_id by (simp add: cong_def)
qed

end


subsection ‹The ring of residues modulo n›

definition Residues_nat :: "nat  nat ring" where
  "Residues_nat m = carrier = {0..<m}, monoid.mult = λx y. (x * y) mod m, one = 1,
                     ring.zero = 0, add = λx y. (x + y) mod m"

locale residues_nat =
  fixes n :: nat and R
  assumes n_gt_1: "n > 1"
  defines "R  Residues_nat n"
begin

lemma carrier_eq [simp]: "carrier R = {0..<n}"
  and mult_eq [simp]: "x R y = (x * y) mod n"
  and add_eq [simp]: "x R y = (x + y) mod n"
  and one_eq [simp]: "𝟭R = 1"
  and zero_eq [simp]: "𝟬R = 0"
  by (simp_all add: Residues_nat_def R_def)

lemma mult_eq': "(⊗R) = (λx y. (x * y) mod n)"
  and add_eq': "(⊕R) = (λx y. (x + y) mod n)"
  by (intro ext; simp)+

sublocale abelian_group R
proof(rule abelian_groupI, goal_cases)
  case (1 x y)
  then show ?case
    using n_gt_1
    by (auto simp: mod_simps algebra_simps simp flip: less_Suc_eq_le)
next
  case (6 x)
  { assume "x < n" "1 < n"
    hence "n - x  {0..<n}" "((n - x) + x) mod n = 0" if "x  0"
      using that by auto
    moreover have "0  {0..<n}" "(0 + x) mod n = 0" if "x = 0"
      using that n_gt_1 by auto
    ultimately have "y{0..<n}. (y + x) mod n = 0"
      by meson
  }

  with 6 show ?case using n_gt_1 by auto
qed (use n_gt_1 in auto simp add: mod_simps algebra_simps)

sublocale comm_monoid R
  using n_gt_1 by unfold_locales (auto simp: mult_ac mod_simps)

sublocale cring R
  by unfold_locales (auto simp: mod_simps algebra_simps)

lemma Units_eq: "Units R = totatives n"
proof safe
  fix x assume x: "x  Units R"
  then obtain y where y: "[x * y = 1] (mod n)"
    using n_gt_1 by (auto simp: Units_def cong_def)
  hence "coprime x n"
    using cong_imp_coprime cong_sym coprime_1_left coprime_mult_left_iff by metis
  with x show "x  totatives n" by (auto simp: totatives_def Units_def intro!: Nat.gr0I)
next
  fix x assume x: "x  totatives n"
  then obtain y where "y < n" "[x * y = 1] (mod n)"
    using coprime_iff_invertible'_nat[of n x] by (auto simp: totatives_def)
  with x show "x  Units R"
    using n_gt_1 by (auto simp: Units_def mult_ac cong_def totatives_less)
qed

sublocale units: residues_mult_nat n "units_of R"
proof unfold_locales
  show "units_of R  Residues_Mult n"
    by (auto simp: units_of_def Units_eq Residues_Mult_def totatives_def Suc_le_eq mult_eq')
qed (use n_gt_1 in auto) 

lemma nat_pow_eq [simp]: "x [^]R (k :: nat) = (x ^ k) mod n"
  using n_gt_1 by (induction k) (auto simp: mod_simps mult_ac)

lemma nat_pow_eq': "([^]R) = (λx k. (x ^ k) mod n)"
  by (intro ext) simp

end


subsection ‹The ring of residues modulo a prime›

locale residues_nat_prime =
  fixes p :: nat and R
  assumes prime_p: "prime p"
  defines "R  Residues_nat p"
begin

sublocale residues_nat p R
  using prime_gt_1_nat[OF prime_p] by unfold_locales (auto simp: R_def)

lemma carrier_eq' [simp]: "totatives p = {0<..<p}"
  using prime_p by (auto simp: totatives_prime)

lemma order_eq: "order (units_of R) = p - 1"
  using prime_p by (simp add: units.order_eq totient_prime)

lemma order_eq' [simp]: "totient p = p - 1"
  using prime_p by (auto simp: totient_prime)

sublocale field R
proof (rule cring_fieldI)
  show "Units R = carrier R - {𝟬R}"
    by (subst Units_eq) (use prime_p in auto simp: totatives_prime›)
qed

lemma residues_prime_cyclic: "x{0<..<p}. {0<..<p} = {y. i. y = x ^ i mod p}"
proof -
  from n_gt_1 have "{0..<p} - {0} = {0<..<p}" by auto
  thus ?thesis using finite_field_mult_group_has_gen by simp
qed

lemma residues_prime_cyclic': "x{0<..<p}. units.ord x = p - 1"
proof -
  from residues_prime_cyclic obtain x
    where x: "x  {0<..<p}" "{0<..<p} = {y. i. y = x ^ i mod p}" by metis
  have "units.ord x = p - 1"
  proof (intro antisym)
    show "units.ord x  p - 1"
      using units.ord_dvd_group_order[of x] x(1) by (auto simp: units.order_eq intro!: dvd_imp_le)
  next
    (* TODO FIXME: a bit ugly; could be simplified if we had a theory of finite cyclic rings *)
    have "p - 1 = card {0<..<p}" by simp
    also have "{0<..<p} = {y. i. y = x ^ i mod p}" by fact
    also have "card   card ((λi. x ^ i mod p) ` {..<units.ord x})"
    proof (intro card_mono; safe?)
      fix j :: nat
      have "j = units.ord x * (j div units.ord x) + (j mod units.ord x)"
        by simp
      also have "x [^]units_of R  = x [^]units_of R (units.ord x * (j div units.ord x))
                   units_of R x [^]units_of R (j mod units.ord x)"
        using x by (subst units.nat_pow_mult) auto
      also have "x [^]units_of R (units.ord x * (j div units.ord x)) =
                   (x [^]units_of R units.ord x) [^]units_of R (j div units.ord x)"
        using x by (subst units.nat_pow_pow) auto
      also have "x [^]units_of R units.ord x = 1"
        using x(1) by (subst units.pow_ord_eq_1) auto
      finally have "x ^ j mod p = x ^ (j mod units.ord x) mod p" using n_gt_1 by simp
      thus "x ^ j mod p  (λi. x ^ i mod p) ` {..<units.ord x}"
        using units.ord_ge_1[of x] x(1) by force
    qed auto
    also have "  card {..<units.ord x}"
      by (intro card_image_le) auto
    also have " = units.ord x" by simp
    finally show "p - 1  units.ord x" .
  qed
  with x show ?thesis by metis
qed

end


subsection -1› in residue rings›

lemma minus_one_cong_solve_weak:
  fixes n x :: nat
  assumes "1 < n" "x  totatives n" "y  totatives n"
    and  "[x = n - 1] (mod n)" "[x * y = 1] (mod n)"
  shows "y = n - 1"
proof -
  define G where "G = Residues_Mult n"
  interpret residues_mult_nat n G
    by unfold_locales (use n > 1 in simp_all add: G_def›)
  have "[x * (n - 1) = x * n - x] (mod n)"
    by (simp add: algebra_simps)
  also have "[x * n - x = (n - 1) * n - (n - 1)] (mod n)"
    using assms by (intro cong_diff_nat cong_mult) auto
  also have "(n - 1) * n - (n - 1) = (n - 1) ^ 2"
    by (simp add: power2_eq_square algebra_simps)
  also have "[(n - 1)2 = 1] (mod n)"
    using assms by (intro square_minus_one_cong_one) auto
  finally have "x * (n - 1) mod n = 1"
    using n > 1 by (simp add: cong_def)
  hence "y = n - 1" 
    using inv_unique'[of x "n - 1"] inv_unique'[of x y] minus_one_in_totatives[of n] assms(1-3,5)
    by (simp_all add: mult_ac cong_def)
  then show ?thesis by simp
qed

lemma coprime_imp_mod_not_zero:
  fixes n x :: nat
  assumes "1 < n" "coprime x n"
  shows "0 < x mod n"
  using assms coprime_0_left_iff nat_dvd_not_less by fastforce

lemma minus_one_cong_solve:
  fixes n x :: nat
  assumes "1 < n"
    and eq: "[x = n - 1] (mod n)" "[x * y = 1] (mod n)"
    and coprime: "coprime x n" "coprime y n"
  shows "[y = n - 1](mod n)"
proof -
  have "0 < x mod n" "0 < y mod n"
    using coprime coprime_imp_mod_not_zero 1 < n by blast+
  moreover have "x mod n < n" "y mod n < n"
    using 1 < n by auto
  moreover have "[x mod n = n - 1] (mod n)" "[x mod n * (y mod n) = 1] (mod n)"
    using eq by auto
  moreover have "coprime (x mod n) n" "coprime (y mod n) n"
    using coprime coprime_mod_left_iff 1 < n by auto
  ultimately have "[y mod n = n - 1] (mod n)"
    using minus_one_cong_solve_weak[OF 1 < n, of "x mod n" "y mod n"]
    by (auto simp: totatives_def)
  then show ?thesis by simp
qed

corollary square_minus_one_cong_one':
  fixes n x :: nat
  assumes "1 < n"
  shows "[(n - 1) * (n - 1) = 1](mod n)"
  using square_minus_one_cong_one[OF assms, of "n - 1"] assms
  by (fastforce simp: power2_eq_square)

end

Theory QuadRes

(*
  File:     Miller_Rabin.thy
  Authors:  Daniel Stüwe

  Some facts about Quadratic Residues that are missing from the library
*)
section ‹Additional Material on Quadratic Residues›
theory QuadRes
imports 
  Jacobi_Symbol
  Algebraic_Auxiliaries
begin

text ‹Proofs are inspired by \cite{Quadratic_Residues}.›

lemma inj_on_QuadRes:
  fixes p :: int
  assumes "prime p"
  shows "inj_on (λx. x^2 mod p) {0..(p-1) div 2}"
proof 
  fix x y :: int
  assume elem: "x  {0..(p-1) div 2}" "y  {0..(p-1) div 2}"

  have * : "abs(a) < p  p dvd a  a = 0" for a :: int
    using dvd_imp_le_int by force

  assume "x2 mod p = y2 mod p"

  hence "[x2 = y2] (mod p)" unfolding cong_def .

  hence "p dvd (x2 - y2)" by (simp add: cong_iff_dvd_diff)

  hence "p dvd (x + y) * (x - y)" 
    by (simp add: power2_eq_square square_diff_square_factored) 
  
  hence "p dvd (x + y)  p dvd (x - y)"
    using ‹prime p by (simp add: prime_dvd_mult_iff) 

  moreover have "p dvd x + y  x + y = 0" "p dvd x - y  x - y = 0" 
           and "0  x" "0  y"
      using elem  
      by (fastforce intro!: * )+
  
  ultimately show "x = y" by auto
qed

lemma QuadRes_set_prime: 
  assumes "prime p" and "odd p"
  shows "{x . QuadRes p x  x  {0..<p}} = {x^2 mod p | x . x  {0..(p-1) div 2}}"
proof(safe, goal_cases)
  case (1 x)
  then obtain y where "[y2 = x] (mod p)" 
    unfolding QuadRes_def by blast

  then have A: "[(y mod p)2 = x] (mod p)" 
    unfolding cong_def
    by (simp add: power_mod)

  then have "[(-(y mod p))2 = x] (mod p)" 
    by simp

  then have B: "[(p - (y mod p))2 = x] (mod p)" 
    unfolding cong_def 
    using minus_mod_self1
    by (metis power_mod)

  have "p = 1 + ((p - 1) div 2) * 2"
    using prime_gt_0_int[OF ‹prime p] ‹odd p
    by simp

  then have C: "(p - (y mod p))  {0..(p - 1) div 2}  y mod p  {0..(p - 1) div 2}"
    using prime_gt_0_int[OF ‹prime p] 
    by (clarsimp, auto simp: le_less)

  then show ?case proof
    show ?thesis if "p - y mod p  {0..(p - 1) div 2}"
      using that B
      unfolding cong_def
      using x  {0..<p} by auto

    show ?thesis if "y mod p  {0..(p - 1) div 2}"
      using that A
      unfolding cong_def
      using x  {0..<p} by auto
  qed
qed (auto simp: QuadRes_def cong_def)

corollary QuadRes_iff: 
  assumes "prime p" and "odd p"
  shows "(QuadRes p x  x  {0..<p})  ( a  {0..(p-1) div 2}. a^2 mod p = x)"
proof -
  have "(QuadRes p x  x  {0..<p})  x  {x. QuadRes p x  x  {0..<p}}"
    by auto
  also note QuadRes_set_prime[OF assms]
  also have "(x  {x2 mod p |x. x  {0..(p - 1) div 2}}) = (a{0..(p - 1) div 2}. a2 mod p = x)"
    by blast
  finally show ?thesis .
qed

corollary card_QuadRes_set_prime:
  fixes p :: int
  assumes "prime p" and "odd p"
  shows "card {x. QuadRes p x  x  {0..<p}} = nat (p+1) div 2"
proof -
  have "card {x. QuadRes p x  x  {0..<p}} = card {x2 mod p | x . x  {0..(p-1) div 2}}"
    unfolding QuadRes_set_prime[OF assms] ..

  also have "{x2 mod p | x . x  {0..(p-1) div 2}} = (λx. x2 mod p) ` {0..(p-1) div 2}"
    by auto

  also have "card ... = card {0..(p-1) div 2}"
    using inj_on_QuadRes[OF ‹prime p] by (rule card_image)

  also have "... = nat (p+1) div 2" by simp

  finally show ?thesis .
qed

corollary card_not_QuadRes_set_prime:
  fixes p :: int
  assumes "prime p" and "odd p"
  shows "card {x. ¬QuadRes p x  x  {0..<p}} = nat (p-1) div 2"
proof -
  have "{0..<p}  {x. QuadRes p x  x  {0..<p}} = {x. QuadRes p x  x  {0..<p}}"
    by blast

  moreover have "nat p - nat (p + 1) div 2 = nat (p - 1) div 2"
    using ‹odd p prime_gt_0_int[OF ‹prime p]
    by (auto elim!: oddE simp: nat_add_distrib nat_mult_distrib)

  ultimately have "card {0..<p} - card ({0..<p}  {x. QuadRes p x  x  {0..<p}}) = nat (p - 1) div 2"
    using card_QuadRes_set_prime[OF assms] and card_atLeastZeroLessThan_int by presburger    

  moreover have "{x. ¬QuadRes p x  x  {0..<p}} = {0..<p} - {x. QuadRes p x  x  {0..<p}}"
    by blast

  ultimately show ?thesis by (auto simp add: card_Diff_subset_Int)
qed

lemma not_QuadRes_ex_if_prime:
  assumes "prime p" and "odd p"
  shows " x. ¬QuadRes p x"
proof -
  have "2 < p" using odd_prime_gt_2_int assms by blast

  then have False if "{x . ¬QuadRes p x  x  {0..<p}} = {}"
    using card_not_QuadRes_set_prime[OF assms]
    unfolding that
    by simp

  thus ?thesis by blast
qed

lemma not_QuadRes_ex:
  "1 < p  odd p  x. ¬QuadRes p x"
proof (induction p rule: prime_divisors_induct)
  case (factor p x)
  then show ?case 
    by (meson not_QuadRes_ex_if_prime QuadRes_def cong_iff_dvd_diff dvd_mult_left even_mult_iff)
qed simp_all

end

Theory Euler_Witness

(*
  File:     Euler_Witness.thy
  Author:   Daniel Stüwe

  Euler Witnesses as used in the Solovay--Strassen primality test
*)
section ‹Euler Witnesses›
theory Euler_Witness
imports 
  Jacobi_Symbol
  Residues_Nat
  QuadRes
begin

text ‹
Proofs are inspired by \cite{solovay_strassen_ori, SolovayStrassenTest, wiki:Solovay_Strassen_test, planetmath:SolovayStrassenTest}.
›
definition "euler_witness a p  [Jacobi a p  a ^ ((p - 1) div 2)] (mod p)"
abbreviation "euler_liar a p  ¬ euler_witness a p"

lemma euler_witness_mod[simp]: "euler_witness (a mod p) p = euler_witness a p"
  unfolding euler_witness_def cong_def
  by (simp add: power_mod)

lemma euler_liar_mod: "euler_liar (a mod p) p = euler_liar a p" by simp

lemma euler_liar_cong:
  assumes "[a = b] (mod p)"
  shows "euler_liar a p = euler_liar b p"
  by (metis assms cong_def euler_witness_mod)

lemma euler_witnessI: 
  "[x ^ ((n - 1) div 2) = a] (mod int n )  [Jacobi x (int n)  a] (mod int n)
      euler_witness x n"
  unfolding euler_witness_def  using cong_trans by blast

lemma euler_witnessI2:
  fixes a b :: int and k :: nat
  assumes "[a  b] (mod k)"
    and "k dvd n" 
    and main: "euler_liar x n  [Jacobi x n = a] (mod k)"
              "euler_liar x n  [x ^ ((n - 1) div 2) = b] (mod k)"
  shows "euler_witness x n"
proof (rule ccontr)
  assume "euler_liar x n"
  then have "[Jacobi x (int n) = x ^ ((n - 1) div 2)] (mod k)"
    using k dvd n cong_dvd_modulus euler_witness_def int_dvd_int_iff by blast

  moreover note main[OF ‹euler_liar x n]
  
  ultimately show False  
    using [a  b] (mod k) and cong_trans cong_sym
    by metis
qed

lemma euler_witness_exists_weak:
  assumes "odd n" "¬prime n" "2 < n"
  shows "a. euler_witness a n  coprime a n"
proof -
  show ?thesis proof (cases "squarefree n")
    case True

    obtain p k where n: "n = p * k" and "1 < p" "p < n" "1 < k" "k < n" "prime p"
      using prime_divisor_exists_strong_nat[of n] ¬prime n 2 < n by auto

    have "coprime p k" using n = p * k and ‹squarefree n
      using squarefree_mult_imp_coprime by blast

    hence "coprime (int p) (int k)" by simp

    have "odd p" using n ‹odd n by simp
    
    then obtain b where "¬QuadRes p b" 
      using not_QuadRes_ex[of "int p"]
      using ‹prime p prime_gt_1_int by auto 

    then have "[b  0] (mod p)" 
      unfolding cong_def QuadRes_def
      by (metis zero_power2)
    
    from binary_chinese_remainder_int[OF ‹coprime (int p) (int k), of b 1] 
    obtain x :: int where x: "[x = b] (mod p)" "[x = 1] (mod k)" by blast

    have "euler_witness x n"
    proof (rule euler_witnessI2[of "-1" 1 k])
      show "[x ^ ((n - 1) div 2) = 1] (mod k)" 
        using [x = 1] (mod k) and cong_def
        using cong_pow by fastforce
    next 
      have "Jacobi x n = Jacobi x p * Jacobi x k" 
        unfolding n 
        using 1 < k 1 < p by fastforce
  
      also note Jacobi_mod_cong[OF [x = b] (mod p)]
      also have "Jacobi b p = Legendre b p" 
        using ‹prime p by (simp add: prime_p_Jacobi_eq_Legendre) 
      also have "... = -1" 
        unfolding Legendre_def 
        using [b  0] (mod p) and ¬ QuadRes p b by auto
  
     also note Jacobi_mod_cong[OF [x = 1] (mod k)] 
  
     finally show "[Jacobi x (int n) = - 1] (mod int k)"
       using 1 < k by auto
    next 
      have "2 < k"
        using ‹odd n and 1 < k unfolding n 
        by(cases "k = 2") auto
        
      then show "[- 1  1] (mod k)" by auto
    next
      show "k dvd n" unfolding n by simp
    qed

    have "coprime x p"
      using [b  0] (mod int p) [x = b] (mod int p) and ‹prime p not_coprime_cong_eq_0[of p x] prime_nat_int_transfer
      by (blast intro: cong_sym cong_trans)

    then have "coprime x n"
      using n [x = 1] (mod int k)
      by (metis coprime_iff_invertible_int coprime_mult_right_iff mult.right_neutral of_nat_mult)

    then show ?thesis 
      using ‹euler_witness x n by blast

  next
    case False
    then obtain p where p: "prime p" "p^2 dvd n" using ‹odd n
      by (metis less_not_refl2 odd_pos squarefree_factorial_semiring) 
    hence "p dvd n" by auto

    from p have Z: "p dvd n div p" by (auto intro: dvd_mult_imp_div simp: power2_eq_square)

    have n: "n = p * (n div p)"
      using p(2) unfolding power2_eq_square by (metis dvd_mult_div_cancel dvd_mult_right)
    
    have "odd p" using p ‹odd n
      by (meson dvd_trans even_power pos2)
      
    then have "2 < p" using prime_ge_2_nat[OF ‹prime p] 
      by presburger  

    define a where a_def: "a = n div p + 1"

    have A: "a^p = (k=0..p. (p choose k) * (n div p)^k)" 
      unfolding binomial a_def
      using atLeast0AtMost by auto
               
    also have B: "... = 1 + (k = 1..p. (p choose k) * (n div p) ^ k)" (is "?A = 1 + ?B")
      by (simp add: sum.atLeast_Suc_atMost)

    also have C: "?B = (n div p) * (k = 1..p. (p choose k) * (n div p) ^ (k - 1))" (is "_ = _ * ?C")
      unfolding sum_distrib_left 
    proof (rule sum.cong)
      fix x
      assume "x  {1..p}"
      hence "0 < x" by simp
      hence "(n div p) ^ x = n div p * (n div p) ^ (x - 1)"
        by (metis mult.commute power_minus_mult)
      thus "(p choose x) * (n div p) ^ x = n div p * ((p choose x) * (n div p) ^ (x-1))"
        by simp
    qed simp

    finally have 1: "a ^ p = 1 + n div p * (k = 1..p. (p choose k) * (n div p) ^ (k - 1))" .

    have D: "p dvd ?C"
    proof (rule dvd_sum, goal_cases)
      fix a
      assume "a  {1..p}"
      show "p dvd (p choose a) * (n div p) ^ (a - 1)"
      proof(cases "a = p")
        note * = dvd_power_le[of _ _ 1, simplified]
        case True
        thus ?thesis 
          using p dvd n div p 2 < p by (auto intro: *)
      next
        case False
        thus ?thesis 
          using a  {1..p} and ‹prime p 
          by (auto intro!: dvd_mult2 prime_dvd_choose)
      qed
    qed
    
    then obtain m where m: "?C = p * m"
      using dvdE by blast

    have "0 < p" using ‹odd p and odd_pos by blast

    then have "?C  0"
      using   sum.atLeast_Suc_atMost[OF Suc_leI]
      by (simp add: Suc_leI sum.atLeast_Suc_atMost)

    then have "m  0" using m by fastforce 

    have "n dvd ?B"
      unfolding C m  using p by (auto simp: power2_eq_square)

    hence "?B mod n = 0" by presburger

    hence "[a^p = 1] (mod n)" 
      unfolding A B cong_def
      by (subst mod_Suc_eq[symmetric, unfolded Suc_eq_plus1_left]) simp

    have "¬ p dvd n - 1"
      using p assms(1)
      by (metis One_nat_def Suc_leI dvd_diffD1 dvd_mult_right not_prime_unit odd_pos power2_eq_square)

    have "[(n div p + 1)  1] (mod n)"
      using 2 < p ‹odd n and ‹prime p p2 dvd n
      using div_mult2_eq n nonzero_mult_div_cancel_left by (force dest!: cong_to_1_nat)

    then have "ord n a  1"
      using 2 < p ‹odd n and ‹prime p p2 dvd n
      using ord_works[of a n] 
      unfolding a_def 
      by auto

    then have "ord n a = p"
      using ord_divides[of a p n] [a ^ p = 1] (mod n) ‹prime p
      using prime_nat_iff by blast

    have "coprime n a"
      using ‹prime p p2 dvd n p dvd n div p n
      unfolding a_def
      by (metis coprime_add_one_right coprime_mult_left_iff dvd_def)

    have "[(n - 1) div 2  0] (mod p)" 
      using ¬ p dvd n - 1 ‹odd n
      by (subst cong_altdef_nat) (auto elim!: oddE)

    then have "[p  (n - 1) div 2] (mod p)"
      using cong_mult_self_right[of 1 p, simplified] cong_sym cong_trans by blast

    then have "[a^((n-1) div 2)  1] (mod n)"
      using [a ^ p = 1] (mod n)
      using order_divides_expdiff[OF ‹coprime n a, of p "(n-1) div 2", symmetric]
      unfolding ‹ord n a = p
      using cong_sym cong_trans
      by blast

    then have "[(int a)^((n-1) div 2)  1] (mod n)"
      unfolding cong_def
      by (metis of_nat_1 of_nat_eq_iff of_nat_mod of_nat_power)

    have "Jacobi a n = Jacobi a (p * int (n div p))"
      using n by (metis of_nat_mult)
      
    also have "... = Jacobi a p * Jacobi a (n div p)"
      using ‹odd n and n = p * (n div p)
      by (subst Jacobi_mult_right) auto

    also have "Jacobi a p = Jacobi 1 p"
      using p dvd n div p
      by (intro Jacobi_mod_cong) (auto simp: cong_iff_dvd_diff a_def)

    also have "... = 1"
      by (simp add: 0 < p)

    also have "Jacobi a (n div p) = Jacobi 1 (n div p)"
      by (rule Jacobi_mod_cong) (simp add: cong_iff_dvd_diff a_def)

    also have "... = 1"
      using p dvd n ‹prime p n > 2
      by (intro Jacobi_1_eq_1) (auto intro!: Nat.gr0I elim!: dvdE)

    finally show ?thesis using [int a ^ ((n - 1) div 2)  1] (mod n) ‹coprime n a
      unfolding euler_witness_def
      by (intro exI[of _ "int a"]) (auto simp: cong_sym coprime_commute)
  qed
qed

lemma euler_witness_exists:
  assumes "odd n" "¬prime n" "2 < n"
  shows "a. euler_witness a n  coprime a n  0 < a  a < n"
proof -
  obtain a where a: "euler_witness a n" "coprime a n"
    using euler_witness_exists_weak assms by blast

  then have "euler_witness (a mod n) n" "coprime (a mod n) n"
    using ‹odd n by (simp_all add: odd_pos)

  moreover have "0 < (a mod n)" 
  proof -
    have "0  (a mod n)" by (rule pos_mod_sign) (use 2 < n in simp)

    moreover have "0  (a mod n)"
    using ‹coprime (a mod n) n  coprime_0_left_iff[of "int n"] 2 < n by auto

    ultimately show ?thesis by linarith
  qed

  moreover have "a mod n < n"
    by (rule pos_mod_bound) (use 2 < n in simp)

  ultimately show ?thesis by blast
qed

lemma euler_witness_exists_nat:
  assumes "odd n" "¬prime n" "2 < n"
  shows "a. euler_witness (int a) n  coprime a n  0 < a  a < n"
  using euler_witness_exists[OF assms]
  using zero_less_imp_eq_int by fastforce

lemma euler_liar_1_p[intro, simp]: "p  0  euler_liar 1 p"
  unfolding euler_witness_def by simp

lemma euler_liar_mult:
  shows "euler_liar y n  euler_liar x n  euler_liar (x*y) n"
  unfolding euler_witness_def
  by (auto simp: power_mult_distrib intro: cong_mult)

lemma euler_liar_mult':
  assumes "1 < n" "coprime y n"
  shows "euler_liar y n  euler_witness x n  euler_witness (x*y) n"
proof(simp add: euler_witness_def power_mult_distrib, rule cong_mult_uneq', goal_cases)
case 1
  then show ?case 
    using Jacobi_eq_0_iff_not_coprime[of n y] Jacobi_values[of y n] and assms
    by auto
qed simp_all

lemma prime_imp_euler_liar:
  assumes "prime p" "2 < p" "0 < x" "x < p"
  shows   "euler_liar x p"
  using assms prime_p_Jacobi_eq_Legendre and euler_criterion
  unfolding euler_witness_def 
  by simp

locale euler_witness_context =
  fixes p :: nat
  assumes p_gt_1: "p > 1" and odd_p: "odd p"
begin

definition G where "G = Residues_Mult p"

sublocale residues_mult_nat p G
  by unfold_locales (use p_gt_1 in simp_all add: G_def›)

definition "H = {x. coprime x p  euler_liar (int x) p  x  {1..<p}}"
         
sublocale H: subgroup H G
proof -
  have subset: "H  carrier G" by (auto simp: H_def totatives_def)
  show "subgroup H G"
  proof (rule group.subgroupI, goal_cases)
    case 1
    then show ?case by (fact is_group)
  next
    case 3
    have "euler_liar 1 p" using p_gt_1
      unfolding euler_witness_def cong_def by simp
    then show ?case
      using p_gt_1 by (auto simp: H_def)
  next
    case (4 x)
    then have "coprime x p" "euler_liar x p" "1  x" "x < p"
      by (auto simp: H_def)

    define y where "y = invG x"

    from x  H› have "x G y = 𝟭G"
      unfolding y_def using subset by (intro r_inv) auto

    hence *: "(x * y) mod p = 1" by (auto simp: y_def)
    hence **: "(int x * int y) mod p = 1"
      by (metis of_nat_1 of_nat_mod of_nat_mult)

    from * have "coprime y p" 
      using p_gt_1 x < p 
      by (auto simp add: coprime_iff_invertible'_nat cong_def mult.commute)

    moreover from 4 have "y  carrier G" 
      using subset unfolding y_def by (intro inv_closed) auto

    hence "1  y" "y < p" using p_gt_1 totatives_less[of y p]
      by (auto simp: totatives_def)    

    moreover have "euler_liar 1 p"
      using p_gt_1 by (intro euler_liar_1_p) auto
    then have "euler_liar (int x * int y) p"
      using ** p_gt_1 by (subst euler_liar_cong[of "int x * int y" 1 p]) (auto simp: cong_def)

    then have "euler_liar y p"
      using ‹coprime x p ‹euler_liar x p and euler_liar_mult'[OF p_gt_1, of x y]
      by (metis coprime_int_iff mult.commute)
      
    ultimately show ?case by (auto simp: H_def simp flip: y_def) 
  next
    case (5 x y)
    then show ?case
      unfolding euler_witness_def H_def
      by (auto intro!: gre1I_nat cong_mult
               simp: coprime_commute coprime_dvd_mult_left_iff
                     nat_dvd_not_less zmod_int power_mult_distrib)
  qed fact+
qed

lemma H_finite: "finite H"
  unfolding H_def by simp

lemma euler_witness_coset:
  assumes "euler_witness x p"
  shows "y  H #>G x  euler_witness y p"
  using assms euler_liar_mult'[OF p_gt_1, of _ x] unfolding r_coset_is_image H_def  
  by (auto simp add: mult.commute of_nat_mod)

lemma euler_liar_coset:
  assumes "euler_liar x p" "x  carrier G"
  shows "y  H #>G x  euler_liar y p"
  using is_group H.rcos_const[of x] assms totatives_less[of x p] p_gt_1
  by (auto simp: H_def totatives_def)

lemma in_cosets_aux:
  assumes "euler_witness x p" "x  carrier G"
  shows "H #>G x  rcosetsG H"
  using assms by (intro rcosetsI) (auto simp: H_def totatives_def)

lemma H_cosets_aux:
  assumes "euler_witness x p"
  shows "(H #>G x)  H = {}"
  using euler_witness_coset assms unfolding H_def by blast

lemma H_rcosets_aux:
  assumes "euler_witness x p" "x  carrier G"
  shows "{H, H #>G x}  rcosetsG H"
  using in_cosets_aux[OF assms] H.subgroup_in_rcosets[OF is_group]
  by blast

lemma H_not_eq_coset:
  assumes "euler_witness x p"
  shows "H  H #>G x"
  using H.one_closed H_def assms(1) euler_witness_coset by blast

lemma finite_cosets_H: "finite (rcosetsG H)"
  using rcosets_part_G[OF H.is_subgroup] 
  by (auto intro: finite_UnionD)

lemma card_cosets_limit:
  assumes "euler_witness x p" "x  carrier G"
  shows "2  card (rcosetsG H)"
  using H_not_eq_coset[OF assms(1)] card_mono[OF finite_cosets_H H_rcosets_aux[OF assms]]
  by simp

lemma card_euler_liars_cosets_limit:
  assumes "¬prime p" "2 < p"
  shows "2  card (rcosetsG H)" "card H < (p - 1) div 2"
proof -
  have "H  rcosetsG H" " H  carrier G"
    by (auto simp: is_group H.subgroup_in_rcosets simp del: carrier_eq)

  obtain a :: nat where "euler_witness a p" "coprime a p" "0 < a" "a < p"
    using odd_p assms euler_witness_exists_nat
    by blast

  moreover have a: "a  carrier G"
    using calculation by (auto simp: totatives_def)

  ultimately show "2  card (rcosetsG H)"
    using card_cosets_limit by blast

  have "finite H"
    by (rule finite_subset[OF H.subset]) auto

  have "finite (H #>G a)"
    by (rule cosets_finite[OF rcosetsI]) (use H.subset a in auto)

  have "H #>G a  rcosetsG H"
    using H.subset rcosetsI a  carrier G› by blast

  then have "card H = card (H #>G a)"
    using card_rcosets_equal H.subset by blast

  moreover have "H  H #>G a  carrier G"
    using rcosets_part_G[OF H.is_subgroup]
    using H.subgroup_in_rcosets[OF is_group] and ‹H #>G a  rcosetsG H›
    by auto

  then have "card H + card (H #>G a)  card (carrier G)"
    using ‹finite H› ‹finite (H #>G a)
    using H_cosets_aux[OF ‹euler_witness a p]
    using H.subset finite_subset card_Un_disjoint
    by (subst card_Un_disjoint[symmetric]) (auto intro: card_mono simp flip: card_Un_disjoint)

  ultimately have "card H  card (carrier G) div 2"
    by linarith

  obtain pa k where pa: "p = pa * k" "1 < pa" "pa < p" "1 < k" "k < p" "prime pa"
    using prime_divisor_exists_strong_nat[OF p_gt_1 ¬ prime p]
    by blast

  hence "¬coprime pa p" by simp

  then have "carrier G  {1..<p}"
    using ¬ prime p pa(2, 3) by (auto simp: totatives_def)

  then have "card (carrier G) < p - 1"
    by (metis card_atLeastLessThan finite_atLeastLessThan psubset_card_mono)

  show "card H < (p - 1) div 2"
    using ‹card H  card (carrier G) div 2 ‹card (carrier G) < p - 1
    using odd_p less_mult_imp_div_less[of "card (carrier G)" "(p - 1) div 2" 2]
    by auto
qed

lemma prime_imp_G_is_H:
  assumes "prime p" "2 < p"
  shows "carrier G = H"
  unfolding H_def using assms prime_imp_euler_liar[of p] totatives_less[of _ p]
  by (auto simp: totatives_def)

end

end

Theory Carmichael_Numbers

(*
  File:     Carmichael_Numbers.thy
  Authors:  Daniel Stüwe

  Definition and basic properties of Carmichael numbers
*)
section ‹Carmichael Numbers›
theory Carmichael_Numbers
imports
  Residues_Nat
begin

text ‹
  A Carmichael number is a composite number n› that Fermat's test incorrectly labels
  as primes no matter which witness a› is chosen (except in the case that a› shares a
  factor with n›). \cite{Carmichael_numbers, wiki:Carmichael_number}
›
definition Carmichael_number :: "nat  bool" where
  "Carmichael_number n  n > 1  ¬prime n  (a. coprime a n  [a ^ (n - 1) = 1] (mod n))"

lemma Carmichael_number_0[simp, intro]: "¬Carmichael_number 0"
  unfolding Carmichael_number_def by simp

lemma Carmichael_number_1[simp, intro]: "¬Carmichael_number 1"
  by (auto simp: Carmichael_number_def)

lemma Carmichael_number_Suc_0[simp, intro]: "¬Carmichael_number (Suc 0)"
  by (auto simp: Carmichael_number_def)

lemma Carmichael_number_not_prime: "Carmichael_number n  ¬prime n"
  by (auto simp: Carmichael_number_def)

lemma Carmichael_number_gt_3: "Carmichael_number n  n > 3"
proof -
  assume *: "Carmichael_number n"
  hence "n > 1" by (auto simp: Carmichael_number_def)
  {
    assume "¬(n > 3)"
    with n > 1 have "n = 2  n = 3" by auto
    with * and Carmichael_number_not_prime[of n] have False by auto
  }
  thus "n > 3" by auto
qed

text ‹
  The proofs are inspired by \cite{Carmichael_numbers, Carmichael_number_square_free}.
›
lemma Carmichael_number_imp_squarefree_aux:
  assumes "Carmichael_number n"
  assumes n: "n = p^r * l" and "prime p" "¬p dvd l"
  assumes "r > 1"
  shows False
proof -
  have "¬ prime n" using ‹Carmichael_number n unfolding Carmichael_number_def by blast

  have * : "[a^(n-1) = 1] (mod n)" if "coprime a n" for a
    using ‹Carmichael_number n that
    unfolding Carmichael_number_def
    by blast

  have "1  n"
    unfolding n using ‹prime p ¬ p dvd l
    by (auto intro: gre1I_nat)

  have "2  n"
  proof(cases "n = 1")
    case True
    then show ?thesis
      unfolding n using 1 < r prime_gt_1_nat[OF ‹prime p]
      by simp
  next
    case False
    then show ?thesis using 1  n by linarith
  qed

  have "p < p^r"
    using prime_gt_1_nat[OF ‹prime p] 1 < r
    by (metis power_one_right power_strict_increasing_iff)

  hence "p < n" using 1  n less_le_trans n by fastforce 

  then have [simp]: "{..n} - {0..Suc 0} = {2..n}" by auto

  obtain a where a: "[a = p + 1] (mod p^r)" "[a = 1] (mod l)"
    using binary_chinese_remainder_nat[of "p^r" l "p + 1" 1] 
    and ‹prime p prime_imp_coprime_nat coprime_power_left_iff ¬p dvd l
    by blast

  hence "coprime a n"
    using lucas_coprime_lemma[of 1 a l] cong_imp_coprime[of "p+1" a "p^r"] 
      and coprime_add_one_left cong_sym 
    unfolding n = p ^ r * l coprime_mult_right_iff coprime_power_right_iff power_one_right
    by blast

  hence "[a ^ (n - 1) = 1] (mod n)"
    using * by blast

  hence "[a ^ (n - 1) = 1] (mod p^r)"
    using n cong_modulus_mult_nat by blast

  hence A: "[a ^ n = a] (mod p^r)"
    using cong_scalar_right[of "a^(n-1)" 1 "p^r" a] 1  n
    unfolding power_Suc2[symmetric]
    by simp

  have "r = Suc (Suc (r - 2))"
    using 1 < r by linarith

  then have "p^r = p^2 * p^(r-2)"
    by (simp add: algebra_simps flip: power_add power_Suc)
   
  hence "[a ^ n = a] (mod p^2)" "[a = p + 1] (mod p^2)"
    using 1 < r A cong_modulus_mult_nat [a = p + 1] (mod p^r)
    by algebra+

  hence 1: "[(p + 1) ^ n = (p + 1)] (mod p^2)"
    by (metis (mono_tags) cong_def power_mod)

  have "[(p + 1) ^ n = (kn. of_nat (n choose k) * p ^ k * 1 ^ (n - k))] (mod p^2)"
    using binomial[of p 1 n] by simp

  also have "(kn. of_nat (n choose k) * p ^ k * 1 ^ (n - k)) =
             (k = 0..1. (n choose k) * p ^ k) + (k {2..n}. of_nat (n choose k) * p ^ k * 1 ^ (n - k))"
    using 2  n finite_atMost[of n]
    by (subst sum.subset_diff[where B = "{0..1}"]) auto

  also have "[(k = 0..1. (n choose k) * p ^ k) = 1] (mod p^2)"
    by (simp add: cong_altdef_nat p ^ r = p2 * p ^ (r - 2) n)

  also have "[(k {2..n}. of_nat (n choose k) * p ^ k * 1 ^ (n - k)) = 0] (mod p^2)"
    by (rule cong_eq_0_I) (clarsimp simp: cong_0_iff le_imp_power_dvd)

  finally have 2: "[(p + 1) ^ n = 1] (mod p^2)" by simp

  from cong_trans[OF cong_sym[OF 1] 2] 
  show ?thesis
    using prime_gt_1_nat[OF ‹prime p]
    by (auto dest: residue_one_dvd[unfolded One_nat_def] simp add: cong_def numeral_2_eq_2)
qed

theorem Carmichael_number_imp_squarefree:
  assumes "Carmichael_number n"
  shows "squarefree n"
proof(rule squarefreeI, rule ccontr)
  fix x :: nat
  assume "x2 dvd n"
  from assms have "n > 0" using Carmichael_number_gt_3[of n] by auto
  from x2 dvd n and 0 < n have "0 < x" by auto

  assume "¬ is_unit x"
  then obtain p where "prime p" "p dvd x"
    using prime_divisor_exists[of x] 0 < x
    by blast

  with x2 dvd n have "p^2 dvd n"
    by auto

  obtain l where n: "n = p ^ multiplicity p n * l"
    using multiplicity_dvd[of p n] by blast

  then have "¬ p dvd l"
    using multiplicity_decompose'[where x = n and p = p]
    using ‹prime p 0 < n
    by (metis nat_dvd_1_iff_1 nat_mult_eq_cancel1 neq0_conv prime_prime_factor_sqrt zero_less_power)

  have "2  multiplicity p n"
    using p2 dvd n 0 < n prime_gt_1_nat[OF ‹prime p]
    by (auto intro!: multiplicity_geI simp: power2_eq_square)
    
  then show False 
    using Carmichael_number_imp_squarefree_aux[OF ‹Carmichael_number n n] ‹prime p ¬ p dvd l
    by auto
qed

corollary Carmichael_not_primepow:
  assumes "Carmichael_number n"
  shows   "¬primepow n"
  using Carmichael_number_imp_squarefree[of n] Carmichael_number_not_prime[of n] assms
        primepow_gt_0_nat[of n] by (auto simp: not_squarefree_primepow)

lemma Carmichael_number_imp_squarefree_alt_weak:
  assumes "Carmichael_number n"
  shows "p l. (n = p * l)  prime p  ¬p dvd l"
proof -
  from assms have "n > 1"
    using Carmichael_number_gt_3[of n] by simp
  have "squarefree n"
    using Carmichael_number_imp_squarefree assms 
    by blast

  obtain p l where "p * l = n" "prime p" "1 < p"
    using assms prime_divisor_exists_strong_nat prime_gt_1_nat
    unfolding Carmichael_number_def by blast

  then have "multiplicity p n = 1"
    using 1 < n ‹squarefree n and multiplicity_eq_zero_iff[of n p] squarefree_factorial_semiring''[of n]
    by auto

  then have "¬p dvd l"
    using 1 < n ‹prime p p * l = n multiplicity_decompose'[of n p]
    by force

  show ?thesis 
    using p * l = n ‹prime p ¬p dvd l
    by blast
qed

theorem Carmichael_number_odd:
  assumes "Carmichael_number n"
  shows   "odd n"
proof (rule ccontr)
  assume "¬ odd n"
  hence "even n" by simp
  from assms have "n  4" using Carmichael_number_gt_3[of n] by simp
  have "[(n - 1) ^ (n - 1) = n - 1] (mod n)"
    using ‹even n and n  4 by (intro odd_pow_cong) auto

  then have "[(n - 1) ^ (n - 1)  1] (mod n)"
    using cong_trans[of 1 "(n - 1) ^ (n - 1)" n "n-1", OF cong_sym] 4  n
    by (auto simp: cong_def)

  moreover have "coprime (n - 1) n"
    using n  4 coprime_diff_one_left_nat[of n] by auto

  ultimately show False
    using assms unfolding Carmichael_number_def by blast
qed

lemma Carmichael_number_imp_squarefree_alt:
  assumes "Carmichael_number n"
  shows "p l. (n = p * l)  prime p  ¬p dvd l  2 < l"
proof -
  obtain p l where [simp]: "(n = p * l)" and "prime p" "¬p dvd l"
    using Carmichael_number_imp_squarefree_alt_weak and assms by blast

  moreover have "odd n" using Carmichael_number_odd and assms by blast

  consider "l = 0  l = 2" | "l = 1" | "2 < l"
    by fastforce

  then have "2 < l" 
  proof cases
    case 1
    then show ?thesis
      using ‹odd n by auto
  next
    case 2
    then show ?thesis
      using n = p * l ‹prime p ‹Carmichael_number n
      unfolding Carmichael_number_def by simp
  qed simp

  ultimately show ?thesis by blast
qed

lemma Carmichael_number_imp_dvd:
  fixes n :: nat
  assumes Carmichael_number: "Carmichael_number n" and "prime p" "p dvd n"
  shows "p - 1 dvd n - 1"
proof -
  have "¬prime n" using Carmichael_number unfolding Carmichael_number_def by blast
  obtain u where "n = p * u" using p dvd n by blast
  have "squarefree n" using Carmichael_number_imp_squarefree assms by blast
  then have "¬p dvd u"
    using ‹prime p not_prime_unit[of p]
    unfolding power2_eq_square squarefree_def n = p * u
    by fastforce

  define R where "R = Residues_nat p"
  interpret residues_nat_prime p R
    by unfold_locales (simp_all only: ‹prime p R_def)
  
  obtain a where a: "a  {0<..<p}" "units.ord a = p - 1"
    using residues_prime_cyclic' ‹prime p by metis
  from a have "a  totatives p" by (auto simp: totatives_prime ‹prime p)

  have "coprime p u"
    using ‹prime p ¬ p dvd u
    by (simp add: prime_imp_coprime_nat) 

  then obtain x where "[x = a] (mod p)" "[x = 1] (mod u)" 
    using binary_chinese_remainder_nat[of p u a 1] by blast

  have "coprime x p"
    using a  totatives p and cong_imp_coprime[OF cong_sym[OF [x = a] (mod p)]]
    by (simp add: coprime_commute totatives_def)
  moreover have "coprime x u" 
    using coprime_1_left and cong_imp_coprime[OF cong_sym[OF [x = 1] (mod u)]] by blast
  ultimately have "coprime x n"
    by (simp add: n = p * u)

  have "[a ^ (n - 1) = x ^ (n - 1)] (mod p)"
    using [x = a] (mod p) by (intro cong_pow) (auto simp: cong_sym_eq)
  also have "[x ^ (n - 1) = 1] (mod n)"
    using Carmichael_number ‹coprime x n unfolding Carmichael_number_def by blast
  then have "[x ^ (n - 1) = 1] (mod p)"
    using n = p * u cong_modulus_mult_nat by blast 
  finally have "ord p a dvd n - 1"
    by (simp add: ord_divides [symmetric])
  also have "ord p a = p - 1"
    using a a  totatives p by (simp add: units.ord_residue_mult_group)
  finally show ?thesis .    
qed

text ‹
  The following lemma is also called Korselt's criterion.
›
lemma Carmichael_numberI:
  fixes n :: nat
  assumes "¬ prime n" "squarefree n" "1 < n" and
          DIV: "p. p  prime_factors n  p - 1 dvd n - 1"
  shows   "Carmichael_number n"
  unfolding Carmichael_number_def
proof (intro assms conjI allI impI)
  fix a :: nat assume "coprime a n"

  have n: "n = (prime_factors n)"
    using prime_factorization_nat and squarefree_factorial_semiring'[of n] 1 < n ‹squarefree n
    by fastforce

  have "x ∈# prime_factorization n  y ∈# prime_factorization n  x  y  coprime x y" for x y
    using in_prime_factors_imp_prime primes_coprime
    by blast

  moreover {
    fix p 
    assume p: "p ∈# prime_factorization n"

    have "¬p dvd a"
      using ‹coprime a n p coprime_common_divisor_nat[of a n p]
      by (auto simp: in_prime_factors_iff)
    with p have "[a ^ (p - 1) = 1] (mod p)"
      by (intro fermat_theorem) auto
    hence "ord p a dvd p - 1"
      by (subst (asm) ord_divides)
    also from p have "p - 1 dvd n - 1"
      by (rule DIV)
    finally have "[a ^ (n - 1) = 1] (mod p)"
      by (subst ord_divides)
  }

  ultimately show "[a ^ (n - 1) = 1] (mod n)"
    using n coprime_cong_prod_nat by metis
qed

theorem Carmichael_number_iff:
  "Carmichael_number n 
     n  1  ¬prime n  squarefree n  (pprime_factors n. p - 1 dvd n - 1)"
proof -
  consider "n = 0" | "n = 1" | "n > 1" by force
  thus ?thesis using Carmichael_numberI[of n] Carmichael_number_imp_dvd[of n]
    by cases (auto simp: Carmichael_number_not_prime Carmichael_number_imp_squarefree)
qed

text ‹
  Every Carmichael number has at least three distinct prime factors.
›
theorem Carmichael_number_card_prime_factors:
  assumes "Carmichael_number n"
  shows   "card (prime_factors n)  3"
proof (rule ccontr)
  from assms have "n > 3"
    using Carmichael_number_gt_3[of n] by simp
  assume "¬(card (prime_factors n)  3)"
  moreover have "card (prime_factors n)  0"
    using assms Carmichael_number_gt_3[of n] by (auto simp: prime_factorization_empty_iff)
  moreover have "card (prime_factors n)  1"
    using assms by (auto simp: one_prime_factor_iff_primepow Carmichael_not_primepow)
  ultimately have "card (prime_factors n) = 2"
    by linarith
  then obtain p q where pq: "prime_factors n = {p, q}" "p  q"
    by (auto simp: card_Suc_eq numeral_2_eq_2)
  hence "prime p" "prime q" by (auto simp: in_prime_factors_iff)

  have "n = (prime_factors n)"
    using assms by (subst squarefree_imp_prod_prime_factors_eq)
                   (auto simp: Carmichael_number_imp_squarefree)
  with pq have n_eq: "n = p * q" by simp

  have "p - 1 dvd n - 1" and "q - 1 dvd n - 1" using assms pq
    unfolding Carmichael_number_iff by blast+
  with ‹prime p ‹prime q n = p * q p  q show False
  proof (induction p q rule: linorder_wlog)
    case (le p q)
    hence "p < q" by auto
    have "[q = 1] (mod q - 1)"
      using prime_gt_1_nat[of q] ‹prime q by (simp add: cong_def le_mod_geq)
    hence "[p * q - 1 = p * 1 - 1] (mod q - 1)"
      using le prime_gt_1_nat[of p] by (intro cong_diff_nat cong_mult) auto
    hence "[p - 1 = n - 1] (mod q - 1)"
      by (simp add: n = p * q cong_sym_eq)
    also have "[n - 1 = 0] (mod q - 1)"
      using le by (simp add: cong_def)
    finally have "(p - 1) mod (q - 1) = 0"
      by (simp add: cong_def)
    also have "(p - 1) mod (q - 1) = p - 1"
      using prime_gt_1_nat[of p] ‹prime p p < q by (intro mod_less) auto
    finally show False
      using prime_gt_1_nat[of p] ‹prime p by simp
  qed (simp_all add: mult.commute)
qed

lemma Carmichael_number_iff':
  fixes n :: nat
  defines "P  prime_factorization n"
  shows "Carmichael_number n 
           n > 1  size P  1  (p∈#P. count P p = 1  p - 1 dvd n - 1)"
  unfolding Carmichael_number_iff
  by (cases "n = 0") (auto simp: P_def squarefree_factorial_semiring' count_prime_factorization)

text ‹
  The smallest Carmichael number is 561, and it was found and proven so by
  Carmichael in 1910~\cite{carmichael1910note}.
›
lemma Carmichael_number_561: "Carmichael_number 561" (is "Carmichael_number ?n")
proof -
  have [simp]: "prime_factorization (561 :: nat) = {#3, 11, 17#}"
    by (rule prime_factorization_eqI) auto
  show ?thesis by (subst Carmichael_number_iff') auto
qed

end

Theory Fermat_Witness

(*
  File:     Fermat_Witness.thy
  Authors:  Daniel Stüwe

  Fermat witnesses as used in Fermat's test
*)
section ‹Fermat Witnesses›
theory Fermat_Witness
  imports Euler_Witness Carmichael_Numbers
begin

definition divide_out :: "'a :: factorial_semiring  'a  'a × nat" where
  "divide_out p x = (x div p ^ multiplicity p x, multiplicity p x)"

lemma fst_divide_out [simp]: "fst (divide_out p x) = x div p ^ multiplicity p x"
  and snd_divide_out [simp]: "snd (divide_out p x) = multiplicity p x"
  by (simp_all add: divide_out_def)


function divide_out_aux :: "'a :: factorial_semiring  'a × nat  'a × nat" where
  "divide_out_aux p (x, acc) =
     (if x = 0  is_unit p  ¬p dvd x then (x, acc) else divide_out_aux p (x div p, acc + 1))"
  by auto
termination proof (relation "measure (λ(p, x, _). multiplicity p x)")
  fix p x :: 'a and acc :: nat
  assume "¬(x = 0  is_unit p  ¬p dvd x)"
  thus "((p, x div p, acc + 1), p, x, acc)  measure (λ(p, x, _). multiplicity p x)"
    by (auto elim!: dvdE simp: multiplicity_times_same)
qed auto

lemmas [simp del] = divide_out_aux.simps

lemma divide_out_aux_correct:
  "divide_out_aux p z = (fst z div p ^ multiplicity p (fst z), snd z + multiplicity p (fst z))"
proof (induction p z rule: divide_out_aux.induct)
  case (1 p x acc)
  show ?case
  proof (cases "x = 0  is_unit p  ¬p dvd x")
    case False
    have "x div p div p ^ multiplicity p (x div p) = x div p ^ multiplicity p x"
      using False
      by (subst dvd_div_mult2_eq [symmetric])
         (auto elim!: dvdE simp: multiplicity_dvd multiplicity_times_same)
    with False show ?thesis using 1
      by (subst divide_out_aux.simps)
         (auto elim: dvdE simp: multiplicity_times_same multiplicity_unit_left
                                not_dvd_imp_multiplicity_0)
  qed (auto simp: divide_out_aux.simps multiplicity_unit_left not_dvd_imp_multiplicity_0)
qed

lemma divide_out_code [code]: "divide_out p x = divide_out_aux p (x, 0)"
  by (simp add: divide_out_aux_correct divide_out_def)

lemma multiplicity_code [code]: "multiplicity p x = snd (divide_out_aux p (x, 0))"
  by (simp add: divide_out_aux_correct)

(* TODO Move *)
lemma multiplicity_times_same_power:
  assumes "x  0" "¬is_unit p" "p  0"
  shows   "multiplicity p (p ^ k * x) = multiplicity p x + k"
  using assms by (induction k) (simp_all add: multiplicity_times_same mult.assoc)

lemma divide_out_unique_nat:
  fixes n :: nat
  assumes "¬is_unit p" "p  0" "¬p dvd m" and "n = p ^ k * m"
  shows   "m = n div p ^ multiplicity p n" and "k = multiplicity p n"
proof -
  from assms have "n  0" by (intro notI) auto
  thus *: "k = multiplicity p n"
    using assms by (auto simp: assms multiplicity_times_same_power not_dvd_imp_multiplicity_0)
  from assms show "m = n div p ^ multiplicity p n"
    unfolding * [symmetric] by auto
qed


context
  fixes a n :: nat 
begin

definition "fermat_liar  [a^(n-1) = 1] (mod n)"
definition "fermat_witness  [a^(n-1)  1] (mod n)"

definition "strong_fermat_liar 
              (k m. odd m  n - 1 = 2^k * m 
                 [a^m = 1](mod n)  (i  {0..< k}. [a^(2 ^ i * m) = n-1] (mod n)))"

definition "strong_fermat_witness  ¬ strong_fermat_liar"

lemma fermat_liar_witness_of_composition[iff]:
  "¬fermat_liar = fermat_witness"
  "¬fermat_witness = fermat_liar"
  unfolding fermat_liar_def and fermat_witness_def
  by simp_all

lemma strong_fermat_liar_code [code]:
  "strong_fermat_liar  
     (let (m, k) = divide_out 2 (n - 1)
      in [a^m = 1](mod n)  (i  {0..< k}. [a^(2 ^ i * m) = n-1] (mod n)))"
  (is "?lhs = ?rhs")
proof (cases "n > 1")
  case True
  define m where "m = (n - 1) div 2 ^ multiplicity 2 (n - 1)"
  define k where "k = multiplicity 2 (n - 1)"
  have mk: "odd m  n - 1 = 2 ^ k * m"
    using True multiplicity_decompose[of "n - 1" 2] multiplicity_dvd[of 2 "n - 1"]
    by (auto simp: m_def k_def)
  show ?thesis
  proof
    assume ?lhs
    hence "[a^m = 1] (mod n)  (i  {0..< k}. [a^(2 ^ i * m) = n-1] (mod n))"
      using True mk by (auto simp: divide_out_def strong_fermat_liar_def)
    thus ?rhs by (auto simp: Let_def divide_out_def m_def k_def)
  next
    assume ?rhs
    thus ?lhs using divide_out_unique_nat[of 2]
      by (auto simp: strong_fermat_liar_def divide_out_def)
  qed
qed (auto simp: strong_fermat_liar_def divide_out_def)


context
  assumes * : "a  {1 ..< n}"
begin

lemma strong_fermat_witness_iff:
  "strong_fermat_witness =
     (k m. odd m  n - 1 = 2 ^ k * m  [a ^ m  1] (mod n)  
            (i{0..<k}. [a ^ (2 ^ i * m)  n-1] (mod n)))"
  unfolding strong_fermat_witness_def strong_fermat_liar_def
  by blast

lemma not_coprime_imp_witness: "¬coprime a n  fermat_witness"
  using * lucas_coprime_lemma[of "n-1" a n]
  by (auto simp: fermat_witness_def)

corollary liar_imp_coprime: "fermat_liar  coprime a n"
  using not_coprime_imp_witness fermat_liar_witness_of_composition by blast

lemma fermat_witness_imp_strong_fermat_witness:
  assumes "1 < n" "fermat_witness"
  shows "strong_fermat_witness"
proof -
  from ‹fermat_witness› have "[a^(n-1)  1] (mod n)"
    unfolding fermat_witness_def .

  obtain k m where "odd m" and n: "n - 1 = 2 ^ k * m"
    using * by (auto intro: multiplicity_decompose'[of "(n-1)" 2])

  moreover have "[a ^ m  1] (mod n)"
    using [a^(n-1)  1] (mod n) n ord_divides by auto

  moreover {
    fix i :: nat
    assume "i  {0..<k}"
    hence "i  k - 1" "0 < k" by auto
    then have "[a ^ (2 ^ i * m)  n - 1] (mod n)" "[a ^ (2 ^ i * m)  1] (mod n)"
    proof(induction i rule: inc_induct)
      case base
        have * :"a ^ (n - 1) = (a ^ (2 ^ (k - 1) * m))2"
          using 0 < k n semiring_normalization_rules(27)[of "2 :: nat" "k - 1"]
          by (auto simp flip: power_even_eq simp add: algebra_simps)
  
        case 1
        from * show ?case 
          using [a^(n-1)  1] (mod n) and square_minus_one_cong_one[OF 1 < n] by auto
       
        case 2
        from * show ?case using n [a^(n-1)  1] (mod n) and square_one_cong_one by metis
    next
      case (step q)
      then have
        IH2: "[a ^ (2 ^ Suc q * m)  1] (mod n)" using 0 < k by blast+

      have * : "a ^ (2 ^ Suc q * m) = (a ^ (2 ^ q * m))2"
        by (auto simp flip: power_even_eq simp add: algebra_simps)

      case 1
      from * show ?case using IH2 and square_minus_one_cong_one[OF 1 < n] by auto

      case 2
      from * show ?case using IH2 and square_one_cong_one by metis
    qed
  }

  ultimately show ?thesis unfolding strong_fermat_witness_iff by blast
qed

corollary strong_fermat_liar_imp_fermat_liar:
  assumes "1 < n" "strong_fermat_liar"
  shows  "fermat_liar" 
    using fermat_witness_imp_strong_fermat_witness assms
    and fermat_liar_witness_of_composition strong_fermat_witness_def
    by blast

lemma euler_liar_is_fermat_liar:
  assumes "1 < n" "euler_liar a n" "coprime a n" "odd n"
  shows "fermat_liar"
proof -
  have "[Jacobi a n = a ^ ((n - 1) div 2)] (mod n)"
    using ‹euler_liar a n unfolding euler_witness_def by simp

  hence "[(Jacobi a n)^2 = (a ^ ((n - 1) div 2))^2] (mod n)"
    by (simp add: cong_pow)

  moreover have "Jacobi a n  {1, -1}"
    using Jacobi_values Jacobi_eq_0_iff_not_coprime[of n] ‹coprime a n 1 < n
    by force

  ultimately have "[1 = (a ^ ((n - 1) div 2))^2] (mod n)"
    using cong_int_iff by force

  also have "(a ^ ((n - 1) div 2))^2 = a ^ (n - 1)"
    using ‹odd n by (simp flip: power_mult)

  finally show ?thesis 
    using cong_sym fermat_liar_def 
    by blast
qed

corollary fermat_witness_is_euler_witness:
  assumes "1 < n" "fermat_witness" "coprime a n" "odd n"
  shows "euler_witness a n"
  using assms euler_liar_is_fermat_liar fermat_liar_witness_of_composition
  by blast

end
end

lemma one_is_fermat_liar[simp]: "1 < n  fermat_liar 1 n"
  using fermat_liar_def by auto

lemma one_is_strong_fermat_liar[simp]: "1 < n  strong_fermat_liar 1 n"
  using strong_fermat_liar_def by auto

lemma prime_imp_fermat_liar:
  "prime p  a  {1 ..< p}  fermat_liar a p"
  unfolding fermat_liar_def
  using fermat_theorem and nat_dvd_not_less by simp

lemma not_Carmichael_numberD:
  assumes "¬Carmichael_number n" "¬prime n"  and "1 < n"
  shows " a  {2 ..< n} . fermat_witness a n  coprime a n"
proof -
  obtain a where "coprime a n" "[a^(n-1)  1] (mod n)"
    using assms unfolding Carmichael_number_def by blast

  moreover from this and 1 < n have "a mod n  {1..<n}"
    by (cases "a = 0") (auto intro! : gre1I_nat)

  ultimately have "a mod n  {1 ..< n}"  "coprime (a mod n) n" "[(a mod n)^(n-1)  1] (mod n)"
    using 1 < n by simp_all

  hence "fermat_witness (a mod n) n"
    using fermat_witness_def by blast

  hence "1  (a mod n)"
    using 1 < n (a mod n)  {1 ..< n} and one_is_fermat_liar fermat_liar_witness_of_composition(1)
    by metis

  thus ?thesis
    using ‹fermat_witness (a mod n) n ‹coprime (a mod n) n (a mod n)  {1..<n}
    by fastforce
qed

proposition prime_imp_strong_fermat_witness:
  fixes p :: nat
  assumes "prime p" "2 < p" "a  {1 ..< p}"
  shows "strong_fermat_liar a p"
proof -
  { fix k m :: nat
    define j where "j  LEAST k. [a ^ (2^k * m) = 1] (mod p)"

    have "[a ^ (p - 1) = 1] (mod p)"
      using fermat_theorem[OF ‹prime p, of a] a  {1 ..< p} by force

    moreover assume "odd m" and n: "p - 1 = 2 ^ k * m"

    ultimately have "[a ^ (2 ^ k * m) = 1] (mod p)" by simp

    moreover assume "[a ^ m  1] (mod p)"
    then have "0 < x" if "[a ^ (2 ^ x * m) = 1] (mod p)" for x
      using that by (auto intro: gr0I)

    ultimately have "0 < j" "j  k" "[a ^ (2 ^ j * m) = 1] (mod p)"
      using LeastI2[of _ k "(<) 0"] Least_le[of _ k] LeastI[of _ k]
      by (simp_all add: j_def)

    moreover from this have "[a ^ (2^(j-1) * m)  1] (mod p)"
      using not_less_Least[of "j - 1" "λk. [a ^ (2^k * m) = 1] (mod p)"]
      unfolding j_def by simp

    moreover have "a ^ (2 ^ (j - 1) * m) * a ^ (2 ^ (j - 1) * m) = a ^ (2 ^ j * m)"
      using 0 < j
      by (simp add: algebra_simps semiring_normalization_rules(27) flip: power2_eq_square power_even_eq)

    ultimately have "(j-1){0..<k} " "[a ^ (2 ^ (j-1) * m) = p - 1] (mod p)"
      using cong_square_alt[OF ‹prime p, of "a ^ (2 ^ (j-1) * m)"]
      by simp_all
  }

  then show ?thesis unfolding strong_fermat_liar_def by blast
qed

lemma ignore_one:
  fixes P :: "_  nat  bool"
  assumes "P 1 n" "1 < n"
  shows "card {a  {1..<n}. P a n} = 1 + card {a. 2  a  a < n  P a n}"
proof -
  have "insert 1 {a. 2  a  a < n  P a n} = {a. 1  a  a < n  P a n}"
    using assms by auto

  moreover have "card (insert 1 {a. 2  a  a < n  P a n}) = Suc (card {a. 2  a  a < n  P a n})"
    using card_insert_disjoint by auto
  
  ultimately show ?thesis by force
qed

text ‹Proofs in the following section are inspired by \cite{Cornwell, MillerRabinTest, lecture_notes}.›

proposition not_Carmichael_number_imp_card_fermat_witness_bound:
  fixes n :: nat
  assumes "¬prime n" "¬Carmichael_number n" "odd n" "1 < n"
  shows "(n-1) div 2 < card {a  {1 ..< n}. fermat_witness a n}"
    and "(card {a. 2  a  a < n  strong_fermat_liar a n}) < real (n - 2) / 2"
    and "(card {a. 2  a  a < n  fermat_liar a n}) < real (n - 2) / 2"
proof -
  define G where "G = Residues_Mult n"
  interpret residues_mult_nat n G
    by unfold_locales (use n > 1 in simp_all only: G_def›)
  define h where "h  λa. a ^ (n - 1) mod n"
  define ker where "ker = kernel G (Gcarrier := h ` carrier G) h"

  have "finite ker" by (auto simp: ker_def kernel_def)
  moreover have "1  ker" using n > 1 by (auto simp: ker_def kernel_def h_def)
  ultimately have [simp]: "card ker > 0"
    by (subst card_gt_0_iff) auto

  have totatives_eq: "totatives n = {k{1..<n}. coprime k n}"
    using totatives_less[of _ n] n > 1 by (force simp: totatives_def)
  have ker_altdef: "ker = {a  {1..<n}. fermat_liar a n}"
    unfolding ker_def fermat_liar_def carrier_eq kernel_def totatives_eq using n > 1
    by (force simp: h_def cong_def intro: coprimeI_power_mod)

  have h_is_hom: "h  hom G G" 
    unfolding hom_def using nat_pow_closed
    by (auto simp: h_def power_mult_distrib mod_simps)
  then interpret h: group_hom G G h
    by unfold_locales

  obtain a where a: "a  {2..<n}" "fermat_witness a n" "coprime a n" 
    using assms power_one not_Carmichael_numberD by blast

  have "h a  1" using a by (auto simp: fermat_witness_def cong_def h_def)
  hence "2  card {1, h a}" by simp
  also have "  card (h ` carrier G)"
  proof (intro card_mono; safe?)
    from n > 1 have "1 = h 1" by (simp add: h_def)
    also have "  h ` carrier G" by (intro imageI) (use n > 1 in auto)
    finally show "1  h ` carrier G" .
  next
    show "h a  h ` carrier G"
      using a by (intro imageI) (auto simp: totatives_def)
  qed auto
  also have " * card ker = order G"
    using homomorphism_thm_order[OF h.group_hom_axioms] by (simp add: ker_def order_def)
  also have "order G < n - 1"
    using totient_less_not_prime[of n] assms by (simp add: order_eq)
  finally have "card ker < (n - 1) div 2"
    using ‹odd n by (auto elim!: oddE)

  have "(n - 1) div 2 < (n - 1) - card ker"
    using ‹card ker < (n - 1) div 2 by linarith
  also have " = card ({1..<n} - ker)"
    by (subst card_Diff_subset) (auto simp: ker_altdef)
  also have "{1..<n} - ker = {a  {1..<n}. fermat_witness a n}"
    by (auto simp: fermat_witness_def fermat_liar_def ker_altdef)
  finally show "(n - 1) div 2 < card {a  {1..<n}. fermat_witness a n}" .

  have "card {a. 2  a  a < n  fermat_liar a n}  card (ker - {1})"
    by (intro card_mono) (auto simp: ker_altdef fermat_liar_def fermat_witness_def)
  also have " = card ker - 1"
    using n > 1 by (subst card_Diff_subset) (auto simp: ker_altdef fermat_liar_def)
  also have " < (n - 2) div 2"
    using ‹card ker < (n - 1) div 2 ‹odd n ‹card ker > 0 by linarith
  finally show *: "card {a. 2  a  a < n  fermat_liar a n} < real (n - 2) / 2"
    by simp

  have "card {a. 2  a  a < n  strong_fermat_liar a n} 
          card {a. 2  a  a < n  fermat_liar a n}"
    by (intro card_mono) (auto intro!: strong_fermat_liar_imp_fermat_liar)
  moreover note *
  ultimately show "card {a. 2  a  a < n  strong_fermat_liar a n} < real (n - 2) / 2"
    by simp
qed


proposition Carmichael_number_imp_lower_bound_on_strong_fermat_witness:
  fixes n :: nat
  assumes Carmichael_number: "Carmichael_number n"
  shows "(n - 1) div 2 < card {a  {1..<n}. strong_fermat_witness a n}"
    and "real (card {a . 2  a  a < n  strong_fermat_liar a n}) < real (n - 2) / 2"
proof -
  from assms have "n > 3" by (intro Carmichael_number_gt_3)
  hence "n - 1  0" "¬is_unit (2 :: nat)" by auto
  obtain k m where "odd m" and n_less: "n - 1 = 2 ^ k * m"
    using multiplicity_decompose'[OF n - 1  0 ¬is_unit (2::nat)] by metis

  obtain p l where n: "n = p * l" and "prime p" "¬ p dvd l" "2 < l"
    using Carmichael_number_imp_squarefree_alt[OF Carmichael_number] 
    by blast

  then have "coprime p l" using prime_imp_coprime_nat by blast

  have "odd n" using Carmichael_number_odd Carmichael_number by simp

  have "2 < n" using n > 3 ‹odd n by presburger

  note prime_gt_1_nat[OF ‹prime p]

  have "2 < p" using ‹odd n n ‹prime p prime_ge_2_nat
               and dvd_triv_left le_neq_implies_less by blast 

  let ?P = "λ k. ( a. coprime a p  [a^(2^k * m) = 1] (mod p))"
  define j where "j  LEAST k. ?P k"

  define H where "H  {a  {1..<n} . coprime a n  ([a^(2^(j-1) * m) = 1] (mod n) 
                                                     [a^(2^(j-1) * m) = n - 1] (mod n))}"

  have k : "a. coprime a n  [a ^ (2 ^ k * m) = 1] (mod n)"
    using Carmichael_number unfolding Carmichael_number_def n_less by blast

  obtain k' m' where "odd m'" and p_less: "p - 1 = 2 ^ k' * m'"
    using 1 < p by (auto intro: multiplicity_decompose'[of "(p-1)" 2])

  have "p - 1 dvd n - 1"
    using Carmichael_number_imp_dvd[OF Carmichael_number ‹prime p] n = p * l
    by fastforce

  then have "p - 1 dvd 2 ^ k' * m"
    unfolding n_less p_less
    using ‹odd m ‹odd m'
      and coprime_dvd_mult_left_iff[of "2^k'" m "2^k"] coprime_dvd_mult_right_iff[of m' "2^k" m]
    by auto

  have k': "a. coprime a p  [a ^ (2 ^ k' * m) = 1] (mod p)"
  proof safe
    fix a
    assume "coprime a p"
    hence "¬ p dvd a" using p_coprime_right_nat[OF ‹prime p] by simp

    have "[a ^ (2 ^ k' * m') = 1] (mod p)"
     unfolding p_less[symmetric]
      using fermat_theorem ‹prime p ¬ p dvd a by blast

    then show "[a ^ (2 ^ k' * m) = 1] (mod p)"
      using p - 1 dvd 2 ^ k' * m
      unfolding p_less n_less
      by (meson dvd_trans ord_divides)
  qed

  have j_prop: "[a^(2^j * m) = 1] (mod p)" if "coprime a p" for a
    using that LeastI[of ?P k', OF k', folded j_def] cong_modulus_mult coprime_mult_right_iff 
    unfolding j_def n by blast

  have j_least: "[a^(2^i * m) = 1] (mod p)" if "coprime a p" "j  i" for a i
  proof -
    obtain c where i: "i = j + c" using le_iff_add[of j i] j  i by blast

    then have "[a ^ (2 ^ i * m) = a ^ (2 ^ (j + c) * m)] (mod p)" by simp
    
    also have "[a ^ (2 ^ (j + c) * m) = (a ^ (2 ^ j * m)) ^ (2 ^ c)] (mod p)"
      by (simp flip: power_mult add: algebra_simps power_add)
    
    also note j_prop[OF ‹coprime a p]
   
    also have "[1 ^ (2 ^ c) = 1] (mod p)" by simp
    
    finally show ?thesis .
  qed

  have neq_p: "[p - 1  1](mod p)"
    using 2 < p and cong_less_modulus_unique_nat[of "p-1" 1 p]
    by linarith

  have "0 < j"
  proof (rule LeastI2[of ?P k', OF k', folded j_def], rule gr0I)
    fix x
    assume "a. coprime a p  [a ^ (2 ^ x * m) = 1] (mod p)"
    then have "[(p - 1) ^ (2 ^ x * m) = 1] (mod p)"
      using coprime_diff_one_left_nat[of p]  prime_gt_1_nat[OF ‹prime p]
      by simp
    
    moreover assume "x = 0"
    hence "[(p-1)^(2^x*m) = p - 1](mod p)"
      using ‹odd m odd_pow_cong[OF _ ‹odd m, of p] prime_gt_1_nat[OF ‹prime p]
      by auto

    ultimately show False
      using [p - 1  1] (mod p) by (simp add: cong_def)
  qed

  then have "j - 1 < j" by simp

  then obtain x where "coprime x p" "[x^(2^(j-1) * m)  1] (mod p)"
    using not_less_Least[of "j - 1" ?P, folded j_def] unfolding j_def by blast

  define G where "G = Residues_Mult n"
  interpret residues_mult_nat n G
    by unfold_locales (use n > 3 in simp_all only: G_def›)

  have H_subset: "H  carrier G" unfolding H_def by (auto simp: totatives_def)

  from n > 3 have n > 1 by simp
  interpret H: subgroup H G
  proof (rule subgroupI, goal_cases)
    case 1
    then show ?case using H_subset .
  next
    case 2
    then show ?case unfolding H_def using 1 < n by force
  next
    case (3 a)
    define y where "y = invG a"

    then have "y  carrier G"
      using H_subset a  H by (auto simp del: carrier_eq)

    then have "1  y" "y < n" "coprime y n"
      using totatives_less[of y n] n > 3 by (auto simp: totatives_def)

    moreover have "[y ^ (2 ^ (j - Suc 0) * m) = Suc 0] (mod n)" 
      if "[y ^ (2 ^ (j - Suc 0) * m)  n - Suc 0] (mod n)"
    proof -
      from a  H have "[a * y = 1] (mod n)"
        using H_subset r_inv[of a] y_def by (auto simp: cong_def)
      hence "[(a * y) ^ (2 ^ (j - 1) * m) = 1 ^ (2 ^ (j - 1) * m)] (mod n)"
        by (intro cong_pow)
      hence "[(a * y) ^ (2 ^ (j - 1) * m) = 1] (mod n)"
        by simp
      hence * : "[a ^ (2 ^ (j - 1) * m) * y ^ (2 ^ (j - 1) * m) = 1] (mod n)"
          by (simp add: power_mult_distrib) 
      from a  H have "1  a" "a < n" "coprime a n"
        unfolding H_def by auto

      have "[a ^ (2 ^ (j - 1) * m) = 1] (mod n)  [a ^ (2 ^ (j - 1) * m) = n - 1] (mod n)"
        using a  H by (auto simp: H_def)
      thus ?thesis
      proof
        note *
        also assume "[a ^ (2 ^ (j - 1) * m) = 1] (mod n)"
        finally show ?thesis by simp
      next
        assume "[a ^ (2 ^ (j - 1) * m) = n - 1] (mod n)"
        then have "[y ^ (2 ^ (j - 1) * m) = n - 1] (mod n)" 
          using minus_one_cong_solve[OF 1 < n] * ‹coprime a n ‹coprime y ncoprime_power_left_iff
          by blast+
        thus ?thesis using that by simp
      qed
    qed

    ultimately show ?case using a  H unfolding H_def y_def by auto
  next
    case (4 a b)
    hence "a  totatives n" "b  totatives n"
      by (auto simp: H_def totatives_def)
    hence "a * b mod n  totatives n"
      using m_closed[of a b] by simp
    hence "a * b mod n  {1..<n}" "coprime (a * b) n"
      using totatives_less[of "a * b" n] n > 3 by (auto simp: totatives_def)

    moreover define x y where "x = a ^ (2 ^ (j - 1) * m)" and "y = b ^ (2 ^ (j - 1) * m)"
    have "[x * y = 1] (mod n)  [x * y = n - 1] (mod n)"
    proof -
      have *: "x mod n  {1, n - 1}" "y mod n  {1, n - 1}"
        using 4 by (auto simp: H_def x_def y_def cong_def)
      have "[x * y = (x mod n) * (y mod n)] (mod n)"
        by (intro cong_mult) auto
      moreover have "((x mod n) * (y mod n)) mod n  {1, n - 1}"
        using * square_minus_one_cong_one'[OF 1 < n] n > 1 by (auto simp: cong_def)
      ultimately show ?thesis using n > 1 by (simp add: cong_def mod_simps)
    qed

    ultimately show ?case by (auto simp: H_def x_def y_def power_mult_distrib)
  qed

  { obtain a where "[a = x] (mod p)" "[a = 1] (mod l)" "a < p * l"
      using binary_chinese_remainder_unique_nat[of p l x 1]
        and ¬ p dvd l ‹prime p prime_imp_coprime_nat
      by auto
  
    moreover have "coprime a p"
      using ‹coprime x p cong_imp_coprime[OF cong_sym[OF [a = x] (mod p)]] coprime_mult_right_iff
      unfolding n by blast
  
    moreover have "coprime a l" 
      using coprime_1_left cong_imp_coprime[OF cong_sym[OF [a = 1] (mod l)]]
      by blast

    moreover from ‹prime p and ‹coprime a p have "a > 0"
      by (intro Nat.gr0I) auto
 
    ultimately have "a  carrier G"
      using 2 < l by (auto intro: gre1I_nat simp: n totatives_def)
  
    have "[a ^ (2^(j-1) * m)  1] (mod p)"
      using [x^(2^(j-1) * m)  1] (mod p) [a = x] (mod p) and cong_trans cong_pow cong_sym
      by blast

    then have "[a ^ (2^(j-1) * m)  1] (mod n)"
      using cong_modulus_mult_nat n by fast

    moreover 
    have "[a ^ (2 ^ (j - Suc 0) * m)  n - 1] (mod n)"
    proof -
      have "[a ^ (2 ^ (j - 1) * m) = 1] (mod l)"
      using cong_pow[OF [a = 1] (mod l)] by auto
  
      moreover have "Suc 0  (n - Suc 0) mod l"
        using n 2 < l ‹odd n
        by (metis mod_Suc_eq mod_less mod_mult_self2_is_0 numeral_2_eq_2 odd_Suc_minus_one zero_neq_numeral)

      then have "[1  n - 1] (mod l)"
          using 2 < l ‹odd n unfolding cong_def by simp

      moreover have "l  Suc 0" using 2 < l by simp

      ultimately have "[a ^ (2 ^ (j - Suc 0) * m)  n - 1] (mod l)"
        by (auto simp add: cong_def n mod_simps dest: cong_modulus_mult_nat)

      then show ?thesis
        using cong_modulus_mult_nat mult.commute n by metis
    qed

    ultimately have "a  H" unfolding H_def by auto
  
    hence "H  carrier (G)"
      using H_subset subgroup.mem_carrier and a  carrier (G) 
      by fast
  }

  have "card H  order G div 2"
    by (intro proper_subgroup_imp_bound_on_card) (use H  carrier G H.is_subgroup in auto)
  also from assms have "¬prime n" by (auto dest: Carmichael_number_not_prime)
  hence "order G div 2 < (n - 1) div 2"
    using totient_less_not_prime[OF ¬ prime n 1 < n] ‹odd n
    by (auto simp add: order_eq elim!: oddE)
  finally have "card H < (n - 1) div 2" .

  {
    { fix a
      assume "1  a" "a < n"
      hence "a  {1..<n}" by simp
  
      assume "coprime a n"
      then have "coprime a p" 
      unfolding n by simp
                 
      assume "[a ^ (2 ^ (j - 1) * m)  1] (mod n)" 
      hence "[a ^ m  1] (mod n)"
        by (metis dvd_trans dvd_triv_right ord_divides)
    
      moreover assume "strong_fermat_liar a n"
  
      ultimately obtain i where "i  {0 ..< k}" "[a^(2^i * m) = n-1](mod n)"
        unfolding strong_fermat_liar_def using ‹odd m n_less by blast
  
      then have "[a ^ (2 ^ i * m) = n - 1] (mod p)"
        unfolding n using cong_modulus_mult_nat by blast
  
      moreover have "[n - 1  1] (mod p)"
      proof(subst cong_altdef_nat, goal_cases)
        case 1
        then show ?case using 1 < n by linarith
      next
        case 2
        have "¬ p dvd 2" using 2 < p by (simp add: nat_dvd_not_less)

        moreover have "2  n" using 1 < n by linarith

        moreover have "p dvd n" using n by simp

        ultimately have "¬ p dvd n - 2" using dvd_diffD1 by blast
        
        then show ?case by (simp add: numeral_2_eq_2)
      qed

      ultimately have "[a ^ (2 ^ i * m)  Suc 0] (mod p)" using cong_sym by (simp add: cong_def)
  
      then have "i < j" using j_least[OF ‹coprime a p, of i] by force
  
      have "[(a ^ (2 ^ Suc i * m)) = 1] (mod n)"
        using square_minus_one_cong_one[OF 1 < n [a^(2^i * m) = n-1](mod n)]
        by (simp add: power2_eq_square power_mult power_mult_distrib)
  
      { assume "i < j - 1"
  
        have "(2 :: nat) ^ (j - Suc 0) = ((2 ^ i) * 2 ^ (j - Suc i))"
          unfolding power_add[symmetric] using i < j - 1 by simp
  
        then have "[a ^ (2 ^ (j - 1) * m) = (a ^ (2 ^ i * m)) ^ (2^(j - 1 - i))] (mod n)"
          by (auto intro!: cong_pow_I simp flip: power_mult simp add: algebra_simps power_add)
          
        also note [a ^ (2 ^ i * m) = n - 1] (mod n)
  
        also have "[(n - 1) ^ (2^(j - 1 - i)) = 1] (mod n) "
          using 1 < n i < j - 1 using even_pow_cong by auto
  
        finally have False
          using [a ^ (2 ^ (j - 1) * m)  1] (mod n)
          by blast
      }

      hence "i = j - 1" using i < j by fastforce
  
      hence "[a ^ (2 ^ (j - 1) * m) = n - 1] (mod n)" using [a^(2^i * m) = n-1](mod n) by simp
    }
  
    hence "{a  {1..<n}. strong_fermat_liar a n}  H"
      using strong_fermat_liar_imp_fermat_liar[of _ n, OF _ 1 < n]  liar_imp_coprime
      by (auto simp: H_def)
  }

  moreover have "finite H" unfolding H_def by auto

  ultimately have strong_fermat_liar_bounded: "card {a  {1..<n}. strong_fermat_liar a n} < (n - 1) div 2 "
    using card_mono[of H] le_less_trans[OF _ ‹card H < (n - 1) div 2] by blast

  moreover {
    have "{1..<n} - {a  {1..<n}. strong_fermat_liar a n} = {a  {1..<n}. strong_fermat_witness a n}"
      using strong_fermat_witness_def by blast
  
    then have "card {a  {1..<n}. strong_fermat_witness a n} = (n-1) - card {a  {1..<n}. strong_fermat_liar a n}"
      using card_Diff_subset[of "{a  {1..<n}. strong_fermat_liar a n}" "{1..<n}"]
      by fastforce
  }
  
  ultimately show "(n - 1) div 2 < card {a  {1..<n}. strong_fermat_witness a n}"
    by linarith

  show "real (card {a . 2  a  a < n  strong_fermat_liar a n}) < real (n - 2) / 2"
    using strong_fermat_liar_bounded ignore_one one_is_strong_fermat_liar 1 < n 
    by simp
qed

corollary strong_fermat_witness_lower_bound:
  assumes "odd n" "n > 2" "¬prime n"
  shows   "card {a. 2  a  a < n  strong_fermat_liar a n} < real (n - 2) / 2"
  using Carmichael_number_imp_lower_bound_on_strong_fermat_witness(2)[of n]
        not_Carmichael_number_imp_card_fermat_witness_bound(2)[of n] assms
  by (cases "Carmichael_number n") auto

end

Theory Generalized_Primality_Test

(*
  File:     Generalized_Primality_Test.thy
  Authors:  Daniel Stüwe, Manuel Eberl

  Generic probabilistic primality test
*)
section ‹A Generic View on Probabilistic Prime Tests›
theory Generalized_Primality_Test
imports 
  "HOL-Probability.Probability" 
  Algebraic_Auxiliaries
begin

definition primality_test :: "(nat  nat  bool)  nat  bool pmf" where
  "primality_test P n = 
    (if n < 3  even n then return_pmf (n = 2) else
      do {
        a  pmf_of_set {2..< n};
        return_pmf (P n a)
      })"

(* TODO Move *)
lemma expectation_of_bool_is_pmf: "measure_pmf.expectation M of_bool = pmf M True"
  by (simp add: integral_measure_pmf_real[where A=UNIV] UNIV_bool)

lemma eq_bernoulli_pmfI:
  assumes "pmf p True = x"
  shows   "p = bernoulli_pmf x"
proof (intro pmf_eqI)
  fix b :: bool
  from assms have "x  {0..1}" by (auto simp: pmf_le_1)
  thus "pmf p b = pmf (bernoulli_pmf x) b"
    using assms by (cases b) (auto simp: pmf_False_conv_True)
qed

text ‹
  We require a probabilistic primality test to never classify a prime as composite.
  It may, however, mistakenly classify composites as primes.
›
locale prob_primality_test =
  fixes P :: "nat  nat  bool" and n :: nat
  assumes P_works: "odd n  2  a  a < n  prime n  P n a"
begin

lemma FalseD:
  assumes false: "False  set_pmf (primality_test P n)"
  shows   "¬ prime n"
proof -
  from false consider "n  2" "n < 3" | "n  2" "even n" | 
      a where "¬ P n a" "odd n" "2  a" "a < n"
    by (auto simp: primality_test_def not_less split: if_splits)

  then show ?thesis proof cases
  case 1
  then show ?thesis 
    by (cases rule: linorder_neqE_nat) (use prime_ge_2_nat[of n] in auto)
  next
    case 2
    then show ?thesis using primes_dvd_imp_eq two_is_prime_nat by blast
  next
    case 3
    then show ?thesis using P_works by blast
  qed 
qed

theorem prime:
  assumes odd_prime: "prime n"
  shows   "primality_test P n = return_pmf True"
proof -
  have "set_pmf (primality_test P n)  {True, False}"
    by auto
  moreover from assms have "False  set_pmf (primality_test P n)"
    using FalseD by auto
  ultimately have "set_pmf (primality_test P n)  {True}"
    by auto
  thus ?thesis
    by (subst (asm) set_pmf_subset_singleton)
qed

end

text ‹
  We call a primality test q›-good for a fixed positive real number q› if the probability
  that it mistakenly classifies a composite as a prime is less than q›.
›
locale good_prob_primality_test = prob_primality_test +
  fixes q :: real
  assumes q_pos: "q > 0"
  assumes composite_witness_bound:
            "¬prime n  2 < n  odd n 
               real (card {a . 2  a  a < n  P n a}) < q * real (n - 2)"
begin

lemma composite_aux:
  assumes "¬prime n"
  shows   "measure_pmf.expectation (primality_test P n) of_bool < q"
  unfolding primality_test_def using assms composite_witness_bound q_pos
  by (clarsimp simp add: pmf_expectation_bind[where A = "{2..< n}"] sum_of_bool_eq_card field_simps Int_def
    simp flip: sum_divide_distrib)

theorem composite:
  assumes "¬prime n"
  shows   "pmf (primality_test P n) True < q"
  using composite_aux[OF assms] by (simp add: expectation_of_bool_is_pmf)

end

end

Theory Fermat_Test

(*
  File:     Fermat_Test.thy
  Authors:  Daniel Stüwe, Manuel Eberl

  The probabilistic prime test known as Fermat's test
*)
section ‹Fermat's Test›
theory Fermat_Test
imports 
  Fermat_Witness
  Generalized_Primality_Test
begin

definition "fermat_test = primality_test (λn a. fermat_liar a n)"

text ‹
  The Fermat test is a good probabilistic primality test on non-Carmichael numbers.
›
locale fermat_test_not_Carmichael_number =
  fixes n :: nat
  assumes not_Carmichael_number: "¬Carmichael_number n  n < 3"
begin

sublocale fermat_test: good_prob_primality_test "λa n. fermat_liar n a" n "1 / 2"
  rewrites "primality_test (λ a n. fermat_liar n a) = fermat_test"
proof -
  show "good_prob_primality_test (λa n. fermat_liar n a) n (1 / 2)"
    using not_Carmichael_number not_Carmichael_number_imp_card_fermat_witness_bound(3)[of n]
          prime_imp_fermat_liar[of n]
    by unfold_locales auto
qed (auto simp: fermat_test_def)

end

lemma not_coprime_imp_fermat_witness:
  fixes n :: nat
  assumes "n > 1" "¬coprime a n"
  shows   "fermat_witness a n"
  using assms lucas_coprime_lemma[of "n - 1" a n]
  by (auto simp: fermat_witness_def)

theorem fermat_test_prime:
  assumes "prime n"
  shows   "fermat_test n = return_pmf True"
proof -
  interpret fermat_test_not_Carmichael_number n
    using assms Carmichael_number_not_prime by unfold_locales auto
  from assms show ?thesis by (rule fermat_test.prime)
qed

theorem fermat_test_composite:
  assumes "¬prime n" "¬Carmichael_number n  n < 3"
  shows   "pmf (fermat_test n) True < 1 / 2"
proof -
  interpret fermat_test_not_Carmichael_number n by unfold_locales fact+
  from assms(1) show ?thesis by (rule fermat_test.composite)
qed

text ‹
  For a Carmichael number $n$, Fermat's test as defined above mistakenly returns `True'
  with probability $(\varphi(n)-1) / (n - 2)$. This probability is close to 1 if n›
  has few and big prime factors; it is not quite as bad if it has many and/or small factors,
  but in that case, simple trial division can also detect compositeness.

  Moreover, Fermat's test only succeeds for a Carmichael number if it happens to guess a
  number that is not coprime to n›. In that case, the fact that we have found
  a number between 2 and n› that is not coprime to n› alone is 
  proof that n› is composite, and indeed we can even find a non-trivial factor
  by computing the GCD. This means that for Carmichael numbers, Fermat's test is essentially
  no better than the very crude method of attempting to guess numbers coprime to n›.

  This means that, in general, Fermat's test is not very helpful for Carmichael numbers.
›
theorem fermat_test_Carmichael_number:
  assumes "Carmichael_number n"
  shows   "fermat_test n = bernoulli_pmf (real (totient n - 1) / real (n - 2))"
proof (rule eq_bernoulli_pmfI)
  from assms have n: "n > 3" "odd n"
    using Carmichael_number_odd Carmichael_number_gt_3 by auto
  from n have "fermat_test n = pmf_of_set {2..<n}  (λa. return_pmf (fermat_liar a n))"
    by (simp add: fermat_test_def primality_test_def)
  also have " = pmf_of_set {2..<n}  (λa. return_pmf (coprime a n))"
    using n assms lucas_coprime_lemma[of "n - 1" _ n]
    by (intro bind_pmf_cong refl) (auto simp: Carmichael_number_def fermat_liar_def)
  also have "pmf  True = (a=2..<n. indicat_real {True} (coprime a n)) / real (n - 2)"
    using n by (auto simp: pmf_bind_pmf_of_set)
  also have "(a=2..<n. indicat_real {True} (coprime a n)) =
               (a | a  {2..<n}  coprime a n. 1)"
    by (intro sum.mono_neutral_cong_right) auto
  also have " = card {a{2..<n}. coprime a n}"
    by simp
  also have "{a{2..<n}. coprime a n} = totatives n - {1}"
    using n by (auto simp: totatives_def order.strict_iff_order[of _ n])
  also have "card  = totient n - 1"
    using n by (subst card_Diff_subset) (auto simp: totient_def)
  finally show "pmf (fermat_test n) True = real (totient n - 1) / real (n - 2)"
    using n by simp
qed

end

Theory Miller_Rabin_Test

(*
  File:     Miller_Rabin.thy
  Authors:  Daniel Stüwe, Manuel Eberl

  The Miller--Rabin primality test.
*)
section ‹The Miller--Rabin Test›
theory Miller_Rabin_Test
imports 
  Fermat_Witness
  Generalized_Primality_Test
begin

definition "miller_rabin = primality_test (λn a. strong_fermat_liar a n)"

text ‹
  The test is actually $\frac{1}{4}$ good, but we only show $\frac{1}{2}$, since the former
  is much more involved.
›
interpretation miller_rabin: good_prob_primality_test "λn a. strong_fermat_liar a n" n "1 / 2"
  rewrites "primality_test (λn a. strong_fermat_liar a n) = miller_rabin"
proof -
  show "good_prob_primality_test (λn a. strong_fermat_liar a n) n (1 / 2)"
    by standard (use strong_fermat_witness_lower_bound prime_imp_strong_fermat_witness in auto)
qed (simp_all add: miller_rabin_def)

end

Theory Solovay_Strassen_Test

(*
  File:     Solovay_Strassen.thy
  Authors:  Daniel Stüwe, Manuel Eberl

  The Solovay--Strassen primality test.
*)
section ‹The Solovay--Strassen Test›
theory Solovay_Strassen_Test
imports 
  Generalized_Primality_Test
  Euler_Witness
begin

definition solovay_strassen_witness :: "nat  nat  bool" where
  "solovay_strassen_witness n a =
     (let x = Jacobi (int a) (int n) in x  0  [x = int a ^ ((n - 1) div 2)] (mod n))"

definition solovay_strassen :: "nat  bool pmf" where
  "solovay_strassen = primality_test solovay_strassen_witness"

lemma prime_imp_solovay_strassen_witness:
  assumes "prime p" "odd p" "a  {2..<p}"
  shows   "solovay_strassen_witness p a"
proof -
  have eq: "Jacobi a p = Legendre a p"
    using prime_p_Jacobi_eq_Legendre assms by simp
  from ‹prime p have "coprime p a"
    by (rule prime_imp_coprime) (use assms in auto)

  show ?thesis unfolding solovay_strassen_witness_def Let_def eq
  proof
    from ‹coprime p a and ‹prime p show "Legendre (int a) (int p)  0"
      by (auto simp: coprime_commute)
  next
    show "[Legendre (int a) (int p) = int a ^ ((p - 1) div 2)] (mod int p)"
      using assms by (intro euler_criterion) auto
  qed
qed
    
lemma card_solovay_strassen_liars_composite:
  fixes n :: nat
  assumes "¬prime n" "n > 2" "odd n"
  shows   "card {a  {2..<n}. solovay_strassen_witness n a} < (n - 2) div 2"
    (is "card ?A < _")
proof -
  interpret euler_witness_context n
    using assms unfolding euler_witness_context_def by simp
  have "card H < (n - 1) div 2"
    by (intro card_euler_liars_cosets_limit(2) assms)
  also from assms have "H = insert 1 ?A"
    by (auto simp: solovay_strassen_witness_def Let_def
                   euler_witness_def H_def Jacobi_eq_0_iff_not_coprime)
  also have "card  = card ?A + 1"
    by (subst card.insert) auto
  finally show "card ?A < (n - 2) div 2"
    by linarith
qed

interpretation solovay_strassen: good_prob_primality_test solovay_strassen_witness n "1 / 2"
  rewrites "primality_test solovay_strassen_witness = solovay_strassen"
proof -
  show "good_prob_primality_test solovay_strassen_witness n (1 / 2)"
  proof
    fix n :: nat assume "¬prime n" "n > 2" "odd n"
    thus "real (card {a. 2  a  a < n  solovay_strassen_witness n a}) < (1 / 2) * real (n - 2)"
      using card_solovay_strassen_liars_composite[of n] by auto
  qed (use prime_imp_solovay_strassen_witness in auto)
qed (simp_all add: solovay_strassen_def)

end