Playing with polynomials

Stream: new members

[Lucas Allen]

Hello, in an effort to play with polynomials in Lean I came up with the following,

import data.polynomial

open polynomial

noncomputable def P : polynomial ℤ := 3 * X^2 + 1
noncomputable def Q : polynomial ℤ := 4 * X^2 + 7 * X + 6

#check P + Q
#check P - Q
#check P * Q
#check P /ₘ Q

#eval P + Q
#eval P - Q
#eval P * Q
#eval P /ₘ Q

The #eval lines give the error "code generation failed, VM does not have code for 'classical.choice'", and the definitions of P and Q demand the noncomputable. Is there some way around the noncomputability of polynomials, so that #eval P + Q gives 7 * X^2 + 7 * X + 7?

[Kenny Lau]

"If you want to calculate, use Python"

[Kenny Lau]

(or use an older version of Lean/mathlib before we decided to make everything noncomputable)

[Kenny Lau]

probably some version in 2018

[Lucas Allen]

Ok.

[Reid Barton]

Or use a future version of Lean/mathlib with the norm_poly tactic

[Kevin Buzzard]

@Lucas Allen can you please leave 2020, you've come to the wrong year. Please recalibrate a couple of years in either direction.

[Mario Carneiro]

I need a telephone booth emoji

[Kevin Buzzard]

Polynomials used to be computable but this caused trouble with reasoning. Now they're noncomputable but apparently this causes problems with computation. Is Reid right that this can be fixed with a tactic?

[Mario Carneiro]

yes

[Mario Carneiro]

just like norm_num but for polynomials

[Reid Barton]

The trick is the kind of "computation" you want is not the same as the noncomputable kind. For example, norm_num works on the reals too

[Kenny Lau]

that there are two meanings of "computation": the Lean meaning, which is that you write a program that does it; the mathematical meaning, which is that you reduce something to a normal form.

[Kevin Buzzard]

There's some trickery involved though, isn't there? Because norm_num doesn't _really_ work with reals, right? It can only make progress if all the reals are rational.

[Reid Barton]

It really works with exprs that represent reals

[Kevin Buzzard]

Aah!

[Kevin Buzzard]

So why don't we teach it that sqrt(x)*sqrt(x)=x :-)

[Reid Barton]

Like this?

import data.real.basic
import tactic.norm_num

example : real.sqrt 5 * real.sqrt 5 = 5 := by norm_num

[Kevin Buzzard]

Thanks!

[Kevin Buzzard]

(I never knew it did that -- I had the wrong mental model)

[Reid Barton]

I think the main obstruction to turning norm_num into something more like a computer algebra system is that it's not obvious (at least to me) what kind of organizational framework is appropriate

[Reid Barton]

along with perhaps a lack of specific applications at the moment

[Kevin Buzzard]

I want to prove that $x^2-4$ is not irreduible without doing any work.

[Reid Barton]

but for any computation you can do in a CAS, provided you can give a proof that the computation is correct (when it succeeds), there's no reason that something like norm_num couldn't be extended to do it

[Kevin Buzzard]

Actually I don't really want to prove this at all, I just want to prove stuff about DVRs etc. I don't have a specific application.

[Reid Barton]

This example is pretty trivial because ring can already do the correctness proof

[Kevin Buzzard]

I guess sometimes I want to prove that y^2=x^3-4x has rank 0 (or 1 or whatever).

[Reid Barton]

so you can "just" do the factorization in an untrusted program

[Kevin Buzzard]

and Cremona's software is not spitting out a formally verifiable thing here

[Reid Barton]

Right, I think you skipped over a large range of difficulties there :upside_down:

[Kevin Buzzard]

Nobody can prove the algorithm terminates, for starters

[Reid Barton]

Well that's not actually a problem

[Kevin Buzzard]

Try telling that to the equation compiler

[Kevin Buzzard]

Of course I see what you mean, tactics are meta

[Reid Barton]

For a procedure like this, there are two relevant outcomes:

• the procedure succeeds in producing an answer along with a correct proof
• all other behavior: nontermination, answer but no certificate, certificate but it doesn't type check somehow

[Reid Barton]

