roots of unity - computing low degree

Stream: maths

[Christopher Hoskin]

In normal mathematics one might reasonably write that the complex cube roots of unity are the set [1]:

$$\{ 1, - 1 / 2 + √3 / 2i, - 1 / 2 - √3 / 2i\}$$

I'm trying to figure out how to formulate the equivalent statement in Mathlib. The (incomplete) best I've been able to come up with so far is:

open Complex

lemma complex_cubic_roots : Set.range ((fun x => Subtype.val x) ∘ ((rootsOfUnityEquivNthRoots ℂ 3 ∘
    (rootsOfUnityCircleEquiv 3 ∘ ZMod.rootsOfUnityAddChar 3)))) = {(1 : ℂ),
    (- (1/2) + √3 / 2 * I : ℂ), (- (1/2) - √3 / 2 * I : ℂ)} := by
  have e0 : ZMod.toCircle (N := 3) 0 = (1 : ℂ) := by
    aesop
  have e1 : ZMod.toCircle (N := 3) 1 = (- (1/2) + √3 / 2 * I : ℂ) := by
    have t0 : (1 : ZMod 3).val = 1 := rfl
    rw [ZMod.toCircle_eq_circleExp, t0, Lean.Grind.CommRing.natCast_one, ← mul_div_assoc, mul_one,
      Circle.coe_exp, ← Complex.cos_add_sin_I]
    norm_cast
    norm_num
  have e2 : ZMod.toCircle (N := 3) 2 = (- (1/2) - √3 / 2 * I : ℂ) := by
    have t0 : (2 : ZMod 3).val = 2 := rfl
    rw [ZMod.toCircle_eq_circleExp, t0, ← mul_div_assoc, mul_comm, ← mul_assoc, Circle.coe_exp,
      ← Complex.cos_add_sin_I]
    norm_cast
    norm_num
    aesop
  rw [ ← e1, ← e2]
  conv_rhs => rw [← e0]
  ext x
  simp only [Set.mem_insert_iff, Set.mem_singleton_iff]
  symm
  constructor <;> intro hx
  · simp
    rcases hx with h1 | ⟨h2|h3⟩
    · exact ⟨0, h1.symm⟩
    · exact ⟨1, h2.symm⟩
    · exact ⟨2, h3.symm⟩
  · simp at hx
    obtain ⟨y, hy⟩ := hx
    rw [← hy]
    rw [Subtype.val_inj, Subtype.val_inj, Subtype.val_inj]
    rw [ZMod.injective_toCircle.eq_iff, ZMod.injective_toCircle.eq_iff,
      ZMod.injective_toCircle.eq_iff]
    sorry

(This makes use of two PRs which are awaiting approval:

• #24416
• #25009
  )

The sorry is skipping the goal:

case h.mpr.intro
e0 : ↑(ZMod.toCircle 0) = 1
e1 : ↑(ZMod.toCircle 1) = -(1 / 2) + ↑√3 / 2 * I
e2 : ↑(ZMod.toCircle 2) = -(1 / 2) - ↑√3 / 2 * I
x : ℂ
y : ZMod 3
hy : ↑(ZMod.toCircle y) = x
⊢ y = 0 ∨ y = 1 ∨ y = 2

I suspect I'm going about this completely the wrong way?

[1] https://en.wikipedia.org/wiki/Root_of_unity#Explicit_expressions_in_low_degrees

[Kevin Buzzard]

How about

import Mathlib

open Complex

example : {z : ℂ | z^3 = 1} = {1, -1+I*√3, -1-I*√3} := by
  sorry

? That looks like a far more natural way to say it.

As for your sorry,

example (y : ZMod 3) : y = 0 ∨ y = 1 ∨ y = 2 := by fin_cases y <;> aesop

[Christopher Hoskin]

Thanks, fin_cases y <;> aesop closes the goal.

Kevin Buzzard said:

That looks like a far more natural way to say it.

Well yes - but I was thinking one would want to take advantage of the results in RingTheory/RootsOfUnity as well?

