Compute roots of polynomials

Stream: new members

[Ludwig Monnerjahn]

I want to show that $h(x) := x^3 - 6/8 x - 1/8$ is irreducible, there I showed that it is sufficient to show that $h$ has no roots. I thought of proving this using the rational roots theorem for $8\cdot h$.
But I cannot figure out how to convert the goal to $h(x) = 0 \leftrightarrow x | 8$ (comment 1) and then $x \in \{1, 1/2, 1/4, 1/8\}$ (comment 2).

Does anyone know how to write the proof according to comments 1-4 or has a better way of structuring the proof?

import mathlib

open Polynomial

noncomputable def H : Polynomial ℚ := X ^ 3 - (C (6/8) * X) - C (1/8)  -- x^3 - 6/8 x - 1/8
noncomputable def H' : Polynomial ℚ := C 8 *X ^ 3 - (C 6 * X) - C 1 -- 8x^3 - 6 x - 1

lemma H_eq_mul_H' : H = C (1/8) * H' := by sorry

lemma H_roots: roots H = 0:= by
  have heq: roots H = roots H' := by sorry
  rw[heq]
  refine Multiset.eq_zero_of_forall_not_mem ?_
  by_contra h
  simp at h
  obtain ⟨a, x, h⟩ := h
  have h₀: ↑(IsFractionRing.den ℚ x) ∣ Polynomial.leadingCoeff H' := by apply den_dvd_of_is_root h
  have h₁: Polynomial.leadingCoeff H' = 8 := by sorry

  --! 1. show H(x) = 0 ↔ x | 8
  --! 2. show x ∈ {1, 1/2, 1/4, 1/8}

  -- 3. show H(z) ≠ 0 for z ∈ {1, 1/2, 1/4, 1/8}
  -- 4. contradiction
  sorry

[Damiano Testa]

I suspect that you may find docs#num_dvd_of_is_root useful. I also think that there is probably the need of more API lemmas around this to actually make proving your lemma easier.

[Ludwig Monnerjahn]

Yay, maybe the constant coefficient is easier, I tried using docs#den_dvd_of_is_root. But it the problems I get stuct remain the same.

[Riccardo Brasca]

Nothing serious, but your example doesn't work (you wrote mathlib instead of Mathlib in the first line).

[Riccardo Brasca]

Anyway let me have a look. My first suggestion is to take H' := IsLocalization.integerNormalization _ H

[Damiano Testa]

I played with this a little, but I have to stop for a bit. I made some progress, but did not manage to finish it.

noncomputable abbrev H : Polynomial ℚ := X ^ 3 - (C (6/8) * X) - C (1/8)  -- x^3 - 6/8 x - 1/8
noncomputable abbrev H' : Polynomial ℚ := C 8 *X ^ 3 - (C 6 * X) - C 1 -- 8x^3 - 6 x - 1

noncomputable abbrev Hz : Polynomial ℤ := C 8 *X ^ 3 - (C 6 * X) - C 1 -- 8x^3 - 6 x - 1

-- this statement is ready for use in an extension of `compute_degree` that should be able to solve `leadingCoeff` goals.
lemma leadingCoeff_eq_of_natDegree_le_of_coeff_eq {R} [Semiring R] {r : R} (r0 : r ≠ 0) {f : R[X]}
    {n : Nat} (dn : f.natDegree ≤ n) (cn : f.coeff n = r) :
    f.leadingCoeff = r := by
  subst cn
  unfold leadingCoeff
  rw [natDegree_eq_of_le_of_coeff_ne_zero dn r0]

theorem leadingCoeffHz : leadingCoeff Hz = 8 := by
  unfold Hz
  apply leadingCoeff_eq_of_natDegree_le_of_coeff_eq (n := 3)
  · compute_degree
  · simp only [coeff_sub, coeff_C_mul, coeff_X_pow, ↓reduceIte, mul_one, coeff_X, coeff_C]
    norm_num
  · norm_num

theorem leadingCoeffH'z : leadingCoeff H' = 8 := by
  unfold H'
  apply leadingCoeff_eq_of_natDegree_le_of_coeff_eq (n := 3)
  · compute_degree
  · simp only [coeff_sub, coeff_C_mul, coeff_X_pow, ↓reduceIte, mul_one, coeff_X, coeff_C]
    norm_num
  · norm_num

lemma H_roots: roots H = 0:= by
  have heq: roots H = roots H' := by sorry
  rw[heq]
  refine Multiset.eq_zero_of_forall_not_mem ?_
  by_contra h
  simp at h
  obtain ⟨a, x, h⟩ := h
  have h₁:= leadingCoeffH'z
  have := den_dvd_of_is_root (p := Hz) (r := x) ?_
  · have lcZ := leadingCoeffHz
    rw [lcZ, (show (8 : ℤ) = 2 ^ 3 from rfl), dvd_prime_pow (Int.prime_two)] at this
    rcases this with ⟨i, hi, h⟩
    have := num_dvd_of_is_root (p := Hz) (r := x) (by simp_all)
    simp [← isUnit_iff_dvd_one, Int.isUnit_iff] at this
    sorry  -- here the numerator and denominator are bounded, but there is some missing API around
           -- `IsFractionRing.num` and `IsFractionRing.den`
  · convert h
    simp
  done

[Riccardo Brasca]

I am also playing with it, and I dont think anymore that using IsLocalization.integerNormalization is a good idea.

[Riccardo Brasca]

From here it should be very easy to finish the list of roots

import Mathlib

open Polynomial nonZeroDivisors IsFractionRing

noncomputable def H : Polynomial ℚ := X ^ 3 - (C (6/8) * X) - C (1/8)
noncomputable def H' : Polynomial ℤ := C 8 * X ^ 3 - C 6 * X - 1