If you can build a procedure which provably always has the first behavior then great, but in practice it's fine to have tactics that can fail, as long as they work when you need them to

[Reid Barton]

since our goal is to compute the ranks of elliptic curves (or whatever), not prove an internal theorem "there exists an algorithm which computes the rank of every elliptic curve"

[Kevin Buzzard]

which is good because that's an open problem

[Mario Carneiro]

FYI, this was a surprise to me as well. norm_num doesn't know anything about square roots, that's all simp's work

[Reid Barton]

aha!

[Bryan Gin-ge Chen]

Returning to the initial question, right now #simp [P,Q] P + Q just returns 3 * X ^ 2 + 1 + (4 * X ^ 2 + 7 * X + 6), but maybe it's possible to add some simp lemmas (or make a simp set?) to make it more "useful".

[Lucas Allen]

Kevin Buzzard said:

Lucas Allen can you please leave 2020, you've come to the wrong year. Please recalibrate a couple of years in either direction.

I think practicing the art of impeccable timing is starting to pay off.

[Lucas Allen]

Bryan Gin-ge Chen said:

Returning to the initial question, right now #simp [P,Q] P + Q just returns 3 * X ^ 2 + 1 + (4 * X ^ 2 + 7 * X + 6), but maybe it's possible to add some simp lemmas (or make a simp set?) to make it more "useful".

I had a go at adding a simp lemma to see if it worked,

@[simp] lemma X_add : (X : polynomial ℤ) + X = 2 * X := sorry

#simp [X] X + X

It gives "monomial 1 1 + monomial 1 1", but I want "2 * X". I also tried #simp [(X : polynomial ℤ)] X + X but this gives "simplify tactic failed to simplify." Is there some way to get this to work?

[Johan Commelin]

Ooh, you don't want to put X in the simpset, rather put X_add between the []

[Lucas Allen]

#simp [X_add] X + X gives "simplify tactic failed to simplify". Do I need to say something about it being a ℤ polynomial?

[Johan Commelin]

Aah, I guess you do

[Johan Commelin]

@[simp] lemma X_add : (X : polynomial ℤ) + X = 2 * X := sorry

#simp [X_add] (X : polynomial ℤ) + X

does this work?

[Lucas Allen]

Yes!

[Lucas Allen]

Why is this an invalid simp lemma?

lemma X_add (m : ℕ) (a b : ℤ) : (a * X : polynomial ℤ) + (b * X) = (a + b) * X :=
begin
  exact (add_mul ↑a ↑b X).symm,
end

[Lucas Allen]

Adding @[simp] gives a red squiggle.

[Johan Commelin]

@Lucas Allen because the statement doesn't use m

[Lucas Allen]

Oh no.

[Lucas Allen]

sorry.

[Lucas Allen]

I get a deterministic timeout with the following

@[simp] lemma X_add (a b : ℤ) : (a * X : polynomial ℤ) + (b * X) = (a + b) * X :=
begin
  exact (add_mul ↑a ↑b X).symm,
end

lemma test : (3 * X : polynomial ℤ) + (5 * X) = 8 * X :=
begin
  apply X_add,
end

I don't know why, and I hope it's not obvious again...

[Johan Commelin]

What happens if you rw instead of apply?

[Johan Commelin]

I would have expected the apply to work though.

[Lucas Allen]

"rewrite tactic failed, did not find instance of the pattern in the target expression
↑?m_1 * X + ↑?m_2 * X
state:
⊢ 3 * X + 5 * X = 8 * X"

no deterministic timeout though.

[Lucas Allen]

I'm wondering if it has something to do with the fact that a * X has type polynomial ℤ, and it's not (a : ℤ) * X.

[Lucas Allen]

But (a : ℤ) * X doesn't work.

[Lucas Allen]

lemma X_add (a b : ℤ) : a * (X : polynomial ℤ) + b * X = (a + b) * X :=gives red squiggles under the *'s. It likes lemma test : 3 * (X : polynomial ℤ) + 5 * (X) = 8 * X := however.

[Johan Commelin]

Try (3 : \Z) instead

[Johan Commelin]

but rather, I think you should use C to create constant polynomials

