Can't define Chebyshev polynomials

Stream: new members

[Abhimanyu Pallavi Sudhir]

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]

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]

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

[Kevin Buzzard]

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]

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

[Kevin Buzzard]

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]

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

[Mario Carneiro]

and get the degree and other things

[Kevin Buzzard]

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]

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

[Mario Carneiro]

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]

then addition and multiplication are easy and degree needs equality

[Mario Carneiro]

yes lists are linked lists which aren't so great

[Mario Carneiro]

you could use buffer for really good VM performance

[Mario Carneiro]

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

[Mario Carneiro]

I think the polynomials become big if else chains

[Mario Carneiro]

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

[Mario Carneiro]

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

[Mario Carneiro]

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

[Mario Carneiro]

just like the naive fib algorithm

[Daniel Keys]

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]

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

[Johan Commelin]

So you will want

noncomputable theory

at the top of your file, just after the imports

[Johan Commelin]

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]

Which is almost as good.

[Johan Commelin]

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]

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

[Daniel Keys]

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

[Rob Lewis]

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

[Kevin Buzzard]

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]

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.