Scaled Polynomial Constant Term Formula

  unfold scaledPolynomial constantTermFormula unscaledPolynomial
  rw [coeff_C_mul]
  -- Convert polynomialFromRealRoots to multiset form using list_foldr_eq_multiset_prod
  have h_eq : polynomialFromRealRoots (realProjectionsList N θ) =
      ((realProjectionsList N θ : Multiset ℝ).map (fun r => X - C r)).prod := by
    unfold polynomialFromRealRoots
    rw [list_foldr_eq_multiset_prod (realProjectionsList N θ)]
  -- Apply Vieta's formula: coefficient 0 is (-1)^N * esymm N
  have h_card : (realProjectionsList N θ : Multiset ℝ).card = N := by
    rw [Multiset.coe_card, card_realProjectionsList]
  rw [h_eq]
  rw [Multiset.prod_X_sub_C_coeff (realProjectionsList N θ : Multiset ℝ) (by rw [h_card]; omega : 0 ≤ (realProjectionsList N θ : Multiset ℝ).card)]
  -- Simplify card - 0 to card, then use h_card to get N
  simp only [Nat.sub_zero, h_card]
  -- esymm N of N elements is the product (when selecting all N elements, only one subset exists)
  have h_esymm : ((realProjectionsList N θ : Multiset ℝ).esymm N) = (realProjectionsList N θ).prod := by
    unfold Multiset.esymm
    let s := (realProjectionsList N θ : Multiset ℝ)
    -- powersetCard N s when card s = N has exactly one element
    have h_card_powerset : (s.powersetCard N).card = 1 := by
      rw [Multiset.card_powersetCard, h_card, Nat.choose_self]
    -- That one element is s itself
    have h_mem : s ∈ s.powersetCard N := by
      rw [Multiset.mem_powersetCard]
      exact ⟨le_refl _, h_card⟩
    -- A multiset with one element containing s must equal {s}
    have h_eq : s.powersetCard N = {s} := by
      rw [Multiset.card_eq_one] at h_card_powerset
      obtain ⟨t, ht⟩ := h_card_powerset
      rw [ht] at h_mem
      rw [Multiset.mem_singleton] at h_mem
      rw [ht, h_mem]
    rw [h_eq, Multiset.map_singleton, Multiset.sum_singleton]
    rfl
  rw [h_esymm]
  ring