[Johan Commelin]

Arguably, we should consider making a \bu X the standard way of multiplying by scalars, instead of C a * X.

[Lucas Allen]

lemma xxx : (3 : ℤ) * (X : polynomial ℤ) + (5 : ℤ) * (X) = (8 : ℤ) * X := gives red squiggles under the *'s type mismatch error.

[Johan Commelin]

((3 : \Z) * X : ...) ?

[Lucas Allen]

I think it was (3 : \Z) * (X : ...)

[Lucas Allen]

It likes the \bu. Everything works now.

@[simp] lemma X_add (a b : ℤ) : a • (X : polynomial ℤ) + b • X = (a + b) • X :=
begin
  exact (add_smul a b X).symm
end

lemma xxx : (3 : ℤ) • (X : polynomial ℤ) + (5 : ℤ) • (X) = (8 : ℤ) • X :=
begin
  apply X_add,
end

#simp [X_add] (3 : ℤ) • (X : polynomial ℤ) + (5 : ℤ) • (X) --(3 + 5) • X

[Stanislas Polu]

This is the best I could come up with but I have a feeling it's not idiomatic (or correct?). Is there a way to eval directly into polynomial R ?

theorem aopsbook_v2_c6_em2
  (f : polynomial ℝ)
  (a : ℝ) :
  polynomial.mod f (polynomial.X - polynomial.C a) = polynomial.C (polynomial.eval a f) :=
begin
  sorry
end

[Johan Commelin]

There is eval\_2, I think. It takes a ring hom as argument

[Johan Commelin]

Maybe

open polynomial

... = aeval f (C a)

[Stanislas Polu]

Well from the docs#polynomial.eval:

polynomial.eval = polynomial.eval₂ (ring_hom.id R)

so it should land in R as well?

[Johan Commelin]

No, if you pass it something different from ring_hom.id R... then you can land in another ring

[Stanislas Polu]

Hmmm aeval f (C a) typechecks.

[Stanislas Polu]

theorem aopsbook_v2_c6_em2
  (f : polynomial ℝ)
  (a : ℝ) :
  polynomial.mod f (polynomial.X - polynomial.C a) = polynomial.aeval f (polynomial.C a) :=
begin
  sorry
end

[Anne Baanen]

Doesn't aeval f (C a) mean "evaluate the constant polynomial C a at the point f? docs#polynomial.aeval

[Johan Commelin]

Ooh, so maybe you need to swap the arguments :see_no_evil:

[Anne Baanen]

