(*
  File:     Gauss_Sums_Auxiliary.thy
  Authors:  Rodrigo Raya, EPFL; Manuel Eberl, TUM
*)
section ‹Auxiliary material›
theory Gauss_Sums_Auxiliary
imports
  Dirichlet_L.Dirichlet_Characters
  Dirichlet_Series.Moebius_Mu
  Dirichlet_Series.More_Totient
begin

subsection ‹Various facts›

lemma sum_div_reduce:
  fixes d :: nat and f :: "nat ⇒ complex"
  assumes "d dvd k" "d > 0" 
  shows "(∑n | n ∈ {1..k} ∧ d dvd n. f n) = (∑c ∈ {1..k div d}. f (c*d))" 
   by (rule sum.reindex_bij_witness[of _ "λk. k * d" "λk. k div d"])
     (use assms in ‹fastforce simp: div_le_mono›)+

lemma prod_div_sub:
  fixes f :: "nat ⇒ complex"
  assumes "finite A" "B ⊆ A" "∀b ∈ B. f b ≠ 0"
  shows "(∏ i ∈ A - B. f i) = ((∏ i ∈ A. f i) div (∏ i ∈ B. f i))"
  using assms
proof (induction "card B" arbitrary: B)
case 0
  then show ?case 
    using infinite_super by fastforce
