Stream: new members

Topic: Can't define Chebyshev polynomials

Abhimanyu Pallavi Sudhir (Jan 12 2019 at 17:34):

Why doesn't this work?

def chebyshev : ℕ → polynomial ℝ
| 0 := polynomial.C 1
| 1 := polynomial.X
| (n + 2) := 2 * polynomial.X * chebyshev (n + 1) - chebyshev n


Lean tells me the definition relies on classical.choice. Indeed putting noncomputable at the front fixes things, but why is it noncomputable?

Kevin Buzzard (Jan 12 2019 at 17:52):

The definition of polynomial X is "functions from nat to X which vanish outside a finite set"; "vanish" means "equals zero", and equality is undecidable in the reals. This might be something to do with it. Try polynomial int and see if it fixes things!

Kevin Buzzard (Jan 12 2019 at 17:54):

NB import analysis.polynomial to make it work, but indeed int fixes things, and is probably the "correct" thing to do.

Kevin Buzzard (Jan 12 2019 at 18:00):

import analysis.polynomial

def chebyshev : ℕ → polynomial ℤ
| 0 := polynomial.C 1
| 1 := polynomial.X
| (n + 2) := 2 * polynomial.X * chebyshev (n + 1) - chebyshev n

#eval chebyshev 17 -- just about works on my desktop
/-
C (65536) * X ^ 17 + C (-278528) * X ^ 15 + C (487424) * X ^ 13 + C (-452608) * X ^ 11 + C (239360) * X ^ 9 + C (-71808) * X ^ 7 + C (11424) * X ^ 5 + C (-816) * X ^ 3 + C (17) * X
-/

#eval chebyshev 19 -- (deterministic) timeout


Mario Carneiro (Jan 12 2019 at 18:06):

you can probably get farther with a list based representation of the polynomials

Kevin Buzzard (Jan 12 2019 at 18:07):

I guess the reason lists aren't used in general is that the leading term might be 0, which causes problems in general; however it will not cause problems here.

Mario Carneiro (Jan 12 2019 at 18:07):

well, it would make it harder to print if the leading term was 0

Mario Carneiro (Jan 12 2019 at 18:07):

and get the degree and other things

Kevin Buzzard (Jan 12 2019 at 18:07):

In general I guess one could define some more computationally efficient but slightly broken "non-zero polynomials" with addition only defined if you can prove that the degrees are unequal, and multiplication OK over an integral semi-domain :-)

Kevin Buzzard (Jan 12 2019 at 18:08):

But even lists are inefficient, right, because they're linked lists in reality?

Mario Carneiro (Jan 12 2019 at 18:09):

If you want to forgo the equality checks, you can define a polynomial as a quotient of representations with some zeros at the end

Mario Carneiro (Jan 12 2019 at 18:09):

then addition and multiplication are easy and degree needs equality

Mario Carneiro (Jan 12 2019 at 18:09):

yes lists are linked lists which aren't so great

Mario Carneiro (Jan 12 2019 at 18:10):

you could use buffer for really good VM performance

Mario Carneiro (Jan 12 2019 at 18:10):

but they are better than functions which is what is currently used

Mario Carneiro (Jan 12 2019 at 18:10):

I think the polynomials become big if else chains

Mario Carneiro (Jan 12 2019 at 18:12):

in fact, they probably aren't even reduced, you get these big thunks for calculating the coefficients

Mario Carneiro (Jan 12 2019 at 18:12):

the advantage of lists is they are a strict data structure, you calculate all the values before recursing

Mario Carneiro (Jan 12 2019 at 18:13):

also -- should have mentioned this before -- chebyshev there has an exponential time implementation

Mario Carneiro (Jan 12 2019 at 18:13):

just like the naive fib algorithm

Daniel Keys (Feb 13 2020 at 18:02):

It seems some things may have changed since this discussion took place. Trying to get started working with polynomials; analysis.polynomial is not a valid import any longer. In the following I try to reproduce Kevin's piece of code listed above in this same topic, but it doesn't work on my computer because of the VM not having code for classical.choice. Am I missing something? Also, the definition of p1 is non-computable; do we have a way to compute (like #eval) a value from it (or any other polynomial on the reals, or even nats) at a given point?

import data.polynomial
import data.real.basic

-- this is noncomputable
def p1 : polynomial ℝ := polynomial.C (2:ℝ)
-- these check but I can't produce a value
#check polynomial.eval ( 5 : ℝ ) ( polynomial.C (2:ℝ) )
#check polynomial.C 1

def chebyshev : ℕ → polynomial ℤ -- fails,  VM: no classical.choice
| 0 := polynomial.C 1
| 1 := polynomial.X
| (n + 2) := 2 * polynomial.X * chebyshev (n + 1) - chebyshev n

#eval chebyshev 3 --fails, VM: no classical.choice


Johan Commelin (Feb 13 2020 at 18:09):

choice will certainly be uncomputable. That's a fact of life that we'll need to deal with.

Johan Commelin (Feb 13 2020 at 18:09):

So you will want

noncomputable theory


at the top of your file, just after the imports

Johan Commelin (Feb 13 2020 at 18:11):

Next, #eval won't work, so you can't simply compute the value. But you could write some simp-lemmas that can do some work for you.

Johan Commelin (Feb 13 2020 at 18:11):

Which is almost as good.

Johan Commelin (Feb 13 2020 at 18:13):

It is unfortunate that it is currently not easy to run the simplifier against an expression like chebyshev 3. You need to write example : chebyshev 3 = the_answer := by simp [chebyshev], which defeats the purpose because you will have to precompute the answer. And then Lean will tell you that it agrees :expressionless:

Rob Lewis (Feb 13 2020 at 18:13):

Maybe we should define a command #simp that behaves like #eval except using the simplifier.

Daniel Keys (Feb 13 2020 at 18:16):

Just wondering how Kevin could #eval chebyshev 17, as he shows above (that is, the older messages on the same topic).

Rob Lewis (Feb 13 2020 at 18:18):

#1391 removed this. It worked on toy examples, but isn't very helpful in practice and causes problems elsewhere.

Kevin Buzzard (Feb 13 2020 at 19:58):

Chris is sitting next to me talking about perhaps having some attribute for computing and then you just unfold a lot until you can use norm_num or something...

Rob Lewis (Feb 13 2020 at 20:27):

I think that's tricky. Do you unfold has_add.add? Sometimes yes, because you want to unfold addition on whatever type you're working on. Sometimes no, because you want norm_num to deal with it.

Last updated: May 09 2021 at 19:11 UTC