Yes, we want instead aeval (C a) f, which is the same as C (aeval a f) = C (eval a f). (See e.g. my PR #8038)

[Stanislas Polu]

Fantastic!

[Stanislas Polu]

Follow up question. I'd like to state the rational root theorem using a non generalized form (polynomial with integer coefficients and consider rational roots).

I therefore need to state that the polynomial has integer coefficient and that a fraction in ℚ is a root of the polynomial but that does not typechecks because I'm no authorized to evaluate the polynomial in ℚ if it's defined as polynomial ℤ. Example:

theorem aopsbook_v2_c6_em6
  (P : polynomial ℤ)
  (p q : ℤ)
  (h₀ : p.gcd q = 1)
  (h₁ : is_root P (p:ℚ/q)) :
  p ∣ P.coeff 0 ∧ q ∣ P.leading_coeff :=
begin
  sorry
end

What exactly does polynomial R mean? Does it mean that the coeffs are in R? Can one define a polynomial whose coefficient are in R1 but codomain/domain R2 ?

[Adam Topaz]

I think docs#polynomial.aeval should work?

[Stanislas Polu]

Rephrasing, I'm surprised that eval : R → polynomial R → R or is_root : polynomial R → R → Prop and not is_root: polynomial R → S → Prop ?

[Stanislas Polu]

@Adam Topaz what would be the syntax to use docs#polynomial.aeval here?

[Stanislas Polu]

(and I would use that to specify that p/q is a root) ?

[Adam Topaz]

One sec I'll try to sketch some code

[Riccardo Brasca]

You can use docs#polynomial.map to interpret a polynomial with integer coefficients as a polynomial with rational coefficients

[Adam Topaz]

Sure, you can do that too, but I think you can just use aeval in this case.

[Stanislas Polu]

Interesting, and out of curiousity how do I define (f : R →+* S) for ℤ and ℚ ?

[Stanislas Polu]

I presume it's fraction_map ℤ ℚ assuming that type checks ? :grimacing:

[Anne Baanen]

Almost, algebra_map ℤ ℚ

[Riccardo Brasca]

But I agree with Adam that aeval Is better

[Anne Baanen]

In mathlib we tend to define an algebra R S instance if there's a canonical map algebra_map R S : R →+* S.

[Adam Topaz]

Something like this is what I had in mind:

import ring_theory.polynomial.basic

#check @polynomial.aeval

example : algebra ℤ ℚ := by apply_instance

variables (f : polynomial ℤ) (p q : ℤ) (h : p.gcd q = 1)

example (hh : polynomial.aeval ((p : ℚ) / q) f  = 0) : true := sorry

[Adam Topaz]

(sorry, the #check and first example were just a sanity check for myself ;) )

[Riccardo Brasca]

We also have docs#int.cast_ring_hom

[Anne Baanen]

docs#fraction_map is more like ring_hom, in that fraction_map ℤ ℚ is the type of maps ℤ → ℚ that display ℚ as the fraction field of ℤ

[Stanislas Polu]

Thanks @Adam Topaz!

[Stanislas Polu]

And thanks @Riccardo Brasca and @Anne Baanen things are getting much clearer to me \o/

[Adam Topaz]

Unfortunately it looks like we can't use dot notation for polynomial.aeval because it gives the evaluation as a ring hom.

[Adam Topaz]

I.e. f.aeval ... doesn't work

[Kevin Buzzard]

Try it in Lean 4!

[Stanislas Polu]

So just to make sure I understand, aeval infers the mapping based on the types of the arguments and the right mapping is defined somewhere for it to pick it up automatically?

[Eric Wieser]

Yes, there's a long list of "canonical" mappings at docs#algebra

[Eric Wieser]

(under > Instances)

[Adam Topaz]

@Stanislas Polu the a in aeval stands for algebra. So it infers the right mapping based on an algebra instance. Fortunately for you, any ring has a unique structure of a Z-algebra and I think this is an actual instance in mathlib. So presumably you should be able to use aeval to evaluate an integral polynomial at an element in any ring.

[Stanislas Polu]

Thanks! Makes sense :+1:

[Stanislas Polu]

Hi! Me again :wave: What would be an elegant way to define the k-th symmetric sum of the roots of a polynomial to prove that it's (-1)^k*a_{n-k}/a_n ?

[Stanislas Polu]

I've been looking for that theorem to no avail, and have no good idea on how to define such a summation?

[Johan Commelin]

@Stanislas Polu Note that there is src/ring_theory/polynomial/symmetric.lean, but I don't think it contains exactly what you want.

[Stanislas Polu]

Thanks! will look into it :+1:

[Stanislas Polu]

(Seems quite unrelated indeed)

[Johan Commelin]

Well, a more fleshed out version of that file would hopefully make it very easy to talk about the k-th symmetric sum, etc...

[Stanislas Polu]

Relatedly how can I describe a sum on a multiset (to describe a sum that depends on the roots of a polynomial)?

[Anne Baanen]

If you have s : multiset α and you have a term f x for each x ∈ s that you want to add up, then you can use (s.map f).sum (see docs#multiset.map, docs#multiset.sum).

[Stanislas Polu]

@Anne Baanen Thanks!

[Stanislas Polu]

Fantastic!

import data.real.basic
import data.polynomial
open polynomial

theorem aopsbook_v2_c6_ex7
  (P : polynomial ℝ)
  (f : ℝ -> ℝ)
  (h₀ : P = (C 3)*X^3 - (C 14)*X^2 + X + C 62)
  (h₁ : ∀ r, f(r) = 1/(r+3)) :
  (P.roots.map f).sum = 83/74 :=
begin
  sorry
end

[Kevin Buzzard]

This does not look true to me. The cubic only has one real root and P.roots looks to me to be the multiset of real roots. Have I made a mistake?