Theory Probabilistic_Prime_Tests.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" and "Jacobi_Symbol".
  It is defined as
  \[\left(\frac{a}{n}\right) =
      \left(\frac{a}{p_1}\right)^{k_1} \ldots \left(\frac{a}{p_l}\right)^{k_l}\]
  where $(\frac{a}{p})$ denotes the Legendre symbol, a› is an integer, n› is an odd natural
  number and $p_1^{k_1}\ldots p_l^{k_l}$ is its prime factorisation.

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

    finally show ?case .
  qed
qed

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

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

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

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

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

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

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

  finally show ?thesis .

qed

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

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

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

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

  then have "odd x" by force

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

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

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

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

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

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

    also note base_case

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

    finally show ?thesis .
  qed
qed simp_all

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

lemma mod_int_wlog [consumes 1, case_names modulo]:
  fixes P :: "int  bool"
  assumes "b > 0"
  assumes "k. 0  k  k < b  n mod b = k  P n"
  shows   "P n"
  using b > 0 assms(2) [of n mod b] by simp

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

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

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

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

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

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

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

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

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

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

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

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

    also note base_case

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

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

    finally show ?thesis .
  qed
qed simp_all

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

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

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

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

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

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

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

lemmas [simp del] = jacobi_code.simps

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

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

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

end