# 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
*)
theory Legendre_Symbol
imports

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

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

unfolding cong_def
by presburger

lemma Suc_mod_eq[simp]: "[Suc a = Suc b] (mod 2) = [a = b] (mod 2)"

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

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]

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

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

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)

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

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

also have "[(n - 1) * (n - 1) = Suc (n * n) - (n + n)] (mod n)"
using power2_diff_nat[of 1 n] 1 < n

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

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)

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+

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

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

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

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"

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

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

*)
imports
Jacobi_Symbol
Algebraic_Auxiliaries
begin

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

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

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

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

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

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

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

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

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

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

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}} = {}"
unfolding that
by simp

thus ?thesis by blast
qed

"1 < p  odd p  x. ¬QuadRes p x"
proof (induction p rule: prime_divisors_induct)
case (factor p x)
then show ?case
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
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

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 ‹prime p prime_gt_1_int by auto

then have "[b  0] (mod p)"
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")

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]

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

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

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

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

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"

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

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

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

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]

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

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

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

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

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)