Theory Euler_Witness

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

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

text ‹
Proofs are inspired by \<^cite>‹"solovay_strassen_ori" and "SolovayStrassenTest" and "wiki:Solovay_Strassen_test" and "planetmath:SolovayStrassenTest"›.
›
definition "euler_witness a p ⟷ [Jacobi a p ≠ a ^ ((p - 1) div 2)] (mod p)"
abbreviation "euler_liar a p ≡ ¬ euler_witness a p"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    then show ?thesis 
      using ‹euler_witness x n› by blast

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

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

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

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

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

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

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

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

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

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

    then have "m ≠ 0" using m by fastforce 

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

    hence "?B mod n = 0" by presburger

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

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

    have "[(n div p + 1) ≠ 1] (mod n)"
      using ‹2 < p› ‹odd n› and ‹prime p› ‹p⇧2 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› ‹p⇧2 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› ‹p⇧2 dvd n› ‹p dvd n div p› n
      unfolding a_def
      by (metis coprime_add_one_right coprime_mult_left_iff dvd_def)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    ultimately show ?thesis by linarith
  qed

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

  ultimately show ?thesis by blast
qed

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

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

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

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

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

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

definition G where "G = Residues_Mult p"

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

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

    define y where "y = inv⇘G⇙ 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 ∈ rcosets⇘G⇙ 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} ⊆ rcosets⇘G⇙ 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 (rcosets⇘G⇙ 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 (rcosets⇘G⇙ 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 (rcosets⇘G⇙ H)" "card H < (p - 1) div 2"
proof -
  have "H ∈ rcosets⇘G⇙ 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 (rcosets⇘G⇙ 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 ∈ rcosets⇘G⇙ 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 ∈ rcosets⇘G⇙ 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