Elementary Symmetric Functions of Rotated Roots Are Invariant

 (realProjectionsList N θ₁ : Multiset ℝ)
    let l2 := (realProjectionsList N θ₂ : Multiset ℝ)
    l1.esymm m = l2.esymm m := by
  intro l1 l2
  induction m using Nat.strong_induction_on with
  | h k IH =>
    cases k with
    | zero => omega
    | succ k' =>
      cases k' with
      | zero =>
        -- N ≥ 2, so sum of cosines at N equally spaced angles = 0
        -- esymm 1 is the sum, and both sums equal 0
        rw [multiset_esymm_one_eq_sum, multiset_esymm_one_eq_sum]
        simp only [l1, l2]
        have hN2 : 2 ≤ N := by omega
        rw [multiset_coe_realProjectionsList_sum, multiset_coe_realProjectionsList_sum]
        rw [sum_cos_roots_of_unity N θ₁ hN2, sum_cos_roots_of_unity N θ₂ hN2]
      | succ k'' =>
        -- Use Newton's identity to express esymm (k''+2) in terms of smaller esymm and power sums
        have h_newton_l1 := multiset_newton_identity l1 (k'' + 2) (by omega)
        have h_newton_l2 := multiset_newton_identity l2 (k'' + 2) (by omega)

        -- The RHS of Newton's identity is the same for l1 and l2
        have h_rhs_eq : (Finset.antidiagonal (k'' + 2)).sum (fun a =>
              if a.1 < k'' + 2 then (-1)^a.1 * l1.esymm a.1 * (l1.map (fun x => x ^ a.2)).sum
              else 0) =
            (Finset.antidiagonal (k'' + 2)).sum (fun a =>
              if a.1 < k'' + 2 then (-1)^a.1 * l2.esymm a.1 * (l2.map (fun x => x ^ a.2)).sum
              else 0) := by
          -- Each term in the sum is equal
          apply Finset.sum_congr rfl
          intro a ha_mem
          -- a ∈ antidiagonal (k'' + 2), so a.1 + a.2 = k'' + 2
          have h_sum : a.1 + a.2 = k'' + 2 := Finset.mem_antidiagonal.mp ha_mem
          split_ifs with ha
          · -- When a.1 < k'' + 2, we need to show the terms are equal
            congr 1
            · -- esymm a.1 is invariant by IH
              by_cases ha0 : a.1 = 0
              · simp only [ha0, pow_zero, one_mul]
                -- esymm 0 = 1 for any multiset
                simp only [Multiset.esymm, Multiset.powersetCard_zero_left]
              · have ha_pos : 0 < a.1 := Nat.pos_of_ne_zero ha0
                -- Need to show a.1 < N
                have ha_N : a.1 < N := by
                  calc a.1
                    _ < k'' + 2 := ha
                    _ = k'' + 1 + 1 := by ring
                    _ < N := hm'
                have h_esymm := IH a.1 ha ha_pos ha_N
                rw [h_esymm]
            · -- Power sum is invariant by powerSumCos_invariant
              simp only [l1, l2]
              rw [multiset_powersum_realProjectionsList, multiset_powersum_realProjectionsList]
              by_cases ha2_zero : a.2 = 0
              · simp [ha2_zero]
              · have ha2_pos : 0 < a.2 := Nat.pos_of_ne_zero ha2_zero
                -- Need to show a.2 < N
                have ha2_N : a.2 < N := by
                  -- From antidiagonal: a.1 + a.2 = k'' + 2 and a.1 < k'' + 2
                  -- So a.2 = k'' + 2 - a.1 > 0, and since k'' + 2 < N, we have a.2 < N
                  calc a.2
                    _ = k'' + 2 - a.1 := by omega
                    _ ≤ k'' + 2 := by omega
                    _ = k'' + 1 + 1 := by ring
                    _ < N := hm'
                exact powerSumCos_invariant N a.2 θ₁ θ₂ hN ha2_pos ha2_N
          · rfl

        -- From Newton's identity: m * esymm m = (-1)^(m+1) * RHS
        -- So if RHS are equal, then m * esymm m are equal
        have h_prod_eq : ((k'' + 2) : ℝ) * l1.esymm (k'' + 2) =
            ((k'' + 2) : ℝ) * l2.esymm (k'' + 2) := by
          have h1 := h_newton_l1
          have h2 := h_newton_l2
          push_cast at h1 h2
          rw [h1, h2, h_rhs_eq]

        -- Divide by (k'' + 2) to get the result
        have h_nonzero : (↑k'' + 2 : ℝ) ≠ 0 := by positivity
        exact mul_left_cancel₀ h_nonzero h_prod_eq