Vanishing Sum of Rotated Powers

 Complex.exp (2 * π * Complex.I / N)
    ∑ k ∈ Finset.range N, (z * ω ^ k) ^ j = 0 := by
  intro ω
  have hN' : N ≠ 0 := Nat.pos_iff_ne_zero.mp hN
  have hω : IsPrimitiveRoot ω N := Complex.isPrimitiveRoot_exp N hN'
  simp_rw [mul_pow]
  by_cases hz : z = 0
  · have hj_ne : j ≠ 0 := Nat.pos_iff_ne_zero.mp hj
    simp [hz, zero_pow hj_ne]
  · rw [← Finset.mul_sum]
    have sum_eq : ∑ k ∈ Finset.range N, (ω ^ k) ^ j = 0 := by
      -- Convert (ω^k)^j to ω^(k*j) = (ω^j)^k
      conv_lhs => arg 2; ext k; rw [← pow_mul, mul_comm k j, pow_mul]
      have h_ne_one : ω ^ j ≠ 1 := by
        intro h_eq
        have h_div : N ∣ j := (hω.pow_eq_one_iff_dvd j).mp h_eq
        have : N ≤ j := Nat.le_of_dvd hj h_div
        omega
      have h_pow_N : (ω ^ j) ^ N = 1 := by
        calc (ω ^ j) ^ N = ω ^ (j * N) := by rw [← pow_mul]
          _ = ω ^ (N * j) := by rw [mul_comm]
          _ = (ω ^ N) ^ j := by rw [← pow_mul]
          _ = 1 ^ j := by rw [hω.pow_eq_one]
          _ = 1 := one_pow j
      have h_geom : (ω ^ j - 1) * ∑ k ∈ Finset.range N, (ω ^ j) ^ k = (ω ^ j) ^ N - 1 :=
        mul_geom_sum (ω ^ j) N
      rw [h_pow_N] at h_geom
      have h_eq_zero : (ω ^ j - 1) * ∑ k ∈ Finset.range N, (ω ^ j) ^ k = 0 := by
        rw [h_geom]; ring
      exact eq_zero_of_ne_zero_of_mul_left_eq_zero (sub_ne_zero_of_ne h_ne_one) h_eq_zero
    simp [sum_eq]