lemma H_H' : (C 8) * H = map (algebraMap ℤ ℚ) H' := by
  rw [H, H']
  simp only [one_div, smul_eq_C_mul, map_ofNat, mul_sub, algebraMap_int_eq, Polynomial.map_sub,
    Polynomial.map_mul, Polynomial.map_ofNat, Polynomial.map_pow, map_X, Polynomial.map_one]
  congr 2
  · rw [← mul_assoc]
    congr
    suffices : (8 : ℚ) * (6 / 8) = 6
    convert congr_arg C this
    simp only [map_mul, map_ofNat, mul_eq_mul_left_iff, OfNat.ofNat_ne_zero, or_false]
    norm_num
  · rw [show (8 : ℚ[X]) = C 8 by simp, ← map_mul]
    simp

lemma H'_degree : H'.natDegree = 3 := by
  rw [H']
  compute_degree
  repeat norm_num

lemma H_degree : H.natDegree = 3 := by
  rw [H]
  compute_degree
  repeat norm_num

lemma H_ne_zero : H ≠ 0 := by
  intro h
  have := H_degree
  simp [h] at this

lemma H'_coeff_zero : H'.coeff 0 = -1 := by simp [H']

lemma H'_leading : H'.leadingCoeff = 8 := by
  rw [leadingCoeff, H'_degree, H']
  simp [coeff_one, coeff_X]

lemma H'_roots_num {a : ℚ} (ha : aeval a H' = 0) : num ℤ a ∣ -1 := by
  convert num_dvd_of_is_root ha
  exact H'_coeff_zero.symm

lemma H'_roots_den {a : ℚ} (ha : aeval a H' = 0) : (den ℤ a : ℤ) ∣ 8 := by
  convert den_dvd_of_is_root ha
  exact H'_leading.symm

lemma roots_H_H' {a : ℚ} (ha : a ∈ H.roots) : aeval a H' = 0 := by
  have : aeval a ((C 8) * H) = 0 := by
    rw [mem_roots H_ne_zero, IsRoot.def] at ha
    simp [ha]
  rwa [H_H', aeval_map_algebraMap] at this

lemma roots_num {a : ℚ} (ha : a ∈ H.roots) : num ℤ a ∣ -1 :=
  H'_roots_num (roots_H_H' ha)

lemma roots_den {a : ℚ} (ha : a ∈ H.roots) : (den ℤ a : ℤ) ∣ 8 :=
  H'_roots_den (roots_H_H' ha)

[Riccardo Brasca]

I have to stop half an hour, but I will try to finish it later. It is anyway a ridiculous pain.

[Damiano Testa]

Riccardo, compute_degree! ~ compute_degree; repeat norm_num.

[Riccardo Brasca]

Do you see a way of making the first lemma immediate? I hoped for something like a ext split simp combo, but it didn't work

[Damiano Testa]

@Riccardo Brasca not really immediate:

lemma H_H' : (C 8) * H = map (algebraMap ℤ ℚ) H' := by
  rw [H, H']
  ext n : 1
  simp [coeff_X, coeff_C, coeff_one]
  split_ifs <;> norm_num

[Riccardo Brasca]

Ah, very nice! I didn't notice ext was going too deep.

[Damiano Testa]

There should be a lot more automation for dealing with explicit polynomials. Mathlib has "general" results, but these kinds of explicit computations appear quite a bit!

[Damiano Testa]

Yes, ext entering Rat is always an issue, in my experience.

[Damiano Testa]

Anyway, going from IsFractionRing.num/den to Rat.num/den is also painful, btw...

[Damiano Testa]

I think that we both got to this point and then had to do something else! :lol:

[Yaël Dillies]

Damiano Testa said:

Yes, ext entering Rat is always an issue, in my experience.

That can be fixed by removing the ext attribute from the relevant lemma. We did the same thing for ℂ recently.

[Riccardo Brasca]

OK, I have a proof of

lemma roots_num {a : ℚ} (ha : a ∈ H.roots) : num ℤ a = 1 ∨ num ℤ a = -1

and

lemma roots_den {a : ℚ} (ha : a ∈ H.roots) :
    (den ℤ a : ℤ) ∈ ({1, -1, 2, -2, 4, -4, 8, -8} : Finset ℤ)

[Riccardo Brasca]

But this still leaves 16 (of course 8 in reality) cases, very annoying

[Riccardo Brasca]

Well, not so much.

import Mathlib

open Polynomial IsFractionRing

noncomputable def H : Polynomial ℚ := X ^ 3 - (C (6/8) * X) - C (1/8)
noncomputable def H' : Polynomial ℤ := C 8 * X ^ 3 - C 6 * X - 1

lemma H_H' : (C 8) * H = map (algebraMap ℤ ℚ) H' := by
  rw [H, H']
  ext n : 1
  simp only [map_ofNat, one_div, coeff_ofNat_mul, coeff_sub, coeff_X_pow, coeff_C_mul, coeff_X,
    mul_ite, mul_one, mul_zero, coeff_C, algebraMap_int_eq, Polynomial.map_sub, Polynomial.map_mul,
    Polynomial.map_ofNat, Polynomial.map_pow, map_X, Polynomial.map_one, coeff_one]
  split_ifs <;> norm_num

lemma H'_degree : H'.natDegree = 3 := by
  rw [H']
  compute_degree!

lemma H_degree : H.natDegree = 3 := by
  rw [H]
  compute_degree!

lemma H_ne_zero : H ≠ 0 := fun h ↦ by simpa [h] using H_degree

lemma H'_coeff_zero : H'.coeff 0 = -1 := by simp [H']

lemma H'_leading : H'.leadingCoeff = 8 := by
  rw [leadingCoeff, H'_degree, H']
  simp [coeff_one, coeff_X]

lemma H'_roots_num {a : ℚ} (ha : aeval a H' = 0) : num ℤ a ∣ -1 := by
  convert num_dvd_of_is_root ha
  exact H'_coeff_zero.symm

lemma H'_roots_den {a : ℚ} (ha : aeval a H' = 0) : (den ℤ a : ℤ) ∣ 8 := by
  convert den_dvd_of_is_root ha
  exact H'_leading.symm

lemma roots_H_H' {a : ℚ} (ha : a ∈ H.roots) : aeval a H' = 0 := by
  have : aeval a ((C 8) * H) = 0 := by
    rw [mem_roots H_ne_zero, IsRoot.def] at ha
    simp [ha]
  rwa [H_H', aeval_map_algebraMap] at this

lemma roots_num {a : ℚ} (ha : a ∈ H.roots) : num ℤ a = 1 ∨ num ℤ a = -1 := by
  have := H'_roots_num (roots_H_H' ha)
  rwa [← Int.natAbs_dvd_natAbs, Int.natAbs_neg, Int.natAbs_one, Nat.dvd_one,
    Int.natAbs_eq_iff, Nat.cast_one] at this

lemma roots_den_abs {a : ℚ} (ha : a ∈ H.roots) :
    (den ℤ a : ℤ).natAbs ∈ ({1, 2, 4, 8} : Finset ℕ) := by
  rw [show {1, 2, 4, 8} = Nat.divisors 8 from rfl, Nat.mem_divisors, ← Int.coe_nat_dvd]
  exact ⟨by simp [H'_roots_den (roots_H_H' ha)], by norm_num⟩

lemma roots_den {a : ℚ} (ha : a ∈ H.roots) :
    (den ℤ a : ℤ) ∈ ({1, -1, 2, -2, 4, -4, 8, -8} : Finset ℤ) := by
  have := roots_den_abs ha
  simp only [Finset.mem_insert, Finset.mem_singleton] at this --can I use `fin_cases`?
  rcases this with (h | h | h | h) <;> rcases Int.natAbs_eq_iff.1 h with (h | h) <;> simp [h]

lemma roots : ¬ ∃ a, a ∈ H.roots := by
  intro ⟨a, ha⟩
  have := roots_den ha
  simp only [Int.reduceNeg, Finset.mem_insert, OneMemClass.coe_eq_one, Finset.mem_singleton] at this
  have num := roots_num ha
  rw [mem_roots H_ne_zero, IsRoot.def, H, eval_sub, eval_sub, eval_pow, eval_X, eval_mul, eval_C,
    eval_X, eval_C, ← mk'_num_den' ℤ a] at ha
  simp only [algebraMap_int_eq, eq_intCast, div_pow, one_div] at ha
  rcases num with (h | h) <;>
  rcases this with (h1 | h1 | h1 | h1 | h1 | h1 | h1 | h1) <;> norm_num [h, h1] at ha

[Riccardo Brasca]

@Ludwig Monnerjahn is the fact the the non existence of a root implies irreducibility in mathlb? I thought we were missing that.

[Riccardo Brasca]

Ah, it's #9697

[Riccardo Brasca]

Here a full proof that H is irreducible.

import Mathlib

open Polynomial IsFractionRing

section

theorem irreducible_iff_lt_natDegree_lt {R} [CommSemiring R] [NoZeroDivisors R]
    {p : R[X]} (hp : p.Monic) (hp1 : p ≠ 1) :
    Irreducible p ↔ ∀ q, Monic q → natDegree q ∈ Finset.Ioc 0 (natDegree p / 2) → ¬ q ∣ p := by
  rw [hp.irreducible_iff_natDegree', and_iff_right hp1]
  constructor
  · rintro h g hg hdg ⟨f, rfl⟩
    exact h f g (hg.of_mul_monic_left hp) hg (mul_comm f g) hdg
  · rintro h f g - hg rfl hdg
    exact h g hg hdg (dvd_mul_left g f)

theorem irreducible_iff_roots_eq_zero_of_degree_le_three {R} [CommRing R] [IsDomain R]
    {p : R[X]} (hp : p.Monic)
    (hp2 : 2 ≤ p.natDegree) (hp3 : p.natDegree ≤ 3) : Irreducible p ↔ p.roots = 0 := by
  rw [irreducible_iff_lt_natDegree_lt hp]; swap
  · rintro rfl; rw [natDegree_one] at hp2; cases hp2
  have hp0 : p ≠ 0 := by rintro rfl; rw [natDegree_zero] at hp2; cases hp2
  simp_rw [show p.natDegree / 2 = 1 from (Nat.div_le_div_right hp3).antisymm
    (by apply Nat.div_le_div_right (c := 2) hp2), show Finset.Ioc 0 1 = {1} from rfl,
    Finset.mem_singleton, Multiset.eq_zero_iff_forall_not_mem, mem_roots hp0, ← dvd_iff_isRoot]
  refine ⟨fun h r ↦ h _ (monic_X_sub_C r) (natDegree_X_sub_C r), fun h q hq hq1 ↦ ?_⟩
  rw [hq.eq_X_add_C hq1, ← sub_neg_eq_add, ← C_neg]; apply h

end

noncomputable def H : Polynomial ℚ := X ^ 3 - (C (6/8) * X) - C (1/8)
noncomputable def H' : Polynomial ℤ := C 8 * X ^ 3 - C 6 * X - 1

lemma H_H' : (C 8) * H = map (algebraMap ℤ ℚ) H' := by
  rw [H, H']
  ext n : 1
  simp only [map_ofNat, one_div, coeff_ofNat_mul, coeff_sub, coeff_X_pow, coeff_C_mul, coeff_X,
    mul_ite, mul_one, mul_zero, coeff_C, algebraMap_int_eq, Polynomial.map_sub, Polynomial.map_mul,
    Polynomial.map_ofNat, Polynomial.map_pow, map_X, Polynomial.map_one, coeff_one]
  split_ifs <;> norm_num

lemma H'_degree : H'.natDegree = 3 := by
  rw [H']
  compute_degree!

lemma H_degree : H.natDegree = 3 := by
  rw [H]
  compute_degree!

lemma H_monic : H.Monic := by
  unfold H
  monicity!

lemma H_ne_zero : H ≠ 0 := fun h ↦ by simpa [h] using H_degree

lemma H'_coeff_zero : H'.coeff 0 = -1 := by simp [H']

lemma H'_leading : H'.leadingCoeff = 8 := by
  rw [leadingCoeff, H'_degree, H']
  simp [coeff_one, coeff_X]

lemma H'_roots_num {a : ℚ} (ha : aeval a H' = 0) : num ℤ a ∣ -1 := by
  convert num_dvd_of_is_root ha
  exact H'_coeff_zero.symm

lemma H'_roots_den {a : ℚ} (ha : aeval a H' = 0) : (den ℤ a : ℤ) ∣ 8 := by
  convert den_dvd_of_is_root ha
  exact H'_leading.symm

lemma roots_H_H' {a : ℚ} (ha : a ∈ H.roots) : aeval a H' = 0 := by
  have : aeval a ((C 8) * H) = 0 := by
    rw [mem_roots H_ne_zero, IsRoot.def] at ha
    simp [ha]
  rwa [H_H', aeval_map_algebraMap] at this

lemma roots_num {a : ℚ} (ha : a ∈ H.roots) : num ℤ a = 1 ∨ num ℤ a = -1 := by
  have := H'_roots_num (roots_H_H' ha)
  rwa [← Int.natAbs_dvd_natAbs, Int.natAbs_neg, Int.natAbs_one, Nat.dvd_one,
    Int.natAbs_eq_iff, Nat.cast_one] at this

lemma roots_den_abs {a : ℚ} (ha : a ∈ H.roots) :
    (den ℤ a : ℤ).natAbs ∈ ({1, 2, 4, 8} : Finset ℕ) := by
  rw [show {1, 2, 4, 8} = Nat.divisors 8 from rfl, Nat.mem_divisors, ← Int.coe_nat_dvd]
  exact ⟨by simp [H'_roots_den (roots_H_H' ha)], by norm_num⟩

lemma roots_den {a : ℚ} (ha : a ∈ H.roots) :
    (den ℤ a : ℤ) ∈ ({1, -1, 2, -2, 4, -4, 8, -8} : Finset ℤ) := by
  have := roots_den_abs ha
  simp only [Finset.mem_insert, Finset.mem_singleton] at this --can I use `fin_cases`?
  rcases this with (h | h | h | h) <;> rcases Int.natAbs_eq_iff.1 h with (h | h) <;> simp [h]

lemma H_roots : ¬ ∃ a, a ∈ H.roots := by
  intro ⟨a, ha⟩
  have := roots_den ha
  simp only [Int.reduceNeg, Finset.mem_insert, OneMemClass.coe_eq_one, Finset.mem_singleton] at this
  have num := roots_num ha
  rw [mem_roots H_ne_zero, IsRoot.def, H, eval_sub, eval_sub, eval_pow, eval_X, eval_mul, eval_C,
    eval_X, eval_C, ← mk'_num_den' ℤ a] at ha
  simp only [algebraMap_int_eq, eq_intCast, div_pow, one_div] at ha
  rcases num with (h | h) <;>
  rcases this with (h1 | h1 | h1 | h1 | h1 | h1 | h1 | h1) <;> norm_num [h, h1] at ha

lemma H_irr : Irreducible H := by
  refine (irreducible_iff_roots_eq_zero_of_degree_le_three H_monic ?_ ?_).2
    (Multiset.eq_zero_of_forall_not_mem (fun a ha ↦ H_roots ⟨a, ha⟩)) <;>
  rw [H_degree]; norm_num

[Damiano Testa]

The lemmas leading up to the computation of the leading coefficients should really have to disappear: the tactic compute_degree essentially already "knows" that the leading coefficient is 8 and that it is non-zero (via norm_num).

I wonder if there is enough interest in actually developing further compute_degree into a polynomial tactic that handles most goals involving explicit polynomials.

[Damiano Testa]

@Riccardo Brasca

lemma H_monic : H.Monic := by
  unfold H
  monicity!

[Riccardo Brasca]

Damiano Testa said:

The lemmas leading up to the computation of the leading coefficients should really have to disappear: the tactic compute_degree essentially already "knows" that the leading coefficient is 8 and that it is non-zero (via norm_num).

I wonder if there is enough interest in actually developing further compute_degree into a polynomial tactic that handles most goals involving explicit polynomials.

I was thinking the same. In my opinion we really want that. This is the usual problem that polynomials are noncomputable, so all these computations are really nontrivial.

[Riccardo Brasca]

monicity?!

[Damiano Testa]

Indeed!

[Riccardo Brasca]

Those tactics are really proliferating

[Damiano Testa]

compute_degree deals with finding leading coefficients and proving inequalities of degrees. Once you have that, you can answer a lot of questions about polynomials. This may take some time, but I can look into automating many of the computations on explicit polynomials into a single tactic.

[Riccardo Brasca]

I think @Antoine Chambert-Loir has some thoughts here. In principle with computable polynomials everything would be much easier, right? (Of course I don't have any idea on how to do it properly)

[Damiano Testa]

I also liked @Yakov Pechersky experiment of encoding polynomials as Lists of coefficients. That would also simplify a lot of these computations.

[Damiano Testa]

Regarding computable polynomials: that would definitely help. In that case, maybe all these tactics would convert into norm_num extensions.

[Riccardo Brasca]

@Cyril Cohen how easy is to do this in Coq? We have the explicit polynomial, with rational coefficients, X ^ 3 - (C (6/8) * X) - C (1/8) (C a is the constant polynomial a) and we want to prove that it has no rational roots using the rational root test. We have the result in mathlib, but applying it in practice is a pain.

[Riccardo Brasca]

(BTW we want a version for the rational root test specific to ℤ and ℚ, without any localization stuff, that takes numerator and denominator as in the Rat API, at least this would divide the number of cases by two)

[Damiano Testa]

And I think that these lemmas

lemma H'_degree : H'.natDegree = 3 :=
lemma H_degree : H.natDegree = 3 :=
lemma H_monic : H.Monic :=
lemma H_ne_zero : H ≠ 0 :=
lemma H'_coeff_zero : H'.coeff 0 = -1 :=
lemma H'_leading : H'.leadingCoeff = 8 :=

should all be hidden by automation (and are probably very close to being invisible already, though not quite).

[Riccardo Brasca]

Making H and H' notation would surely help.

[Damiano Testa]

Or abbrev at least.

[Damiano Testa]

Yaël Dillies said:

Damiano Testa said:

Yes, ext entering Rat is always an issue, in my experience.

That can be fixed by removing the ext attribute from the relevant lemma. We did the same thing for ℂ recently.

What do people think about this? I think that I never wanted ext to break a rational number into its numerator and its denominator. Should I remove the instanceext attribute?

[Riccardo Brasca]

Let's try and see what happens

[Yakov Pechersky]

Funny that I'm mentioned in this thread, because iirc I'm also the one that added the ext attribute to Rat

[Damiano Testa]

The ext attribute is in Std.

[Damiano Testa]

docs#Rat.ext

[Riccardo Brasca]

ah, then I don't know

[Damiano Testa]

Yakov Pechersky said:

Funny that I'm mentioned in this thread, because iirc I'm also the one that added the ext attribute to Rat

You were really meant to chime in here! :smile:

[Damiano Testa]

Ok, I won't push for removing the ext attribute from Rat: it is easy enough to circumvent the issue by using ext : whaterver.

[Yakov Pechersky]

I think it should be removed

[Damiano Testa]

Ok, I have never made a PR to Std and this seems like a good time to learn!

[Yakov Pechersky]

Regarding list polys, if the construction is built up, then the tactic framework would be to traverse the expression and construct the list. Then construct the proof that the list.to_poly is equal to the original. And then hopefully the obvious lemmas either compute directly or via norm_num like proof construction via VM calculation

[Damiano Testa]

The expression traversal already happens for compute_degree, I think that it would be relatively easy meta-extract the "correct" list of coefficients. After that a little of ext; simp should prove that the list representation is accurate.

[Damiano Testa]

If I understand correctly, you would then need a few lemmas to connect leadingCoeff, degree, natDegree, Monic,... to the corresponding List properties, right?

[Yakov Pechersky]

Yes, exactly. And the traversal being in the existing tactic is a boon!

[Damiano Testa]

I will focus on extending the "old" approach for now, but I suspect that the list perspective is better in the long run.

[Ludwig Monnerjahn]

@Riccardo Brasca Thank you so much. I was stuck on the proof for days.

[Riccardo Brasca]

If you look at the proof I didn't do anything fancy, I just split the work in several small pieces.

[Riccardo Brasca]

BTW, is there a rcases syntax allowing to do these lines

rcases num with (h | h) <;>
rcases this with (h1 | h1 | h1 | h1 | h1 | h1 | h1 | h1)

in one?

[Damiano Testa]

Defining a "minimal effort" polynomial tactic makes it possible to inline most of the "obvious" lemmas above:

-- using `notation` to avoid having to `unfold` all the time.
notation "H" => (X ^ 3 - (C (6/8) * X) - C (1/8) : Polynomial ℚ)
notation "H'" => (C 8 * X ^ 3 - C 6 * X - 1 : Polynomial ℤ)

lemma H'_degree : natDegree H' = 3 := by polynomial!
lemma H_degree : natDegree H = 3 := by polynomial!
lemma H_monic : Monic H := by polynomial!
lemma H'_leading : leadingCoeff H' = 8 := by polynomial!
lemma H_ne_zero : H ≠ 0 := by
  apply_fun leadingCoeff; exact ne_of_eq_of_ne (by polynomial!) (by norm_num)

[Damiano Testa]

#10387: feel free to comment and add feature requests!

[Miguel Marco]

I don't know if this really help in this case, but it is related: there are factorization algorithms formally verified in Isabelle.

Maybe it would be worth trying to port them to Lean (although I guess that would be a big task). That could provide both a tactic that factorizes polynomials over the rationals, and a proof of the correcteness of the result.

[Junyan Xu]

I just did an experimental implementation that allows proving (C 8) * H = map (algebraMap ℤ ℚ) H' almost by reflection:

import Mathlib.Data.Polynomial.Eval

open scoped Polynomial
open Polynomial (C X)

variable {R} [Semiring R] (p q : R[X])

class Polynomial.ComputableRepr where
  coeff : ℕ → R
  degreeBound : ℕ
  coeff_eq ⦃n : ℕ⦄ : p.coeff n = coeff n
  support_subset : p.support ⊆ Finset.range degreeBound

open Polynomial.ComputableRepr

instance : (0 : R[X]).ComputableRepr where
  coeff _ := 0
  degreeBound := 0
  coeff_eq := by simp
  support_subset := by simp

-- Can't make this work:
/- noncomputable instance (m : ℕ) [m.AtLeastTwo] : (OfNat.ofNat m : R[X]).ComputableRepr where
  coeff n := if n = 0 then m else 0
  degreeBound := 1
  coeff_eq := by sorry
  support_subset := sorry
example : (2 : ℤ[X]).ComputableRepr := inferInstance -- fails -/

instance (r : R) : (C r).ComputableRepr where
  coeff n := if n = 0 then r else 0
  degreeBound := 1
  coeff_eq := by simp [Polynomial.coeff_C]
  support_subset := by rw [← Polynomial.monomial_zero_left]; apply Polynomial.support_monomial'

instance : (X : R[X]).ComputableRepr where
  coeff n := if n = 1 then 1 else 0
  degreeBound := 2
  coeff_eq _ := by simp_rw [Polynomial.coeff_X, eq_comm]
  support_subset := by
    rw [← Polynomial.monomial_one_one_eq_X]
    exact (Polynomial.support_monomial' _ _).trans (by simp)

variable [p.ComputableRepr]

-- TODO: replace assumptions of `support_smul` and here by `[SMulZeroClass S R]`
instance {S} [Monoid S] [DistribMulAction S R] (s : S) : (s • p).ComputableRepr where
  coeff n := s • coeff p n
  degreeBound := degreeBound p
  coeff_eq _ := by rw [Polynomial.coeff_smul, coeff_eq]
  support_subset := (Polynomial.support_smul _ _).trans support_subset

namespace Polynomial.ComputableRepr

variable {p} in
theorem coeff_eq_zero_of_le {n : ℕ} (hn : degreeBound p ≤ n) : p.coeff n = 0 :=
  Polynomial.not_mem_support_iff.mp fun h ↦ hn.not_lt (Finset.mem_range.mp <| support_subset h)

@[reducible] def ofEq (h : p = q) : q.ComputableRepr where
  coeff := coeff p
  degreeBound := degreeBound p
  coeff_eq _ := by rw [← coeff_eq, h]
  support_subset := h ▸ support_subset

def degreeAux (P : ℕ → Prop) [DecidablePred P] (n : ℕ) : WithBot ℕ :=
  if n = 0 then ⊥ else if P n then n else degreeAux P (n - 1)

variable [DecidableEq R]

def degree : WithBot ℕ := degreeAux (coeff p · ≠ 0) (degreeBound p - 1)
def natDegree : ℕ := (degree p).unbot' 0
def leadingCoeff : R := coeff p (natDegree p)
abbrev Monic : Prop := leadingCoeff p = 1

lemma degree_eq : p.degree = degree p := sorry
lemma natDegree_eq : p.natDegree = natDegree p := by rw [natDegree, ← degree_eq]; rfl
lemma leadingCoeff_eq : p.leadingCoeff = leadingCoeff p := by
  rw [leadingCoeff, ← coeff_eq, ← natDegree_eq]; rfl
lemma monic_iff : p.Monic ↔ Monic p := by rw [Monic, ← leadingCoeff_eq]; rfl

instance : Decidable p.Monic := decidable_of_iff _ (monic_iff p).symm

end Polynomial.ComputableRepr

noncomputable instance : (1 : R[X]).ComputableRepr := .ofEq (C 1) _ (map_one C)
noncomputable instance (m : ℕ) : (m : R[X]).ComputableRepr := .ofEq _ _ (Polynomial.C_eq_nat_cast m)

variable [q.ComputableRepr]

instance : (p + q).ComputableRepr where
  coeff n := coeff p n + coeff q n
  degreeBound := max (degreeBound p) (degreeBound q)
  coeff_eq _ := by simp_rw [Polynomial.coeff_add, coeff_eq]
  support_subset := Polynomial.support_add.trans <|
    (Finset.union_subset_union support_subset support_subset).trans fun _ _ ↦ by simp_all

def Nat.withBotAdd : ℕ → ℕ → ℕ
  | 0, _ => 0
  | _, 0 => 0
  | m + 1, n + 1 => m + n + 1

@[simps] def WithBot.natSuccIso : WithBot ℕ ≃o ℕ where
  toFun m := m.recBotCoe 0 (· + 1)
  invFun m := m.rec ⊥ (↑)
  left_inv m := by cases m <;> rfl
  right_inv m := by cases m <;> rfl
  map_rel_iff' {m n} := by
    cases m; · simp [recBotCoe]
    cases n; · simpa [recBotCoe] using WithBot.bot_lt_coe _
    simp [recBotCoe]

def WithBot.natSucc (n : WithBot ℕ) : ℕ := WithBot.natSuccIso n

lemma Polynomial.natSucc_degree (hp : p ≠ 0) : p.degree.natSucc = p.natDegree + 1 := by
  rw [natDegree, ← WithBot.coe_unbot _ (degree_eq_bot.not.mpr hp)]; rfl

lemma WithBot.natSucc_add (m n : WithBot ℕ) : m.natSucc.withBotAdd n.natSucc = (m + n).natSucc := by
  cases m <;> cases n <;> rfl

theorem Polynomial.support_subset_range : p.support ⊆ Finset.range p.degree.natSucc := fun n h ↦ by
  by_cases hp : degree p = ⊥
  · simp_all
  · rw [natSucc_degree p (degree_eq_bot.not.mp hp), Finset.mem_range_succ_iff]
    contrapose! h
    exact not_mem_support_iff.mpr (Polynomial.coeff_eq_zero_of_natDegree_lt h)

open scoped BigOperators

instance : (p * q).ComputableRepr where
  coeff n := ∑ i in Finset.range (min (degreeBound p) (n + 1)), coeff p i * coeff q (n - i)
  degreeBound := (degreeBound p).withBotAdd (degreeBound q)
  coeff_eq n := by
    simp_rw [← coeff_eq, Polynomial.coeff_mul, Finset.Nat.sum_antidiagonal_eq_sum_range_succ_mk]
    rw [← Finset.sum_range_add_sum_Ico _ (min_le_right _ _)]
    convert add_zero _
    refine Finset.sum_eq_zero fun n h ↦ ?_
    rw [Finset.mem_Ico, min_le_iff, ← not_lt, or_iff_left h.2.not_le, ← Finset.mem_range] at h
    rw [Polynomial.not_mem_support_iff.mp fun hn ↦ h.1 (support_subset hn), zero_mul]
  support_subset := (p * q).support_subset_range.trans <| Finset.range_mono <|
    (WithBot.natSuccIso.monotone <| p.degree_mul_le q).trans <| by
      rw [← WithBot.natSucc, ← WithBot.natSucc_add]
      sorry

noncomputable instance (priority := mid) : (n : ℕ) → (p ^ n).ComputableRepr :=
  Nat.rec (ofEq 1 _ (pow_zero _).symm) fun n _ ↦ ofEq (p * p ^ n) _ (pow_succ _ _).symm

instance (m : ℕ) : (X ^ m : R[X]).ComputableRepr where
  coeff n := if n = m then 1 else 0
  degreeBound := m + 1
  coeff_eq _ := by simp_rw [Polynomial.coeff_X_pow]
  support_subset := by
    rw [← Polynomial.monomial_one_right_eq_X_pow]
    exact (Polynomial.support_monomial' _ _).trans (by simp)

noncomputable instance {S} [Semiring S] (f : R →+* S) : (p.map f).ComputableRepr where
  coeff n := f (coeff p n)
  degreeBound := degreeBound p
  coeff_eq _ := by rw [Polynomial.coeff_map, coeff_eq]
  support_subset := (Polynomial.support_map_subset _ _).trans support_subset

section

variable {S} [Ring S] (p q : S[X]) [p.ComputableRepr] [q.ComputableRepr]

instance : (-p).ComputableRepr where
  coeff n := - coeff p n
  degreeBound := degreeBound p
  coeff_eq _ := by rw [Polynomial.coeff_neg, coeff_eq]
  support_subset := by rw [Polynomial.support_neg]; exact support_subset

-- `coeff p n - coeff q n` more efficient?
noncomputable instance : (p - q).ComputableRepr := .ofEq (p + -q) _ (sub_eq_add_neg _ _)

end

theorem Polynomial.ComputableRepr.eq_iff : p = q ↔
    ∀ i ∈ Finset.range (max (degreeBound p) (degreeBound q)), coeff p i = coeff q i :=
  ⟨fun h _ _ ↦ by simp_rw [← coeff_eq, h], fun h ↦ ext fun n ↦ by
    simp only [Finset.mem_range, lt_max_iff] at h
    by_cases hn : ?_
    · simp_rw [coeff_eq]; exact h n hn
    · simp_rw [not_or, not_lt] at hn; rw [coeff_eq_zero_of_le hn.1, coeff_eq_zero_of_le hn.2]⟩

instance [DecidableEq R] : Decidable (p = q) :=
  decidable_of_iff _ (Polynomial.ComputableRepr.eq_iff p q).symm

example : (X + X : ℤ[X]) = (2 * X : ℤ[X]) := by
  rw [← Nat.cast_eq_ofNat] -- have to do this
  decide

noncomputable def H : ℚ[X] := X ^ 3 - (C (6/8) * X) - C (1/8)
noncomputable def H' : ℤ[X] := C 8 * X ^ 3 - (C 6 * X) - C 1

example : DecidableEq ℚ := inferInstance -- okay
-- example : (6 / 8 : ℚ) * 8 = 6 := by decide -- doesn't work because operations on ℚ are irreducible
-- see https://leanprover.zulipchat.com/#narrow/stream/113489-new-members/topic/how.20to.20explicitly.20rfl.20for.20.E2.84.9A.3F/near/406961128

lemma H_H' : (C 8) * H = H'.map (algebraMap ℤ ℚ) := by rw [H, H']; apply of_decide_eq_true; rfl
-- `of_decide_eq_true <| by rw [H, H']; rfl` doesn't work. Maybe the rewritten Decidable instance is not good?

lemma H'_degree : H'.natDegree = 3 := by rw [H', natDegree_eq]; rfl
lemma H_degree : H.natDegree = 3 := by rw [H, natDegree_eq]; rfl
lemma H_monic : H.Monic := by rw [H, monic_iff]; rfl -- decide doesn't work (because it's in ℚ)
lemma H_ne_zero : H ≠ 0 := by rw [H]; apply of_decide_eq_true; rfl
lemma H'_coeff_zero : H'.coeff 0 = -1 := by rw [H', coeff_eq]; rfl
lemma H'_leading : H'.leadingCoeff = 8 := by rw [H', leadingCoeff_eq]; rfl

set_option profiler true
example :
    (1 + X + X ^ 2 - X ^ 5 - X ^ 6 - 2 * X ^ 7 - X ^ 8 - X ^ 9 + X ^ 12 + X ^ 13 + X ^ 14 +
        X ^ 15 + X ^ 16 + X ^ 17 - X ^ 20 - X ^ 22 - X ^ 24 - X ^ 26 - X ^ 28 + X ^ 31 + X ^ 32 +
        X ^ 33 + X ^ 34 + X ^ 35 + X ^ 36 - X ^ 39 - X ^ 40 - 2 * X ^ 41 - X ^ 42 - X ^ 43 +
        X ^ 46 + X ^ 47 + X ^ 48 : ℤ[X]).degree = 48 := by
  rw [← Nat.cast_eq_ofNat (R := ℤ[X]) (n := 2), degree_eq]; rfl
/- typeclass inference of Polynomial.ComputableRepr took 159ms
   elaboration took 314ms -/

example :
    (1 + X + X ^ 2 - X ^ 5 - X ^ 6 - 2 * X ^ 7 - X ^ 8 - X ^ 9 + X ^ 12 + X ^ 13 + X ^ 14 +
        X ^ 15 + X ^ 16 + X ^ 17 - X ^ 20 - X ^ 22 - X ^ 24 - X ^ 26 - X ^ 28 + X ^ 31 + X ^ 32 +
        X ^ 33 + X ^ 34 + X ^ 35 + X ^ 36 - X ^ 39 - X ^ 40 - 2 * X ^ 41 - X ^ 42 - X ^ 43 +
        X ^ 46 + X ^ 47 + X ^ 48 : ℚ[X]).degree = 48 := by
  rw [← Nat.cast_eq_ofNat (R := ℚ[X]) (n := 2), degree_eq]; rfl
/- typeclass inference of Polynomial.ComputableRepr took 172ms
   elaboration took 352ms -/

[Junyan Xu]

It turns out there's an irreducibility issue of operations in Rat which blocks decide. However rfl (the tactic, not the term) still works, and you can use docs#of_decide_eq_true to convert any goal to an equality goal and apply rfl, and all DecidableEq instances work this way. Now I can also compute natDegree, degree, leadingCoeff and Monic by reflection (see updated code above).

[Damiano Testa]

This is very cool! I'm not at a computer now, but I wonder if this is the way to go, instead of the tactic.

Can you try your approach on the tests for compute_degree?

I'm especially curious about performance on the slowest: the computations with the 105 cyclotomic polynomial. If I remember correctly, the tactic takes about 3 seconds to compute the degree.

[Junyan Xu]

I think my method can't handle the tests with [Nontrivial R], since it requires [DecidableEq R], so your tactics are certainly more powerful in that aspect!

I tested the cyclotomic 105 example after adding a dedicated instance for (X^n).ComputableRepr, and it also takes about 3 seconds for ℤ[X] and ℚ[X]. (Profiler is reporting a lower number, and I'm not sure how to get it to report actual kernel checking time.) I updated the code above, and there are two remaining sorries which I'll take care tomorrow :)

[Junyan Xu]

Okay, my previous implementation of degree/degreeAux is probably slightly buggy, but everything is now sorry free at https://gist.github.com/alreadydone/2dca4fde11fb2e9be7f8a10b59216b3f
Some proofs may be golfable with appropriate lemmas.
For now I'll focus on other stuff so I'm not adding support for other ways of constructing polynomials, but you're welcome to experiment with it! I wonder where else this technique can be applied to. (One place is Finsupp addition which also produces non-computable support which makes equality non-decidable.)
It would be nice if we could make the ofNat instance work.
To show the 105th cyclotomic polynomial is of degree 48 over an arbitrary nontrivial ring, one way would be proving a lemma that if coeff p (degreeBound p - 1) is nonzero then degree is equal to degreeBound p - 1 (in this case, degreeBound p already computes to 49). However this coefficient computes (defeq) to 0 + 0 + 0 + ... + 0 + 1, so you may need to apply some (add_zero _).trans in term mode ...

[Damiano Testa]

It seems that your approach leverages typeclass inference to compute degrees, while the tactic uses... well meta-programming.

So, another possible performance issue might be that the tactic discharges some goals using norm_num, while your approach would resort to... decide? That could lead to slower proofs, especially if the coefficients are big (so, possibly an artificial case, for the moment!).

[Damiano Testa]

Anyway, I think that this is a very good experiment!

[Junyan Xu]

Thanks! Yeah I think typeclass inference is a nice way to match syntax, and it's used to build up a computable representation of an expression involving operations for which instances are provided; I don't know how its performance compare to tactics. I think there are more people familiar with typeclasses than there are with meta-programming, so my approach might have a lower entry barrier :)

With Lean 4's GMP arithmetic, I think there's probably no advantage of norm_num over rfl for arithmetic in Nat, Int, and Rat. norm_num extensions (primality, Legendre symbols, etc.) are nice but I think they could also be implemented by as functions (like degree, leadingCoeff etc. in the Polynomial.ComputableRepr namespace in my code).

Non-computability of Finsupp has been a discussion topic a lot before, but most recently @Geoffrey Irving's question prompted me to identify main difficulties of polynomial computability: a degree bound is necessary to offer a finite range for the summation in polynomial multiplication and for equality testing, so it's built into my ComputableRepr together with a proof that it's actually a bound. I think a natural next step is to compute in AdjoinRoot using this approach, which provides a way to address his original question.

In previous discussions people seemed to think it's necessary to redefine the Finsupp/Polynomial types to gain computability, but my approach here shows that:

1. we could make an auxiliary type and "implement" computable operations there, and make no modification to the primary type at all.
2. the auxiliary type doesn't have to be "isomorphic" to the primary type.
3. there doesn't even need to be a function from the primary type to the auxiliary type; the mapping could depend on the syntactic form of the term of the primary type via typeclass inference, if we make the auxiliary type a class. (For example, X - X, 1 - 1 and 0 are mapped to non-equal ComputableRepr terms with degreeBound = 2, 1 and 0 respectively.)
   (I tried to put the @[csimp] attribute on degree_eq, leadingCoeff_eq etc. but got an error; probably Lean doesn't like dependence on the instance argument.)

@Yakov Pechersky's #9973 (Zulip) attempts to use List as the auxiliary type, but the toList there doesn't use typeclass inference to match the syntax; instead it uses natDegree which relies on a computable support to be computable, so it won't actually provide a computable list when applied to expressions outside of trivial ones. I'm also worried that a list representation may not be efficient for computing the degree of X^10000+1, for example.

I'd also like to hear opinions from @Mario Carneiro, @Kyle Miller, and @Eric Wieser. Feel free to move the discussion to another topic/stream as appropriate.

[Mario Carneiro]

With Lean 4's GMP arithmetic, I think there's probably no advantage of norm_num over rfl for arithmetic in Nat, Int, and Rat.

I think this is only true for a specific enumerated set of functions. For functions which are not built in norm_num still has the edge, certainly anything defined by recursion.

[Yakov Pechersky]

The values in these type class instances seem prone to non defeq diamonds since they contain data

[Junyan Xu]

These instances are just for generating computable representations. If Lean finds two different instances/representations for the same expression (potentially not prop-eq), they may differ in efficiency, but will still compute the same correct result when used to compute coefficients, degree, etc.

(The way to compute is to rewrite using one of the lemmas coeff_eq, degree_eq, natDegree_eq, leadingCoeff_eq, monic_iff under the Polynomial.ComputableRepr namespace, and then use rfl or decide (if not dealing with Rat). Lean may not find the instance that you intend when doing the rewrite, but any instance works.)

[Damiano Testa]

Yaël Dillies said:

Damiano Testa said:

Yes, ext entering Rat is always an issue, in my experience.

That can be fixed by removing the ext attribute from the relevant lemma. We did the same thing for ℂ recently.

The ext attribute on Rat has just been removed from Std. It will take a little bit of time to trickle down to Mathlib. but, according to a local build that I made, there should be only one instance in Mathlib where the ext attribute on Rat is used.

[Kim Morrison]

@Damiano Testa (or anyone else here) #10954 is the long overdue bump for Mathlib, that is now broken because of the missing Rat.ext. Could you please fix this asap? This is on the critical path to being able to release v4.7.0rc1 this week.

[Kim Morrison]

(Ideally PRs to Std or Lean that break Mathlib should include a clear reference in the PR description to the location of the patch branch, so that others can help deploy the patch once it is needed.)

[Riccardo Brasca]

You mean making #10954 compile?

[Riccardo Brasca]

I can do it now

[Kim Morrison]

Great, thanks!

[Kim Morrison]

(This could have been done as soon as the Rat.ext removal landed in Std. I've been behind on bumping, but everyone is allowed to do this. :-)

[Kim Morrison]

(@Riccardo Brasca, I just pushed a change to #10954, so we are not updating aesop at the same time. If you've already pull, pull again please. :-)

[Kim Morrison]

Ugh, there are multiple problems, because we haven't bumped for a while.

[Kim Morrison]

I will take care of the Int and List stuff.

[Eric Wieser]

Can we remove Rat.ext in mathlib before bumping to reduce the amount of debt, or is that not how attribute removal works?

[Kim Morrison]

I'm not sure what you mean, sorry.

[Eric Wieser]

Rat.ext was de-@[ext]d in a version of Std that we haven't bumped to yet, right? Can we locally de-attribute it in Mathlib to deal with the fallout without having to fix everything else in the Std bump at the same time?

[Kim Morrison]

Oh, I see. I think we need to just do it. :-) I need the entire bump anyway asap, so I can continue work on the nightly-testing branch.

[Kim Morrison]

#10954 is looking okay now other than Mathlib/Data/NNRat/Defs.lean

[Eric Wieser]

Mainly the suggestion was to allow the fixes to be parallelized with someone else:)

[Kim Morrison]

I think we're parallelizing by pushing to the same branch.

[Riccardo Brasca]

I am fixing /Data/NNRat/Defs.lean, without pulling

[Kim Morrison]

Great, thanks.

[Kim Morrison]

I've fixed everything that is not downstream of that.

[Riccardo Brasca]

If the goal was simply to fix it, it's done (I didn't do anything fancy, I just applied manually the lemma that was called by ext).

[Kim Morrison]

Great, thanks!

[Riccardo Brasca]

I am also fixing the later error in sum of four squares

[Riccardo Brasca]

it builds locally

[Kim Morrison]

Thanks @Riccardo Brasca!

[Riccardo Brasca]

There is a linter error about Std.BitVec.toFin_natCast in Data/BitVec/Lemmas

[Riccardo Brasca]

I am not sure how to fix it, currently simp transforms

toFin (n : BitVec w) = n

into

Fin.ofNat' n _ = n

and I am not sure if this is a good idea

[Damiano Testa]

Sorry about getting late here: I had built Mathlib locally with the not-yet removal of the ext attribute on Std and had found (I think) a single issue. However, I did not really know where I was supposed to make the fix, even after Scott explained some of the details of how to deal with changes to std. I'll try to clean up after myself in the future.