[Kevin Buzzard]

Sure but why would you want to load them all into the statement thus making the lemma unreadable and probably also unusable?

[Jz Pan]

Christopher Hoskin said:

but I was thinking one would want to take advantage of the results in RingTheory/RootsOfUnity as well?

Then add another lemma which state that RootsOfUnity when converted to Set R is {z : R | z^n = 1}.

[Christopher Hoskin]

I had another go and came up with this:

open Complex

set_option maxHeartbeats 500000 in
example : {z : ℂ | z^3 = 1} = {1, -(1/2)+√3/2*I, -(1/2)-√3/2*I} := by
  ext x
  simp only [Set.mem_setOf_eq, Set.mem_insert_iff, Set.mem_singleton_iff]
  rw [← Polynomial.mem_nthRoots zero_lt_three]
  let ω₀ := (rootsOfUnityCircleEquiv 3 ∘ ZMod.rootsOfUnityAddChar 3) 0
  let ω₁ := (rootsOfUnityCircleEquiv 3 ∘ ZMod.rootsOfUnityAddChar 3) 1
  let ω₂ := (rootsOfUnityCircleEquiv 3 ∘ ZMod.rootsOfUnityAddChar 3) 2
  have f1 : (Set.univ : Set (Fin 3)) = {0, 1, 2} := by
    ext y
    fin_cases y <;> aesop
  have r1 : (Set.univ : Set (rootsOfUnity 3 ℂ)) = {ω₀, ω₁, ω₂} := by
    rw [← Set.range_eq_univ.mpr (surjective_rootsOfUnityCircleEquiv_comp_rootsOfUnityAddChar 3),
      ← Set.image_univ, f1]
    ext x
    rw [Set.mem_image]
    constructor
    · intro hx
      obtain ⟨i, hi1, hi2⟩ := hx
      aesop
    · intro hx
      simp at hx
      rcases hx with h1 | ⟨h2|h3⟩
      · exact ⟨0, by aesop⟩
      · exact ⟨1, by aesop⟩
      · exact ⟨2, by aesop⟩
  have e0 : (rootsOfUnityEquivNthRoots ℂ 3) ω₀ = (1 : ℂ) := by
    aesop
  have e1 : (rootsOfUnityEquivNthRoots ℂ 3) ω₁ = (- (1/2) + √3 / 2 * I : ℂ) := by
    simp_all only [Function.comp_apply, AddChar.map_zero_eq_one, map_one,
      rootsOfUnityEquivNthRoots_apply, OneMemClass.coe_one, Units.val_one,
      rootsOfUnityCircleEquiv_apply, ZMod.rootsOfUnityAddChar_val, one_div, ω₁]
    obtain ⟨val_1, property_1⟩ := ω₁
    simp_all only [_root_.mem_rootsOfUnity]
    have t0 : (1 : ZMod 3).val = 1 := rfl
    rw [ZMod.toCircle_eq_circleExp, t0, Lean.Grind.CommRing.natCast_one, ← mul_div_assoc, mul_one,
      Circle.coe_exp, ← Complex.cos_add_sin_I]
    norm_cast
    norm_num
  have e2 : (rootsOfUnityEquivNthRoots ℂ 3) ω₂ = (- (1/2) - √3 / 2 * I : ℂ) := by
    simp_all only [Function.comp_apply, AddChar.map_zero_eq_one, map_one,
      rootsOfUnityEquivNthRoots_apply, OneMemClass.coe_one, Units.val_one,
      rootsOfUnityCircleEquiv_apply, ZMod.rootsOfUnityAddChar_val, one_div, ω₂]
    obtain ⟨val_2, property_2⟩ := ω₂
    simp_all only [_root_.mem_rootsOfUnity]
    have t0 : (2 : ZMod 3).val = 2 := rfl
    rw [ZMod.toCircle_eq_circleExp, t0, ← mul_div_assoc, mul_comm, ← mul_assoc, Circle.coe_exp,
      ← Complex.cos_add_sin_I]
    norm_cast
    norm_num
    aesop
  rw [← e2, ← e1]
  symm
  constructor
  · intro hx
    simp only [Nat.ofNat_pos, Polynomial.mem_nthRoots]
    rcases hx with h1 | ⟨h2|h3⟩
    · aesop
    · rw [h2]
      obtain ⟨val, property⟩ := ω₁
      simp only [rootsOfUnityEquivNthRoots_apply]
      exact (mem_rootsOfUnity' 3 val).mp property
    · rw [h3]
      obtain ⟨val, property⟩ := ω₂
      simp only [rootsOfUnityEquivNthRoots_apply]
      exact (mem_rootsOfUnity' 3 val).mp property
  · intro hx
    let ω := (rootsOfUnityEquivNthRoots ℂ 3).invFun ⟨x,hx⟩
    have ex : (rootsOfUnityEquivNthRoots ℂ 3) ω = x := by
      aesop
    have w1 : ω ∈ ({ω₀, ω₁, ω₂} : Set (rootsOfUnity 3 ℂ)) := by
      rw [← r1]
      trivial
    rw [← ex]
    rw [Subtype.val_inj, Subtype.val_inj]
    rw [(rootsOfUnityEquivNthRoots ℂ 3).injective.eq_iff,
      (rootsOfUnityEquivNthRoots ℂ 3).injective.eq_iff]
    have ee0 : ↑((rootsOfUnityEquivNthRoots ℂ 3) ω) = (1 : ℂ) ↔ ω = ω₀ := by
      constructor
      · intro hw
        have ee1 : (↑((rootsOfUnityEquivNthRoots ℂ 3) ω₀) : ℂ) =
            ↑((rootsOfUnityEquivNthRoots ℂ 3) ω) := by
          rw [hw, e0]
        rw [Subtype.val_inj, (rootsOfUnityEquivNthRoots ℂ 3).injective.eq_iff] at ee1
        rw [ee1]
      · intro hw
        rw [hw, e0]
    rw [ee0]
    exact w1

Doubtless still full of muddle :(

[Michael Stoll]

import Mathlib

open Complex

example : {z : ℂ | z^3 = 1} = {1, -(1/2)+√3/2*I, -(1/2)-√3/2*I} := by
  have H (z : ℂ) : z ^ 3 - 1 =
      (z - 1) * (z - (-(1 / 2) + √3 / 2 * I)) * (z - (-(1 / 2) - √3 / 2 * I)) := by
    ring_nf
    rw [I_sq, ← ofReal_pow, Real.sq_sqrt zero_le_three, ofReal_ofNat]
    ring
  ext1 z
  simp only [Set.mem_setOf_eq, Set.mem_insert_iff, Set.mem_singleton_iff]
  rw [← sub_eq_zero, H]
  simp only [mul_eq_zero, sub_eq_zero, or_assoc]

(golfed a bit more...)

[Kim Morrison]

The have could even be a separate lemma.

[Kim Morrison]

I really hope we can get this down to

have ... : ... := by grind
grind

later this year... :-) It doesn't seem implausible.

For now the best one can do is

example : {z : ℂ | z^3 = 1} = {1, -(1/2)+√3/2*I, -(1/2)-√3/2*I} := by
  have H (z : ℂ) : z ^ 3 - 1 =
      (z - 1) * (z - (-(1 / 2) + √3 / 2 * I)) * (z - (-(1 / 2) - √3 / 2 * I)) := by
    ring_nf
    rw [I_sq, ← ofReal_pow, Real.sq_sqrt (show 0 ≤ 3 by positivity)]
    push_cast
    ring
  ext1 z
  simp only [Set.mem_setOf_eq]
  grind +ring [sub_eq_zero, one_div, mul_eq_zero, Set.mem_insert_iff, Set.mem_singleton_iff]

It will take some special handling for grind to understand Set.mem_setOf_eq, but it is warranted. Grind doesn't understand fractions quite yet (soon!), so it can't do anything inside the have.

[Kevin Buzzard]

As ever I'll note that Hariharan and Macbeth have got a beefed-up norm_num which knows $i^2=-1$.

[Michael Stoll]

There's still ↑√3 ^ 2 = 3 to deal with (where the coercion is from Real to Complex)... (which norm_num doesn't seem to do right now)

[Kevin Buzzard]

Aah yes, @Sidharth Hariharan did you also deal with this? (sorry to ping you during exams)

[Michael Stoll]

actually,

import Mathlib

open Complex

example : I ^ 2 = -1 := by norm_num

does work (with current Mathlib in the web editor).

[Michael Stoll]

And, I guess, what we would like is that ring/ring_nf uses these features of norm_num as well, so that we can prove the have just by ring (or, in the future, by grind +ring I assume).

[Ruben Van de Velde]

I suspect that norm_num only works because its internal simp solves it

[Christopher Hoskin]

So it seems I was barking up the wrong tree with rootsOfUnity? Would it be fair to say that the material in RingTheory/RootsOfUnity is of theoretical interest, but not particularly useful for calculations?

[Michael Stoll]

Ruben Van de Velde said:

I suspect that norm_num only works because its internal simp solves it

Indeed, norm_num1 doesn't.

[Christopher Hoskin]

Jz Pan said:

Then add another lemma which state that RootsOfUnity when converted to Set R is {z : R | z^n = 1}.

Sorry, this is probably a rather obvious question, but how do I do that:

variable [CommMonoid R]

#check ((rootsOfUnity k R) : Set Rˣ)

#check ((rootsOfUnity k R) : Set R)

#check ((rootsOfUnity k R) : Set Rˣ) is okay, but #check ((rootsOfUnity k R) : Set R) gives me:

type mismatch
  rootsOfUnity k R
has type
  Subgroup Rˣ : Type u_4
but is expected to have type
  Set R : Type u_4

[Riccardo Brasca]

There is no coercion that does it automatically, but you can just take the image under the natural map (I am on mobile so writing the code is a little painful, sorry).

[Christopher Hoskin]

So, Units.val '' (rootsOfUnity k R) = {z : R | z^k = 1} ?

[Riccardo Brasca]

Yes, that would work. You can also use ((\up) : R \to ...) (It's probably the same thing, but it makes clear we are talking about the coercion)

[Riccardo Brasca]

It's even possible you don't need to specify the domain and the codomain

[Sidharth Hariharan]

Kevin Buzzard said:

Aah yes, Sidharth Hariharan did you also deal with this? (sorry to ping you during exams)

Ours was mainly designed to compute normal forms in $\mathbb{C}$ but this goes back to a discussion we had at London Learning Lean a few months back on whether similar ideas can be applied in other algebraic extensions (ie if norm_num solves numerical expressions in $K \subset L$ then it should be able to normalise expressions in $L$ by reducing them component-wise and applying whatever algebraic relations exist between elements of a $K$-basis of $L$)

[Edison Xie]

Sidharth Hariharan said:

Ours was mainly designed to compute normal forms in $\mathbb{C}$)

I think the tactic right now doesn’t handle this (just tried), but if we allow norm_num to work on simplify our normal form maybe our norm_numI will be able to prove this, this is another issue Heather and I discussed in Xena last week

[Sidharth Hariharan]

Michael Stoll said:

actually,

import Mathlib

open Complex

example : I ^ 2 = -1 := by norm_num

does work (with current Mathlib in the web editor).

I’m actually pleasantly surprised by this. Unfortunately more complicated expressions in $\mathbb{C}$ that only involve numerals and I aren’t currently handled by the existing norm_num. I’m still playing with it a little but ours can currently handle all sorts of things, including inversion and complex conjugation

[Kevin Buzzard]

Remember norm_num uses simp and it's simp which is solving this, not norm_num1

[Sidharth Hariharan]

Ah right

[Christopher Hoskin]

Riccardo Brasca said:

You can also use ((\up) : R \to ...) (It's probably the same thing, but it makes clear we are talking about the coercion)

Thanks. lemma test : ((↑) : Rˣ → _) '' (rootsOfUnity k R) = {z : R | z^k = 1} := sorry works. Replacing Rˣ with _ results in a "cannot coerce x to type R" error.

[Christopher Hoskin]

Jz Pan said:

Then add another lemma which state that RootsOfUnity when converted to Set R is {z : R | z^n = 1}.

I've done a PR for this lemma #25215

[Kenny Lau]

we have Zsqrtd d for $\Bbb Z[\sqrt d]$, but do we have $\Bbb Z\left[\frac{1+\sqrt{1+4k}}2\right]$? This seems useful.

[Kenny Lau]

if there's not, i don't mind doing it

[Kenny Lau]

I should probably ask this at #Is there code for X?

[Christopher Hoskin]

More generally this works in a Field where 2 is non-zero and -3 has a square root:

open Complex

lemma cyclotomic_polynomial_3_roots {K : Type*} [Field K] [NeZero (2 : K)] {z s : K}
    (hs : s * s = -3) : z ^ 2 + z + 1 = 0 ↔ z = -(1 / 2) + s / 2 ∨ z = -(1 / 2) - s / 2 := by
  suffices 1 * (z * z) + 1 * z + 1 = 0 ↔ z = -(1 / 2) + s / 2 ∨ z = -(1 / 2) - s / 2 by
    rw [← this]; ring_nf
  rw [quadratic_eq_zero_iff one_ne_zero (by rw [hs, discrim]; norm_num)]
  ring_nf

lemma cubic_roots_of_unity {K : Type*} [Field K] [NeZero (2 : K)] {s : K} (hs : s * s = -3) :
  {z : K | z^3 = 1} = {1, -(1 / 2) + s / 2, -(1 / 2) - s / 2} := by
  have H (z : K) : z ^ 3 - 1 = (z - 1) * (z ^ 2 + z + 1) := by ring
  ext1 z
  simp only [Set.mem_setOf_eq, Set.mem_insert_iff, Set.mem_singleton_iff]
  rw [← sub_eq_zero, H, ← cyclotomic_polynomial_3_roots hs, mul_eq_zero, sub_eq_zero]

example : {z : ℂ | z ^ 3 = 1} = {1, -(1 / 2) + √3 / 2 * I, -(1 / 2) - √3 / 2 * I} := by
  have hs : (√3 * I) * (√3 * I) = -3 := by
    ring_nf
    rw [I_sq, ← ofReal_pow, Real.sq_sqrt zero_le_three, mul_neg, mul_one,  ofReal_ofNat]
  rw [cubic_roots_of_unity hs]
  ring_nf

[Christopher Hoskin]

And, for completeness:

lemma cyclotomic_polynomial_3_ne_zero_of_sq_ne {K : Type*} [Field K] [NeZero (2 : K)] {z : K}
    (h : ∀ s : K, s^2 ≠ -3) : z ^ 2 + z + 1 ≠ 0 := by
  suffices 1 * (z * z) + 1 * z + 1 ≠ 0 by
    rw[one_mul, one_mul, ← sq] at this
    exact this
  exact quadratic_ne_zero_of_discrim_ne_sq (fun s => by rw [discrim]; ring; exact (h s).symm) _

lemma cubic_roots_of_unity_of_sq_ne {K : Type*} [Field K] [NeZero (2 : K)]
    (h : ∀ s : K, s^2 ≠ -3) : {z : K | z^3 = 1} = {1} := by
  have H (z : K) : z ^ 3 - 1 = (z - 1) * (z ^ 2 + z + 1) := by ring
  ext1 z
  simp only [Set.mem_setOf_eq, Set.mem_insert_iff, Set.mem_singleton_iff]
  rw [← sub_eq_zero, H, mul_eq_zero_iff_right (cyclotomic_polynomial_3_ne_zero_of_sq_ne h),
    sub_eq_zero]

example : {z : ℝ | z ^ 3 = 1} = {1} := by
  rw [cubic_roots_of_unity_of_sq_ne]
  intro s
  by_contra hc
  have e2 : s=0 := by
    rw [← sq_nonpos_iff]
    simp_all only [Left.neg_nonpos_iff, Nat.ofNat_nonneg]
  rw [e2, zero_pow two_ne_zero] at hc
  simp_all only [zero_eq_neg, OfNat.ofNat_ne_zero]

[Christopher Hoskin]

I'm now wondering what this looks like for polynomials. So far I have:

import Mathlib

open Polynomial

variable {R : Type*} [Field R]

variable {a b c : R}

/-- Roots of a quadratic equation. -/
theorem Polynomial.quadratic_eq_zero_iff [NeZero (2 : R)] (ha : a ≠ 0) {s : R}
    (h : discrim a b c = s * s) (x : R) :
    IsRoot (a • X ^ 2 + b • X + C c) x ↔ x = (-b + s) / (2 * a) ∨ x = (-b - s) / (2 * a) := by
  rw [IsRoot.def]
  simp only [eval_add, eval_smul, eval_pow, eval_X, smul_eq_mul, eval_C]
  rw [sq]
  rw [_root_.quadratic_eq_zero_iff ha h]

theorem Polynomial.quadratic_eq_zero_iff_of_discrim_eq_zero [NeZero (2 : R)] (ha : a ≠ 0)
    (h : discrim a b c = 0) (x : R) : IsRoot (a • X ^ 2 + b • X + C c) x ↔ x = -b / (2 * a) := by
  rw [IsRoot.def]
  simp only [eval_add, eval_smul, eval_pow, eval_X, smul_eq_mul, eval_C]
  rw [sq]
  rw [_root_.quadratic_eq_zero_iff_of_discrim_eq_zero ha h]

/-- A quadratic has no root if its discriminant has no square root. -/
theorem Polynomial.quadratic_ne_zero_of_discrim_ne_sq (h : ∀ s : R, discrim a b c ≠ s^2) (x : R) :
    ¬ IsRoot (a • X ^ 2 + b • X + C c) x := by
  rw [IsRoot.def]
  simp only [eval_add, eval_smul, eval_pow, eval_X, smul_eq_mul, eval_C]
  rw [← ne_eq, sq]
  exact _root_.quadratic_ne_zero_of_discrim_ne_sq h _

set_option linter.style.multiGoal false
theorem Polynomial.quadratic_roots_of_discrim_eq_zero [NeZero (2 : R)] [DecidableEq R] (ha : a ≠ 0)
    {z s : R} (h : discrim a b c = s * s) :
    z ∈ (a • X ^ 2 + b • X + C c).roots ↔ z = (-b + s) / (2 * a) ∨ z = (-b - s) / (2 * a) := by
  rw [Polynomial.mem_roots]
  rw [Polynomial.quadratic_eq_zero_iff ha h]
  have hc : (a • X ^ 2 + b • X + C c).coeff 2 = a := by
    simp [coeff_add, coeff_C_mul, coeff_smul, coeff_X_of_ne_one (Nat.add_one_add_one_ne_one),
      coeff_C_ne_zero (n:=2) ((Nat.zero_ne_add_one 1).symm)]
  rw [← hc] at ha
  by_contra hx
  have h : (a • X ^ 2 + b • X + C c).coeff 2 = (0 : R[X]).coeff 2 := by
    rw [hx]
  have h' : (a • X ^ 2 + b • X + C c).coeff 2 = 0 := by
    simp_all only [coeff_zero, ne_eq, not_true_eq_false]
  rw [hc] at h'
  exact ha h

theorem Polynomial.quadratic_roots [NeZero (2 : R)] [DecidableEq R] (ha : a ≠ 0) {s : R}
    (h : discrim a b c = s * s) :
    (a • X ^ 2 + b • X + C c).roots = {(-b + s) / (2 * a), (-b - s) / (2 * a)} := by
  sorry

Not sure where I'd go next (not very familiar with multisets).