next
  case (Suc n)
  then show ?case 
  proof -
    obtain B' x where decomp: "B = B' ∪ {x} ∧ x ∉ B'" 
      using card_eq_SucD[OF Suc(2)[symmetric]] insert_is_Un by auto
    then have B'card: "card B' = n" using Suc(2) 
      using Suc.prems(2) assms(1) finite_subset by fastforce
    have "prod f (A - B) = prod f ((A-B') - {x})"
      by (simp add: decomp,subst Diff_insert,simp)
    also have "… = (prod f (A-B')) div f x"  
      using prod_diff1[of "A-B'" f x] Suc decomp by auto
    also have "… = (prod f A div prod f B') div f x"
      using Suc(1)[of B'] Suc(3) B'card decomp
            Suc.prems(2) Suc.prems(3) by force
    also have "… = prod f A div (prod f B' * f x)" by auto
    also have "… = prod f A div prod f B"
      using decomp Suc.prems(2) assms(1) finite_subset by fastforce
    finally show ?thesis by blast
  qed
qed 

lemma linear_gcd:
  fixes a b c d :: nat
  assumes "a > 0" "b > 0" "c > 0" "d > 0"
  assumes "coprime a c" "coprime b d" 
  shows "gcd (a*b) (c*d) = (gcd a d) * (gcd b c)"
  using assms
proof -
  define q1 :: nat where "q1 = a div gcd a d" 
  define q2 :: nat where "q2 = c div gcd b c" 
  define q3 :: nat where "q3 = b div gcd b c"
  define q4 :: nat where "q4 = d div gcd a d"
  
  have "coprime q1 q2" "coprime q3 q4" 
    unfolding q1_def q2_def q3_def q4_def
  proof -
    have "coprime (a div gcd a d) c" 
      using ‹coprime a c› coprime_mult_left_iff[of "a div gcd a d" "gcd a d" c]
            dvd_mult_div_cancel[OF gcd_dvd1, of a b] by simp
    then show "coprime (a div gcd a d) (c div gcd b c)"
      using coprime_mult_right_iff[of "a div gcd a d" "gcd b c" "c div gcd b c"]
          dvd_div_mult_self[OF gcd_dvd2[of b c]] by auto    
    have "coprime (b div gcd b c) d" 
      using ‹coprime b d› coprime_mult_left_iff[of "b div gcd b c" "gcd b c" d]
            dvd_mult_div_cancel[OF gcd_dvd1, of a b] by simp
    then show "coprime (b div gcd b c) (d div gcd a d)"
      using coprime_mult_right_iff[of "b div gcd b c" "gcd a d" "d div gcd a d"]
          dvd_div_mult_self[OF gcd_dvd2[of b c]] by auto 
  qed
  moreover have "coprime q1 q4" "coprime q3 q2"
    unfolding q1_def q2_def q3_def q4_def 
    using assms div_gcd_coprime by blast+
  ultimately have 1: "coprime (q1*q3) (q2*q4)" 
    by simp 
  have "gcd (a*b) (c*d) = (gcd a d) * (gcd b c) * gcd (q1*q3) (q2*q4)"
    unfolding q1_def q2_def q3_def q4_def
    by (subst gcd_mult_distrib_nat[of "gcd a d * gcd b c"],
       simp add: field_simps,
       simp add: mult.left_commute semiring_normalization_rules(18))
  from this 1 show "gcd (a*b) (c*d) = (gcd a d) * (gcd b c)" by auto
qed  

lemma reindex_product_bij:
  fixes a b m k :: nat
  defines "S ≡ {(d1,d2). d1 dvd gcd a m ∧ d2 dvd gcd k b}"
  defines "T ≡ {d. d dvd (gcd a m) * (gcd k b)}"
  defines "f ≡ (λ(d1,d2). d1 * d2)"
  assumes "coprime a k"
  shows "bij_betw f S T"
  unfolding bij_betw_def
proof
  show inj: "inj_on f S"
    unfolding f_def 
  proof -
    {fix d1 d2 d1' d2'
    assume "(d1,d2) ∈ S" "(d1',d2') ∈ S"
    then have dvd: "d1 dvd gcd a m" "d2 dvd gcd k b" 
              "d1' dvd gcd a m" "d2' dvd gcd k b"
      unfolding S_def by simp+
    assume "f (d1,d2) = f (d1',d2')" 
    then have eq: "d1 * d2 = d1' * d2'" 
      unfolding f_def by simp
    from eq dvd have eq1: "d1 = d1'" 
      by (simp,meson assms coprime_crossproduct_nat coprime_divisors)
    from eq dvd have eq2: "d2 = d2'"
      using assms(4) eq1 by auto
    from eq1 eq2 have "d1 = d1' ∧ d2 = d2'" by simp} 
   then show "inj_on (λ(d1, d2). d1 * d2) S" 
    using S_def f_def by (intro inj_onI,blast)
  qed
  show surj: "f ` S = T" 
  proof -
    {fix d
      have "d dvd (gcd a m) * (gcd k b)
       ⟷ (∃d1 d2. d = d1*d2 ∧ d1 dvd gcd a m ∧ d2 dvd gcd k b)"
        using division_decomp mult_dvd_mono by blast}
      then show ?thesis
        unfolding f_def S_def T_def image_def
        by auto
  qed
qed

lemma p_div_set:
  shows "{p. p ∈prime_factors a ∧ ¬ p dvd N} = 
         ({p. p ∈prime_factors (a*N)} - {p. p ∈prime_factors N})"
  (is "?A = ?B") 
proof 
  show "?A ⊆ ?B" 
  proof (simp)
    { fix p 
      assume as: "p ∈# prime_factorization a" "¬ p dvd N"
      then have 1: "p ∈ prime_factors (a * N)" 
      proof -
        from in_prime_factors_iff[of p a] as
        have "a ≠ 0" "p dvd a" "prime p" by simp+
        have "N ≠ 0" using ‹¬ p dvd N› by blast
        have "a * N ≠ 0" using ‹a ≠ 0› ‹N ≠ 0› by auto
        have "p dvd a*N" using ‹p dvd a› by simp
        show ?thesis 
          using ‹a*N ≠ 0› ‹p dvd a*N› ‹prime p› in_prime_factors_iff by blast
      qed
      from as have 2: "p ∉ prime_factors N" by blast
      from 1 2 have "p ∈ prime_factors (a * N) - prime_factors N"
       by blast 
    }
    then show "{p. p ∈# prime_factorization a ∧ ¬ p dvd N}
               ⊆ prime_factors (a * N) - prime_factors N" by blast
  qed

  show "?B ⊆ ?A"
  proof (simp)
    { fix p 
      assume as: "p ∈ prime_factors (a * N) - prime_factors N"
      then have 1: "¬ p dvd N" 
      proof -
        from as have "p ∈ prime_factors (a * N)" "p ∉ prime_factors N"
          using DiffD1 DiffD2 by blast+
        then show ?thesis by (simp add: in_prime_factors_iff)        
      qed
      have 2: "p ∈# prime_factorization a" 
      proof -
        have "p dvd (a*N)" "prime p" "a*N ≠ 0" using in_prime_factors_iff as by blast+
        have "p dvd a" using ‹¬ p dvd N› prime_dvd_multD[OF ‹prime p› ‹p dvd (a*N)›] by blast
        have "a ≠ 0" using ‹a*N ≠ 0› by simp
        show ?thesis using in_prime_factors_iff ‹a ≠ 0› ‹p dvd a› ‹prime p› by blast
      qed
      from 1 2 have "p ∈ {p. p ∈# prime_factorization a ∧ ¬ p dvd N}" by blast
    }
    then show "prime_factors (a * N) - prime_factors N
               ⊆ {p. p ∈# prime_factorization a ∧ ¬ p dvd N}" by blast
  qed
qed

lemma coprime_iff_prime_factors_disjoint:
  fixes x y :: "'a :: factorial_semiring"
  assumes "x ≠ 0" "y ≠ 0"
  shows "coprime x y ⟷ prime_factors x ∩ prime_factors y = {}"
proof 
  assume "coprime x y"
  have False if "p ∈ prime_factors x" "p ∈ prime_factors y" for p
  proof -
    from that assms have "p dvd x" "p dvd y"
      by (auto simp: prime_factors_dvd)
    with ‹coprime x y› have "p dvd 1"
      using coprime_common_divisor by auto
    with that assms show False by (auto simp: prime_factors_dvd)
  qed
  thus "prime_factors x ∩ prime_factors y = {}" by auto
next
  assume disjoint: "prime_factors x ∩ prime_factors y = {}"
  show "coprime x y"
  proof (rule coprimeI)
    fix d assume d: "d dvd x" "d dvd y"
    show "is_unit d"
    proof (rule ccontr)
      assume "¬is_unit d"
      moreover from this and d assms have "d ≠ 0" by auto
      ultimately obtain p where p: "prime p" "p dvd d"
        using prime_divisor_exists by auto
      with d and assms have "p ∈ prime_factors x ∩ prime_factors y"
        by (auto simp: prime_factors_dvd)
      with disjoint show False by auto
    qed
  qed
qed

lemma coprime_cong_prime_factors:
  fixes x y :: "'a :: factorial_semiring_gcd"
  assumes "x ≠ 0" "y ≠ 0" "x' ≠ 0" "y' ≠ 0"
  assumes "prime_factors x = prime_factors x'"
  assumes "prime_factors y = prime_factors y'"
  shows   "coprime x y ⟷ coprime x' y'"
  using assms by (simp add: coprime_iff_prime_factors_disjoint)

lemma moebius_prod_not_coprime:
  assumes "¬ coprime N d"
  shows "moebius_mu (N*d) = 0"
proof -
  from assms obtain l where l_form: "l dvd N ∧ l dvd d ∧ ¬ is_unit l" 
    unfolding coprime_def by blast
  then have "l * l dvd N * d" using mult_dvd_mono by auto
  then have "l⇧2 dvd N*d" by (subst power2_eq_square,blast) 
  then have "¬ squarefree (N*d)" 
    unfolding squarefree_def coprime_def using l_form by blast
  then show "moebius_mu (N*d) = 0" 
    using moebius_mu_def by auto 
qed

text‹Theorem 2.18›

(* TODO Place in corresponding theory *)
lemma sum_divisors_moebius_mu_times_multiplicative:
  fixes f :: "nat ⇒ 'a :: {comm_ring_1}"
  assumes "multiplicative_function f" and "n > 0"
  shows   "(∑d | d dvd n. moebius_mu d * f d) = (∏p∈prime_factors n. 1 - f p)"
proof -
  define g where "g = (λn. ∑d | d dvd n. moebius_mu d * f d)"
  define g' where "g' = dirichlet_prod (λn. moebius_mu n * f n) (λn. if n = 0 then 0 else 1)"
  interpret f: multiplicative_function f by fact
  have "multiplicative_function (λn. if n = 0 then 0 else 1 :: 'a)"
    by standard auto
  interpret multiplicative_function g' unfolding g'_def
    by (intro multiplicative_dirichlet_prod multiplicative_function_mult
              moebius_mu.multiplicative_function_axioms assms) fact+

  have g'_primepow: "g' (p ^ k) = 1 - f p" if "prime p" "k > 0" for p k
  proof -
    have "g' (p ^ k) = (∑i≤k. moebius_mu (p ^ i) * f (p ^ i))"
      using that by (simp add: g'_def dirichlet_prod_prime_power)
    also have "… = (∑i∈{0, 1}. moebius_mu (p ^ i) * f (p ^ i))"
      using that by (intro sum.mono_neutral_right) (auto simp: moebius_mu_power')
    also have "… = 1 - f p"
      using that by (simp add: moebius_mu.prime)
    finally show ?thesis .
  qed

  have "g' n = g n"
    by (simp add: g_def g'_def dirichlet_prod_def)
  also from assms have "g' n = (∏p∈prime_factors n. g' (p ^ multiplicity p n))"
   by (intro prod_prime_factors) auto
  also have "… = (∏p∈prime_factors n. 1 - f p)"
    by (intro prod.cong) (auto simp: g'_primepow prime_factors_multiplicity)
  finally show ?thesis by (simp add: g_def)
qed

lemma multiplicative_ind_coprime [intro]: "multiplicative_function (ind (coprime N))"
  by (intro multiplicative_function_ind) auto

lemma sum_divisors_moebius_mu_times_multiplicative_revisited:
  fixes f :: "nat ⇒ 'a :: {comm_ring_1}"
  assumes "multiplicative_function f" "n > 0" "N > 0"
  shows   "(∑d | d dvd n ∧ coprime N d. moebius_mu d * f d) =
           (∏p∈{p. p ∈ prime_factors n ∧ ¬ (p dvd N)}. 1 - f p)"
proof -
  have "(∑d | d dvd n ∧ coprime N d. moebius_mu d * f d) =
          (∑d | d dvd n. moebius_mu d * (ind (coprime N) d * f d))"
    using assms by (intro sum.mono_neutral_cong_left) (auto simp: ind_def)
  also have "… = (∏p∈prime_factors n. 1 - ind (coprime N) p * f p)"
    using assms by (intro sum_divisors_moebius_mu_times_multiplicative)
                   (auto intro: multiplicative_function_mult)
  also from assms have "… = (∏p | p ∈ prime_factors n ∧ ¬(p dvd N). 1 - f p)"
    by (intro prod.mono_neutral_cong_right)
       (auto simp: ind_def prime_factors_dvd coprime_commute dest: prime_imp_coprime)
  finally show ?thesis .
qed

subsection ‹Neutral element of the Dirichlet product›

definition "dirichlet_prod_neutral n = (if n = 1 then 1 else 0)" for n :: nat

lemma dirichlet_prod_neutral_intro: 
  fixes S :: "nat ⇒ complex" and f :: "nat ⇒ nat ⇒ complex"
  defines "S ≡ (λ(n::nat). (∑k | k ∈ {1..n} ∧ coprime k n. (f k n)))" 
  shows "S(n) = (∑k ∈ {1..n}. f k n * dirichlet_prod_neutral (gcd k n))"
proof -
  let ?g = "λk. (f k n)* (dirichlet_prod_neutral (gcd k n))"
  have zeros: "∀k ∈ {1..n} - {k. k ∈ {1..n} ∧ coprime k n}. ?g k = 0"
  proof
    fix k 
    assume "k ∈ {1..n} - {k ∈ {1..n}. coprime k n}"
    then show "(f k n) * dirichlet_prod_neutral (gcd k n) = 0" 
      by (simp add: dirichlet_prod_neutral_def[of "gcd k n"] split: if_splits,presburger) 
  qed
  
  have "S n = (∑k | k ∈ {1..n} ∧ coprime k n. (f k n))"
    by (simp add: S_def)
  also have "… = sum ?g {k. k ∈ {1..n} ∧ coprime k n}"
    by (simp add: dirichlet_prod_neutral_def split: if_splits)    
  also have "… = sum ?g {1..n}"
    by (intro sum.mono_neutral_left, auto simp add: zeros)
  finally show ?thesis by blast 
qed

lemma dirichlet_prod_neutral_right_neutral:
 "dirichlet_prod f dirichlet_prod_neutral n = f n " if "n > 0" for f :: "nat ⇒ complex" and n 
proof -
  {fix d :: nat
    assume "d dvd n"
    then have eq: "n = d ⟷ n div d = 1"
      using div_self that dvd_mult_div_cancel by force
    have "f(d)*dirichlet_prod_neutral(n div d) = (if n = d then f(d) else 0)" 
      by (simp add: dirichlet_prod_neutral_def eq)}
  note summand = this

  have "dirichlet_prod f dirichlet_prod_neutral n = 
          (∑d | d dvd n. f(d)*dirichlet_prod_neutral(n div d))"
    unfolding dirichlet_prod_def by blast
  also have "… = (∑d | d dvd n. (if n = d then f(d) else 0))"
    using summand by simp
  also have "… = (∑d | d = n. (if n = d then f(d) else 0))"
    using that by (intro sum.mono_neutral_right, auto)
  also have "… = f(n)"  by simp
  finally show ?thesis by simp 
qed

lemma dirichlet_prod_neutral_left_neutral:
 "dirichlet_prod dirichlet_prod_neutral f n = f n " 
 if "n > 0" for f :: "nat ⇒ complex" and n
  using dirichlet_prod_neutral_right_neutral[OF that, of f] 
        dirichlet_prod_commutes[of f dirichlet_prod_neutral]
  by argo

corollary I_right_neutral_0:
  fixes f :: "nat ⇒ complex" 
  assumes "f 0 = 0"
  shows "dirichlet_prod f dirichlet_prod_neutral n = f n"
  using assms dirichlet_prod_neutral_right_neutral by (cases n, simp, blast)

subsection ‹Multiplicative functions›

lemma mult_id: "multiplicative_function id" 
  by (simp add: multiplicative_function_def)

lemma mult_moebius: "multiplicative_function moebius_mu"
  using Moebius_Mu.moebius_mu.multiplicative_function_axioms
  by simp

lemma mult_of_nat: "multiplicative_function of_nat" 
  using multiplicative_function_def of_nat_0 of_nat_1 of_nat_mult by blast

lemma mult_of_nat_c: "completely_multiplicative_function of_nat" 
  by (simp add: completely_multiplicative_function_def)

lemma completely_multiplicative_nonzero:
  fixes f :: "nat ⇒ complex"
  assumes "completely_multiplicative_function f" 
          "d ≠ 0"
          "⋀p. prime p ⟹ f(p) ≠ 0" 
  shows "f(d) ≠ 0"
  using assms(2)
proof (induction d rule: nat_less_induct)
  case (1 n)
  then show ?case 
  proof (cases "n = 1")
    case True
    then show ?thesis 
      using assms(1)
      unfolding completely_multiplicative_function_def by simp
  next
    case False
    then obtain p where 2:"prime p ∧ p dvd n" 
      using prime_factor_nat by blast
    then obtain a where 3: "n = p * a"  "a ≠ 0" 
      using 1 by auto
    then have 4: "f(a) ≠ 0" using 1 
      using 2 prime_nat_iff by fastforce
    have 5: "f(p) ≠ 0" using assms(3) 2 by simp
    from 3 4 5 show ?thesis
      by (simp add: assms(1) completely_multiplicative_function.mult)
  qed
qed

lemma multipl_div: 
  fixes m k d1 d2 :: nat and f :: "nat ⇒ complex"
  assumes "multiplicative_function f" "d1 dvd m" "d2 dvd k" "coprime m k"
  shows "f ((m*k) div (d1*d2)) = f(m div d1) * f(k div d2)"
  using assms
  unfolding multiplicative_function_def
  using assms(1) multiplicative_function.mult_coprime by fastforce

lemma multipl_div_mono: 
  fixes m k d :: nat and f :: "nat ⇒ complex"
  assumes "completely_multiplicative_function f" 
          "d dvd k" "d > 0" 
          "⋀p. prime p ⟹ f(p) ≠ 0" 
  shows "f (k div d) = f(k) div f(d)"
proof -
  have "d ≠ 0" using assms(2,3) by auto
  then have nz: "f(d) ≠ 0" using assms(1,4) completely_multiplicative_nonzero by simp

  from assms(2,3) obtain a where div: "k = a * d " by fastforce
  have "f (k div d) = f ((a*d) div d)" using div by simp
  also have "… = f(a)" using assms(3) div by simp
  also have "… = f(a)*f(d) div f(d)" using nz by auto
  also have "… = f(a*d) div f(d)" 
    by (simp add: div assms(1) completely_multiplicative_function.mult)
  also have "… = f (k) div f(d)" using div by simp
  finally show ?thesis by simp
qed

lemma comp_to_mult: "completely_multiplicative_function f ⟹
       multiplicative_function f"
  unfolding completely_multiplicative_function_def
            multiplicative_function_def by auto

end

(*
  File:     Periodic_Arithmetic.thy
  Authors:  Rodrigo Raya, EPFL; Manuel Eberl, TUM

  Periodic arithmetic functions
*)
section ‹Periodic arithmetic functions›
theory Periodic_Arithmetic
imports
  Complex_Main
  "HOL-Number_Theory.Cong"
begin

definition 
  "periodic_arithmetic f k = (∀n. f (n+k) = f n)" 
  for n :: int and k :: nat and f :: "nat ⇒ complex"

lemma const_periodic_arithmetic: "periodic_arithmetic (λx. y) k"
  unfolding periodic_arithmetic_def by blast

lemma add_periodic_arithmetic:
  fixes f g :: "nat ⇒ complex"
  assumes "periodic_arithmetic f k"
  assumes "periodic_arithmetic g k"
  shows "periodic_arithmetic (λn. f n + g n) k"
  using assms unfolding periodic_arithmetic_def by simp

lemma mult_periodic_arithmetic:
  fixes f g :: "nat ⇒ complex"
  assumes "periodic_arithmetic f k"
  assumes "periodic_arithmetic g k"
  shows "periodic_arithmetic (λn. f n * g n) k"
  using assms unfolding periodic_arithmetic_def  by simp

lemma scalar_mult_periodic_arithmetic:
  fixes f :: "nat ⇒ complex" and a :: complex
  assumes "periodic_arithmetic f k"
  shows "periodic_arithmetic (λn. a * f n) k"
  using mult_periodic_arithmetic[OF const_periodic_arithmetic[of a k] assms(1)] by simp

lemma fin_sum_periodic_arithmetic_set:
  fixes f g :: "nat ⇒ complex" 
  assumes "∀i∈A. periodic_arithmetic (h i) k"
  shows "periodic_arithmetic (λn. ∑i ∈ A. h i n) k"
  using assms by (simp add: periodic_arithmetic_def)

lemma mult_period:
  assumes "periodic_arithmetic g k"
  shows "periodic_arithmetic g (k*q)"
  using assms
proof (induction q)
  case 0 then show ?case unfolding periodic_arithmetic_def by simp
next
  case (Suc m)
  then show ?case 
    unfolding periodic_arithmetic_def 
  proof -
   { fix n 
     have "g (n + k * Suc m) = g (n + k +  k * m)"
       by (simp add: algebra_simps)
     also have "… = g(n)" 
       using Suc.IH[OF Suc.prems] assms
       unfolding periodic_arithmetic_def by simp
     finally have "g (n + k * Suc m) = g(n)" by blast
   }
    then show "∀n. g (n + k * Suc m) = g n" by auto
  qed   
qed

lemma unique_periodic_arithmetic_extension:
  assumes "k > 0"
  assumes "∀j<k. g j = h j"
  assumes "periodic_arithmetic g k" and "periodic_arithmetic h k"
  shows "g i = h i"
proof (cases "i < k")
  case True then show ?thesis using assms by simp
next
  case False then show ?thesis 
  proof -
    have "k * (i div k) + (i mod k) = i ∧ (i mod k) < k" 
      by (simp add: assms(1) algebra_simps)
    then obtain q r where euclid_div: "k*q + r = i ∧ r < k"
      using mult.commute by blast
    from assms(3) assms(4) 
    have  "periodic_arithmetic g (k*q)" "periodic_arithmetic h (k*q)" 
      using mult_period by simp+
    have "g(k*q+r) = g(r)" 
      using ‹periodic_arithmetic g (k*q)› unfolding periodic_arithmetic_def 
      using add.commute[of "k*q" r] by presburger
    also have "… = h(r)" 
      using euclid_div assms(2) by simp
    also have "… = h(k*q+r)"
      using ‹periodic_arithmetic h (k*q)› add.commute[of "k*q" r]
      unfolding periodic_arithmetic_def by presburger
    also have "… = h(i)" using euclid_div by simp
    finally show "g(i) = h(i)" using euclid_div by simp
  qed
qed
                  
lemma periodic_arithmetic_sum_periodic_arithmetic:
  assumes "periodic_arithmetic f k"
  shows "(∑l ∈ {m..n}. f l) = (∑l ∈ {m+k..n+k}. f l)"
  using periodic_arithmetic_def assms 
  by (intro sum.reindex_bij_witness
         [of "{m..n}" "λl. l-k" "λl. l+k" "{m+k..n+k}" f f])
      auto

lemma mod_periodic_arithmetic:
  fixes n m :: nat
  assumes "periodic_arithmetic f k"
  assumes "n mod k = m mod k"
  shows "f n = f m"
proof -
  obtain q where 1: "n = q*k+(n mod k)"   
     using div_mult_mod_eq[of n k,symmetric] by blast 
  obtain q' where 2: "m = q'*k+(m mod k)"
     using div_mult_mod_eq[of m k,symmetric] by blast
  from 1 have "f n = f (q*k+(n mod k))" by auto
  also have "… = f (n mod k)"
    using mult_period[of f k q] assms(1) periodic_arithmetic_def[of f "k*q"]
    by (simp add: algebra_simps,subst add.commute,blast)
  also have "… = f (m mod k)" using assms(2) by auto
  also have "… = f (q'*k+(m mod k))"
    using mult_period[of f k q'] assms(1) periodic_arithmetic_def[of f "k*q'"]
    by (simp add: algebra_simps,subst add.commute,presburger)
  also have "… = f m" using 2 by auto
  finally show "f n = f m" by simp
qed

lemma cong_periodic_arithmetic:
  assumes "periodic_arithmetic f k" "[a = b] (mod k)"
  shows   "f a = f b"
  using assms mod_periodic_arithmetic[of f k a b] by (auto simp: cong_def)

lemma cong_nat_imp_eq:
  fixes m :: nat
  assumes "m > 0" "x ∈ {a..<a+m}" "y ∈ {a..<a+m}" "[x = y] (mod m)"
  shows   "x = y"
  using assms
proof (induction x y rule: linorder_wlog)
  case (le x y)
  have "[y - x = 0] (mod m)"
    using cong_diff_iff_cong_0_nat cong_sym le by blast
  thus "x = y"
    using le by (auto simp: cong_def)
qed (auto simp: cong_sym)

lemma inj_on_mod_nat:
  fixes m :: nat
  assumes "m > 0"
  shows   "inj_on (λx. x mod m) {a..<a+m}"
proof
  fix x y assume xy: "x ∈ {a..<a+m}" "y ∈ {a..<a+m}" and eq: "x mod m = y mod m"
  from ‹m > 0› and xy show "x = y"
    by (rule cong_nat_imp_eq) (use eq in ‹simp_all add: cong_def›)
qed

lemma bij_betw_mod_nat_atLeastLessThan:
  fixes k d :: nat
  assumes "k > 0"
  defines "g ≡ (λi. nat ((int i - int d) mod int k) + d)"
  shows   "bij_betw (λi. i mod k) {d..<d+k} {..<k}"
  unfolding bij_betw_def
proof
  show inj: "inj_on (λi. i mod k) {d..<d + k}"
    by (rule inj_on_mod_nat) fact+
  have "(λi. i mod k) ` {d..<d + k} ⊆ {..<k}"
    by auto
  moreover have "card ((λi. i mod k) ` {d..<d + k}) = card {..<k}"
    using inj by (subst card_image) auto
  ultimately show "(λi. i mod k) ` {d..<d + k} = {..<k}"
    by (intro card_subset_eq) auto
qed

lemma periodic_arithmetic_sum_periodic_arithmetic_shift:
  fixes k d :: nat
  assumes "periodic_arithmetic f k" "k > 0" "d > 0"
  shows "(∑l ∈ {0..k-1}. f l) = (∑l ∈ {d..d+k-1}. f l)"
proof -
  have "(∑l ∈ {0..k-1}. f l) = (∑l ∈ {0..<k}. f l)"
    using assms(2) by (intro sum.cong) auto
  also have "… = (∑l ∈ {d..<d+k}. f (l mod k))"
    using assms(2) 
    by (simp add: sum.reindex_bij_betw[OF bij_betw_mod_nat_atLeastLessThan[of k d]] 
                  lessThan_atLeast0)
  also have "… = (∑l ∈ {d..<d+k}. f l)"
    using mod_periodic_arithmetic[of f k] assms(1) sum.cong
    by (meson mod_mod_trivial)
  also have "… = (∑l ∈ {d..d+k-1}. f l)"
    using assms(2,3) by (intro sum.cong) auto
  finally show ?thesis by auto
qed

lemma self_bij_0_k:
  fixes a k :: nat
  assumes "coprime a k" "[a*i = 1] (mod k)" "k > 0" 
  shows "bij_betw (λr. r*a mod k) {0..k-1} {0..k-1}"
  unfolding bij_betw_def
proof
  show "inj_on (λr. r*a mod k) {0..k-1}"
  proof -
    {fix r1 r2
    assume in_k: "r1 ∈ {0..k-1}" "r2 ∈ {0..k-1}"
    assume as: "[r1*a = r2*a] (mod k)"
    then have "[r1*a*i = r2*a*i] (mod k)" 
      using cong_scalar_right by blast
    then have "[r1 = r2] (mod k)" 
      using cong_mult_rcancel_nat as assms(1) by simp
    then have "r1 = r2" using in_k
      using assms(3) cong_less_modulus_unique_nat by auto}
    note eq = this
    show ?thesis unfolding inj_on_def 
      by (safe, simp add: eq cong_def)
  qed
  define f where "f = (λr. r * a mod k)"
  show "f ` {0..k - 1} = {0..k - 1} "
    unfolding image_def
  proof (standard)
    show "{y. ∃x∈{0..k - 1}. y = f x} ⊆ {0..k - 1}" 
    proof -
      {fix y
      assume "y ∈ {y. ∃x∈{0..k - 1}. y = f x}" 
      then obtain x where "y = f x" by blast
      then have "y ∈ {0..k-1}"
        unfolding f_def
        using Suc_pred assms(3) lessThan_Suc_atMost by fastforce}
      then show ?thesis by blast
    qed
    show "{0..k - 1} ⊆ {y. ∃x∈{0..k - 1}. y = f x}"
    proof -
      { fix x 
        assume ass: "x ∈ {0..k-1}"
        then have "x * i mod k ∈ {0..k-1}"
        proof -
          have "x * i mod k ∈ {0..<k}" by (simp add: assms(3))
          have "{0..<k} = {0..k-1}" using Suc_diff_1 assms(3) by auto
          show ?thesis using ‹x * i mod k ∈ {0..<k}› ‹{0..<k} = {0..k-1}› by blast
        qed          
        then have "f (x * i mod k) = x"
        proof -
          have "f (x * i mod k) = (x * i mod k) * a mod k"
            unfolding f_def by blast
          also have "… = (x*i*a) mod k" 
            by (simp add: mod_mult_left_eq) 
          also have "… = (x*1) mod k" 
            using assms(2) 
            unfolding cong_def 
            by (subst mult.assoc, subst (2) mult.commute,
               subst mod_mult_right_eq[symmetric],simp) 
          also have "… = x" using ass assms(3) by auto
          finally show ?thesis .
        qed
        then have "x ∈ {y. ∃x∈{0..k - 1}. y = f x}" 
          using ‹x * i mod k ∈ {0..k-1}› by force
      }
      then show ?thesis by blast 
    qed
  qed
qed

lemma periodic_arithmetic_homothecy:
  assumes "periodic_arithmetic f k"
  shows   "periodic_arithmetic (λl. f (l*a)) k"
  unfolding periodic_arithmetic_def
proof 
  fix n
  have "f ((n + k) * a) = f(n*a+k*a)" by (simp add: algebra_simps)
  also have "… = f(n*a)" 
    using mult_period[OF assms] unfolding periodic_arithmetic_def by simp
  finally show "f ((n + k) * a) = f (n * a)" by simp
qed

theorem periodic_arithmetic_remove_homothecy:
  assumes "coprime a k" "periodic_arithmetic f k" "k > 0" 
  shows "(∑l=1..k. f l) = (∑l=1..k. f (l*a))" 
proof -
  obtain i where inv: "[a*i = 1] (mod k)"
    using assms(1) coprime_iff_invertible_nat[of a k] by auto
  from this self_bij_0_k assms
  have bij: "bij_betw (λr. r * a mod k) {0..k - 1} {0..k - 1}" by blast
  
  have "(∑l = 1..k. f(l)) = (∑l = 0..k-1. f(l))"
    using periodic_arithmetic_sum_periodic_arithmetic_shift[of f k 1] assms by simp
  also have "… = (∑l = 0..k-1. f(l*a mod k))"
    using sum.reindex_bij_betw[OF bij,symmetric] by blast
  also have "… = (∑l = 0..k-1. f(l*a))"
    by (intro sum.cong refl) (use mod_periodic_arithmetic[OF assms(2)] mod_mod_trivial in blast)
  also have "… = (∑l = 1..k. f(l*a))"
    using periodic_arithmetic_sum_periodic_arithmetic_shift[of "(λl. f(l*a))" k 1]
          periodic_arithmetic_homothecy[OF assms(2)] assms(3) by fastforce  
  finally show ?thesis by blast     
qed

end

(*
  File:     Complex_Roots_Of_Unity.thy
  Authors:  Rodrigo Raya, EPFL; Manuel Eberl, TUM

  Complex roots of unity (exp(2iπ/n)) and sums over them.
*)
theory Complex_Roots_Of_Unity
imports
  "HOL-Analysis.Analysis"
  Periodic_Arithmetic
begin

section ‹Complex roots of unity›

definition 
  "unity_root k n = cis (2 * pi * of_int n / of_nat k)"

lemma 
  unity_root_k_0 [simp]: "unity_root k 0 = 1" and
  unity_root_0_n [simp]: "unity_root 0 n = 1"
  unfolding unity_root_def by simp+

lemma unity_root_conv_exp: 
  "unity_root k n = exp (of_real (2*pi*n/k) * 𝗂)"
  unfolding unity_root_def 
  by (subst cis_conv_exp,subst mult.commute,blast)

lemma unity_root_mod: 
  "unity_root k (n mod int k) = unity_root k n"
proof (cases "k = 0")
  case True then show ?thesis by simp
next
  case False
  obtain q :: int where q_def: "n = q*k + (n mod k)" 
    using div_mult_mod_eq[symmetric] by blast
  have "n / k = q + (n mod k) / k"
  proof (auto simp add: divide_simps False)
    have "real_of_int n = real_of_int (q*k + (n mod k))"
      using q_def by simp
    also have "… = real_of_int q * real k + real_of_int (n mod k)"
      using of_int_add of_int_mult by simp
    finally show "real_of_int n = real_of_int q * real k + real_of_int (n mod k)" 
      by blast
  qed
  then have "(2*pi*n/k) = 2*pi*q + (2*pi*(n mod k)/k)"
    using False by (auto simp add: field_simps)
  then have "(2*pi*n/k)*𝗂 = 2*pi*q*𝗂 + (2*pi*(n mod k)/k)*𝗂" (is "?l = ?r1 + ?r2")
    by (auto simp add: algebra_simps)
  then have "exp ?l = exp ?r2"
    using exp_plus_2pin by (simp add: exp_add mult.commute)
  then show ?thesis 
    using unity_root_def unity_root_conv_exp by simp
qed

lemma unity_root_cong:
  assumes "[m = n] (mod int k)"
  shows   "unity_root k m = unity_root k n"
proof -
  from assms have "m mod int k = n mod int k"
    by (auto simp: cong_def)
  hence "unity_root k (m mod int k) = unity_root k (n mod int k)"
    by simp
  thus ?thesis by (simp add: unity_root_mod)
qed

lemma unity_root_mod_nat: 
  "unity_root k (nat (n mod int k)) = unity_root k n"
proof (cases k)
  case (Suc l)
  then have "n mod int k ≥ 0" by auto
  show ?thesis 
    unfolding int_nat_eq 
    by (simp add: ‹n mod int k ≥ 0› unity_root_mod)
qed auto

lemma unity_root_eqD:
 assumes gr: "k > 0"
 assumes eq: "unity_root k i = unity_root k j"
 shows "i mod k = j mod k"
proof - 
  let ?arg1 = "(2*pi*i/k)* 𝗂"
  let ?arg2 = "(2*pi*j/k)* 𝗂"
  from eq unity_root_conv_exp have "exp ?arg1 = exp ?arg2" by simp
  from this exp_eq 
  obtain n :: int where "?arg1 = ?arg2 +(2*n*pi)*𝗂" by blast
  then have e1: "?arg1 - ?arg2 = 2*n*pi*𝗂" by simp
  have e2: "?arg1 - ?arg2 = 2*(i-j)*(1/k)*pi*𝗂"
    by (auto simp add: algebra_simps)
  from e1 e2 have "2*n*pi*𝗂 = 2*(i-j)*(1/k)*pi*𝗂" by simp
  then have "2*n*k*pi*𝗂 = 2*(i-j)*pi*𝗂"
    by (simp add: divide_simps ‹k > 0›)(simp add: field_simps)
  then have "2*n*k = 2*(i-j)"
    by (meson complex_i_not_zero mult_cancel_right of_int_eq_iff of_real_eq_iff pi_neq_zero)
  then have "n*k = i-j" by auto
  then show ?thesis by Groebner_Basis.algebra
qed

lemma unity_root_eq_1_iff:
  fixes k n :: nat
  assumes "k > 0" 
  shows "unity_root k n = 1 ⟷ k dvd n"
proof -
  have "unity_root k n = exp ((2*pi*n/k) * 𝗂)"
    by (simp add: unity_root_conv_exp)
  also have "exp ((2*pi*n/k)* 𝗂) = 1 ⟷ k dvd n"
    using complex_root_unity_eq_1[of k n] assms
    by (auto simp add: algebra_simps)
  finally show ?thesis by simp
qed

lemma unity_root_pow: "unity_root k n ^ m = unity_root k (n * m)"
  using unity_root_def
  by (simp add: Complex.DeMoivre mult.commute algebra_split_simps(6))

lemma unity_root_add: "unity_root k (m + n) = unity_root k m * unity_root k n"
  by (simp add: unity_root_conv_exp add_divide_distrib algebra_simps exp_add)

lemma unity_root_uminus: "unity_root k (-m) = cnj (unity_root k m)"
  unfolding unity_root_conv_exp exp_cnj by simp

lemma inverse_unity_root: "inverse (unity_root k m) = cnj (unity_root k m)" 
  unfolding unity_root_conv_exp exp_cnj by (simp add: field_simps exp_minus)

lemma unity_root_diff: "unity_root k (m - n) = unity_root k m * cnj (unity_root k n)"
  using unity_root_add[of k m "-n"] by (simp add: unity_root_uminus)

lemma unity_root_eq_1_iff_int:
  fixes k :: nat and n :: int
  assumes "k > 0" 
  shows "unity_root k n = 1 ⟷ k dvd n"
proof (cases "n ≥ 0")
  case True
  obtain n' where "n = int n'" 
    using zero_le_imp_eq_int[OF True] by blast
  then show ?thesis 
    using unity_root_eq_1_iff[OF ‹k > 0›, of n'] of_nat_dvd_iff by blast
next
  case False
  then have "-n ≥ 0" by auto
  have "unity_root k n = inverse (unity_root k (-n))"
    unfolding inverse_unity_root by (simp add: unity_root_uminus)
  then have "(unity_root k n = 1) = (unity_root k (-n) = 1)" 
    by simp
  also have "(unity_root k (-n) = 1) = (k dvd (-n))"
    using unity_root_eq_1_iff[of k "nat (-n)",OF ‹k > 0›] False 
          int_dvd_int_iff[of k "nat (-n)"] nat_0_le[OF ‹-n ≥ 0›] by auto
  finally show ?thesis by simp
qed

lemma unity_root_eq_1 [simp]: "int k dvd n ⟹ unity_root k n = 1"
  by (cases "k = 0") (auto simp: unity_root_eq_1_iff_int)

lemma unity_periodic_arithmetic:
  "periodic_arithmetic (unity_root k) k" 
  unfolding periodic_arithmetic_def 
proof
  fix n
  have "unity_root k (n + k) = unity_root k ((n+k) mod k)"
    using unity_root_mod[of k] zmod_int by presburger
  also have "unity_root k ((n+k) mod k) = unity_root k n" 
    using unity_root_mod zmod_int by auto
  finally show "unity_root k (n + k) = unity_root k n" by simp
qed

lemma unity_periodic_arithmetic_mult:
  "periodic_arithmetic (λn. unity_root k (m * int n)) k"
  unfolding periodic_arithmetic_def
proof 
  fix n
  have "unity_root k (m * int (n + k)) = 
        unity_root k (m*n + m*k)"
    by (simp add: algebra_simps)
  also have "… = unity_root k (m*n)"
    using unity_root_mod[of k "m * int n"] unity_root_mod[of k "m * int n + m * int k"] 
          mod_mult_self3 by presburger
  finally show "unity_root k (m * int (n + k)) =
             unity_root k (m * int n)" by simp
qed

lemma unity_root_periodic_arithmetic_mult_minus:
  shows "periodic_arithmetic (λi. unity_root k (-int i*int m)) k" 
  unfolding periodic_arithmetic_def
proof 
  fix n 
  have "unity_root k (-(n + k) * m) = cnj (unity_root k (n*m+k*m))" 
    by (simp add: ring_distribs unity_root_diff unity_root_add unity_root_uminus)
  also have "… = cnj (unity_root k (n*m))"
    using mult_period[of "unity_root k" k m] unity_periodic_arithmetic[of k]
    unfolding periodic_arithmetic_def by presburger
  also have "… = unity_root k (-n*m)"
    by (simp add: unity_root_uminus)
  finally show "unity_root k (-(n + k) * m) = unity_root k (-n*m)"
    by simp
qed

lemma unity_div:
 fixes a :: int and d :: nat
 assumes "d dvd k"
 shows "unity_root k (a*d) = unity_root (k div d) a" 
proof -
  have 1: "(2*pi*(a*d)/k) = (2*pi*a)/(k div d)"
    using Suc_pred assms by (simp add: divide_simps, fastforce)   
  have "unity_root k (a*d) = exp ((2*pi*(a*d)/k)* 𝗂)"
    using unity_root_conv_exp by simp
  also have "… = exp (((2*pi*a)/(k div d))* 𝗂)"
    using 1 by simp
  also have "… = unity_root (k div d) a"
    using unity_root_conv_exp by simp
  finally show ?thesis by simp
qed

lemma unity_div_num:
  assumes "k > 0" "d > 0" "d dvd k"
  shows "unity_root k (x * (k div d)) = unity_root d x"
  using assms dvd_div_mult_self unity_div by auto


section ‹Geometric sums of roots of unity›

text‹
  Apostol calls these `geometric sums', which is a bit too generic. We therefore decided
  to refer to them as `sums of roots of unity'.
›
definition "unity_root_sum k n = (∑m<k. unity_root k (n * of_nat m))"

lemma unity_root_sum_0_left [simp]: "unity_root_sum 0 n = 0" and
      unity_root_sum_0_right [simp]: "k > 0 ⟹ unity_root_sum k 0 = k" 
  unfolding unity_root_sum_def by simp_all

text ‹Theorem 8.1›
theorem unity_root_sum:
  fixes k :: nat and n :: int
  assumes gr: "k ≥ 1"
  shows "k dvd n ⟹ unity_root_sum k n = k"
    and "¬k dvd n ⟹ unity_root_sum k n = 0"
proof -
  assume dvd: "k dvd n"
  let ?x = "unity_root k n"
  have unit: "?x = 1" using dvd gr unity_root_eq_1_iff_int by auto
  have exp: "?x^m = unity_root k (n*m)" for m using unity_root_pow by simp
  have "unity_root_sum k n = (∑m<k. unity_root k (n*m))" 
    using unity_root_sum_def by simp 
  also have "… = (∑m<k. ?x^m)" using exp by auto
  also have "… = (∑m<k. 1)" using unit by simp
  also have "… = k" using gr by (induction k, auto)
  finally show "unity_root_sum k n = k" by simp
next
  assume dvd: "¬k dvd n"
  let ?x = "unity_root k n"
  have "?x ≠ 1" using dvd gr unity_root_eq_1_iff_int by auto
  have "(?x^k - 1)/(?x - 1) = (∑m<k. ?x^m)"
    using geometric_sum[of ?x k, OF ‹?x ≠ 1›] by auto
  then have sum: "unity_root_sum k n = (?x^k - 1)/(?x - 1)"
    using unity_root_sum_def unity_root_pow by simp
  have "?x^k = 1" 
    using gr unity_root_eq_1_iff_int unity_root_pow by simp
  then show "unity_root_sum k n = 0" using sum by auto
qed

corollary unity_root_sum_periodic_arithmetic: 
 "periodic_arithmetic (unity_root_sum k) k"
  unfolding periodic_arithmetic_def
proof 
  fix n 
  show "unity_root_sum k (n + k) = unity_root_sum k n"
    by (cases "k = 0"; cases "k dvd n") (auto simp add: unity_root_sum)
qed

lemma unity_root_sum_nonzero_iff:
  fixes r :: int
  assumes "k ≥ 1" and "r ∈ {-k<..<k}"
  shows "unity_root_sum k r ≠ 0 ⟷ r = 0"
proof
  assume "unity_root_sum k r ≠ 0"
  then have "k dvd r" using unity_root_sum assms by blast
  then show "r = 0" using assms(2) 
    using dvd_imp_le_int by force
next
  assume "r = 0"
  then have "k dvd r" by auto
  then have "unity_root_sum k r = k" 
    using assms(1) unity_root_sum by blast
  then show "unity_root_sum k r ≠ 0" using assms(1) by simp
qed

end