Symbolic Chebyshev polynomials

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

Apply:

• attachment: trac_9706_chebyshev_collection.patch
• attachment: 9706_review.patch
• attachment: trac_9706_new_base_classes.patch
• attachment: trac_9706-review-ts.patch

Depends on #864 Depends on #9640 Depends on #10018 Depends on #11868 Depends on #15422

CC: @fredrik-johansson @sagetrac-fstan @kcrisman

Component: symbolics

Keywords: orthogonal polynomials, symbolics

Author: Stefan Reiterer

Branch/Commit: u/vbraun/chebyshev @ bb227b6

Reviewer: Burcin Erocal, Travis Scrimshaw, Stefan Reiterer, Jeroen Demeyer

Issue created by migration from https://trac.sagemath.org/ticket/9706

[jdemeyer]

comment:42

Sorry to spoil the party, but this is a regression:

sage: K.<a> = NumberField(x^3-x-1)
sage: chebyshev_T(10^3,a)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-17-aa97c56dd147> in <module>()
----> 1 chebyshev_T(Integer(10)**Integer(3),a)

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.BuiltinFunction.__call__ (sage/symbolic/function.cpp:8126)()

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:5279)()

TypeError: cannot coerce arguments: no canonical coercion from Number Field in a with defining polynomial x^3 - x - 1 to Symbolic Ring

(and yes: Chebyshev polynomials are important in number theory)

[jdemeyer]

comment:44

I would also like to point out that PARI is faster at evaluating Chebyshev polynomials:

sage: timeit('''chebyshev_T(10^5,2)''')
5 loops, best of 3: 270 ms per loop
sage: timeit('''pari('polchebyshev (10^5,1,2)')''')
625 loops, best of 3: 447 µs per loop

sage: timeit('''chebyshev_T(10^5, Mod(2,1009))''')
5 loops, best of 3: 208 ms per loop
sage: timeit('''pari('polchebyshev(10^5, 1, Mod(2,1009))')''')
625 loops, best of 3: 11.5 µs per loop

We should definately use PARI to evaluate Chebyshev polynomials for certain types of input.

[jdemeyer]

comment:45

Another regression:

sage: parent(chebyshev_T(10^2, RIF(2)))
Real Field with 53 bits of precision

[jdemeyer]

comment:46

The "Clenshaw method" uses a very naive method of evaluating the recursion which needs O(n) steps, while there is a much faster method (which compute T_2n and U_2n in function of T_n and U_n) which only needs O(log(n)) steps.

Even this is totally feasible:

sage: timeit('''pari('polchebyshev(10^10, 1, Mod(2,1009))')''')
625 loops, best of 3: 16.3 µs per loop

[jdemeyer]

comment:47

This is also bad:

sage: R.<x> = QQ[]
sage: parent(chebyshev_T(5, x))
Symbolic Ring

[jdemeyer]

comment:48

Another suggestion: why not make the base class more general and call it SymbolicPolynomial? I think that's a more natural class of functions, which could (in the future) also include cyclotomic polynomials for example.

[jdemeyer]

comment:49

I propose the logic for evaluating chebyshev_T(n, x) should be:

1. if x is symbolic, then use the method of the current patch.
2. if x is not symbolic, try evaluation using PARI.
3. if conversion to PARI fails (for example for RIF), use an efficient O(log(n)) recursion algorithm.

[fredrik-johansson]

comment:50

Jeroen: do you know a reference for the recursion pari uses?

[jdemeyer]

comment:51

Replying to @fredrik-johansson:

Jeroen: do you know a reference for the recursion pari uses?

No, but it's pretty straight-forward (think doubling formulas for cos and sin).

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:52

Thank you all for the input! I think it is still a good idea that I only implement chebychev polys for now, since there are a lot of improvements out there.

@jdemeyer

@Bugs I will look into this. And yes I'm aware that Chebyshev Polynomials are important to number theory since there are quite interesting generalizations on general fields.

@Clenshaw: PARI is a good hint, I will look into this. And I already think I know how to generalize the clenshaw method, to get O(log N). The reason for this quite naive choice, was that the method can be applied to all ortho polys. Nevertheless I will adapt it on Chebychev Polys since there are more possibilities since we can use trigonometric formulas. I think the benefit of implementing it directly in sage is that there is less trouble if one wants to use more general data types, since there are no type casts. I will try to find an optimal way for this. (Maybe an additional switch)

@SymbolicPolynomial: I don't think this is a good idea because ortho polys are quite special even among the polynomials. But If you really would like to have a SymbolicPolynomial class I would propose to introduce the SymbolicPolynomial class, and derive the OrthogonalPolynomials from that class. I make the following suggestion: I will finish the OrthogonalPolynomials with the current design. And then open an new ticket where we discuss the design of a general polynomial parent class. Fortunately, such design changes are very easy to implement in Python, and I don't see any big problem in introducing an intermediate class. But if you want to introduce such a class there are some major decisions to make: -Where do we put this class? (such a general class should not belong to orthogonoal_polys.py ) -What should all SymbolicPolynomials have in common? -What should they have concerning other general polynomials?

But it would be really good to add the ortho polys, and I really want to finish this task.

[jdemeyer]

comment:53

SymbolicPolynomial: I don't think this is a good idea because ortho polys are quite special even among the polynomials.

What's special about orthogonal polynomials from a computer algebra point of view? I can tell you that "symbolic polynomials" are special because you generally want to be able to evaluate them for any ring element (as opposed to other symbolic functions, which often only make sense in real or complex fields).

I make the following suggestion: I will finish the OrthogonalPolynomials with the current design.

Sure...

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:54

Replying to @jdemeyer:

SymbolicPolynomial: I don't think this is a good idea because ortho polys are quite special even among the polynomials.

What's special about orthogonal polynomials from a computer algebra point of view? I can tell you that "symbolic polynomials" are special because you generally want to be able to evaluate them for any ring element (as opposed to other symbolic functions, which often only make sense in real or complex fields).

It makes a difference in an OO-Design point of view, because you can apply a whole bunch of evaluation techniques and tricks due to the three term recursion (e.g. clenshaw method or the eval_numpy method are not needed for symbolic polynomials, since there are no methods for that). Thats the reason why I suggested to derive them from a base class on top to avoid redundant methods. That would be a clean solution. And I'm already thinking on some features, which I want to implement in future version, which are very specific to ortho polys. And if you want to introduce such a general class, you should not put it into orthogonal_polys.py, but somewhere else, were it fits better into the class hierarchy.

Please don't get me wrong, I'm not stating, that I think it is a bad idea to introduce such an abstract base class, but If you want a clean OO Design, it is not done by a simple renaming.

[jdemeyer]

comment:55

Evaluating chebyshev_T(n,x) can be done as

(Matrix(2,2,[x,x^2-1,1,x])^n)[0,0]

While chebyshev_U(n-1,x) equals

(Matrix(2,2,[x,x^2-1,1,x])^n)[1,0]

These can be evaluated with O(log(n)) operations.

[jdemeyer]

comment:56

Replying to @sagetrac-maldun:

And I'm already thinking on some features, which I want to implement in future version, which are very specific to ortho polys.

That alone is a very good to keep the OrthogonalPolynomial class.

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:57

Replying to @jdemeyer:

Evaluating chebyshev_T(n,x) can be done as

(Matrix(2,2,[x,x^2-1,1,x])^n)[0,0]

While chebyshev_U(n-1,x) equals

(Matrix(2,2,[x,x^2-1,1,x])^n)[1,0]

These can be evaluated with O(log(n)) operations.

Thanks for the hint. I have also an own idea to implement this. My implementation should be optimal with respect to the flop count, but yours could be faster since matrix multiplication and powers are well optimized. I will compare both methods and use the faster one.

For reference: The method which I mean is based on the generalized recursion formula (originating from the cosine addition theorem): T_{n+m} + T{n-m} = 2 T_n T_m

For T_N one can now use the binary representation of N to recursively build T_N in O(log N) time

[jdemeyer]

comment:58

This should be a ValueError:

sage: chebyshev_T(1/2,0)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-7-830f13ad2f0d> in <module>()
----> 1 chebyshev_T(Integer(1)/Integer(2),Integer(0))

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.BuiltinFunction.__call__ (sage/symbolic/function.cpp:8126)()                                                                                                                                                               

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:5531)()                                                                                                                                                                      

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _eval_(self, *args)
    489                 if not is_Expression(args[-1]):
    490                     try:
--> 491                         return self._evalf_(*args)
    492                     except AttributeError:
    493                         pass

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _evalf_(self, *args, **kwds)
    606             precision = step_parent.prec()
    607         except AttributeError:
--> 608             precision = RR.prec()
    609 
    610         from sage.libs.mpmath.all import call as mpcall

NameError: global name 'RR' is not defined

[jdemeyer]

comment:59

This should be ArithmeticError (I guess), since deriving w.r.t. the index simply isn't defined:

sage: var('n,x')
(n, x)
sage: chebyshev_T(n,x).diff(n)
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-14-a23f5209eb49> in <module>()
----> 1 chebyshev_T(n,x).diff(n)

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.derivative (sage/symbolic/expression.cpp:16561)()

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/misc/derivative.so in sage.misc.derivative.multi_derivative (sage/misc/derivative.c:2715)()

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._derivative (sage/symbolic/expression.cpp:16951)()

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _derivative_(self, *args, **kwds)
    714         diff_param = kwds['diff_param']
    715         if diff_param == 0:
--> 716             raise NotImplementedError("derivative w.r.t. to the index is not supported yet")
    717 
    718         return args[0]*chebyshev_U(args[0]-1,args[1])

NotImplementedError: derivative w.r.t. to the index is not supported yet

Also: doctest your exceptions.

[jdemeyer]

comment:60

I personally like to define chebyshev_T(-n,x) = chebyshev_T(n,x) and chebyshev_U(-n,x) = -chebyshev_U(n-2,x). This makes sense from a number-theoretic point of view and is also consistent with chebyshev_T(n,cos(x)) = cos(n*x) and chebyshev_U(n-1,cos(x)) = sin(n*x)/sin(x).

[jdemeyer]

comment:61

For real/complex fields, you should use

chebyshev_T(n,x) = ((x+sqrt(x^2-1))^n + (x-sqrt(x^2-1))^n)/2
chebyshev_U(n-1,x) = ((x+sqrt(x^2-1))^n - (x-sqrt(x^2-1))^n)/(2*sqrt(x^2-1))

[fredrik-johansson]

comment:62

The doubling recursion formulas should be better also for real/complex fields, I think: computing nth powers is basically the same amount of work, and it's best to avoid the square roots (especially for real |x| < 1).

[jdemeyer]

comment:63

Replying to @fredrik-johansson:

The doubling recursion formulas should be better also for real/complex fields, I think: computing nth powers is basically the same amount of work, and it's best to avoid the square roots (especially for real |x| < 1).

The matrix algorithm does seem more numerically stable (checked by using both algorithms inside RIF). So it's easy then, if there is one algorithm which is obviously the best.

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:64

I found an interesting paper on this topic: http://www.mathematik.uni-kassel.de/~koepf/cheby.pdf

Maybe this could give some input to the discussion.

It states that for an expanded representation the approach I have initally chosen (series expansion) should be the best for the symbolic evaluation, since this is somehow that what the user expects from ohter CAS. However, the proposed recursive/symbolic method would be interesting too. (p16f), since it gives a more compact form for large n. Maybe a switch for n >= 100?

Considering rational numbers as input, the recursion formula works best, due to this paper. This should be considered too.

@Matrix Multiplication: It's also my favorite, but I will compare first.

[fredrik-johansson]

comment:65

Here is an implementation of the divide-and-conquer algorithm that doesn't require caching (it only makes one recursive call). It might look even nicer if one rewrites it in iterative form. I think it's also equivalent to the code in pari. It should be faster than matrix powering by a constant factor, just like the analogous Fibonacci number algorithms.

def chebyshev_t(n, x):
    # returns (T(n,x), T(n-1,x)), or (T(n,x), _) if both=False
    def recur(n, x, both=False):
        if n == 0:
            return 1, x
        if n == 1:
            return x, 1
        a, b = recur((n+1)//2, x, both or n % 2)
        if n % 2 == 0:
            return 2*a^2 - 1, both and 2*a*b - x
        else:
            return 2*a*b - x, both and 2*b^2 - 1
    return recur(n, x, False)[0]

Come to think of it, it might even be useful to publicly expose a method that returns both T(n,x) and T(n-1,x) simultaneously.

Similar code for U (using same algorithm as Pari):

def chebyshev_u(n, x):
    def recur(n, x, both=False):
        if n == 0:
            return 1, both and 2*x
        if n == 1:
            return 2*x, both and 4*x^2-1
        a, b = recur((n-1)//2, x, True)
        if n % 2 == 0:
            return (b+a)*(b-a), both and 2*b*(x*b-a)
        else:
            return 2*a*(b-x*a), both and (b+a)*(b-a)
    return recur(n, x, False)[0]

Edit: streamlined the code.

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:66

I think your recursive implementation is very good. If you try to implement it iteratively, you have to consider some cases (current in row even/odd and next in row even/odd), and the code gets quite ugly in my opinion. And I think Knuth is right. One should prefer readable code over over optimzed "faster" code...

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:67

this would be a functioning iterative algorithm:

def chebyshev_t(n,x):

    if n == 0:
        return 1
    elif n == 1:
        return x
    elif n == 2:
        return 2*x**2-1
    else:
        T_c = x
        T_p = 2*x**2 -1

        for k in range(floor(log(n,2)),0,-1):
            T_p_old = T_p
            T_c_old = T_c
            if (n//2**(k-1)) % 2 == 0: # next is even
                T_p = 2*T_p_old*T_c_old - x
                T_c = 2*T_c_old**2 - 1
            elif (n//2**(k-1)) % 2  == 1: # next is odd
                T_p = 2*T_p_old**2 - 1
                T_c = 2*T_p_old*T_c_old - x

        # Cases for output
        if log(n - 1,2) in ZZ or log(n-2,2) in ZZ:
            return T_c

        elif n % 2 == 0:
            if n//2 % 2 == 0:                   
                return T_c
            else:
                return T_p
        elif n % 2 == 1:
            if n//2 % 2 == 0:                   
                return T_p
            else:
                return T_c

I made it shorter. What is preferable? recursive or iterative? Normaly iterative, but in this case it is longer due to the indices battles ...

Edit: I also measured time: recursive is faster, even if I do some optimization (e.g. xrange)

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:68

Replying to @jdemeyer:

Sorry to spoil the party, but this is a regression:

sage: K.<a> = NumberField(x^3-x-1)
sage: chebyshev_T(10^3,a)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-17-aa97c56dd147> in <module>()
----> 1 chebyshev_T(Integer(10)**Integer(3),a)

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.BuiltinFunction.__call__ (sage/symbolic/function.cpp:8126)()

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:5279)()

TypeError: cannot coerce arguments: no canonical coercion from Number Field in a with defining polynomial x^3 - x - 1 to Symbolic Ring

(and yes: Chebyshev polynomials are important in number theory)

I had now a look on these regressions: The problem originates from the fact, tat the ortho polys are now symbolic functions instead of simple function call in maxima. the symbolic functions have a lot more mechanisms concerning type coercions.

The only way around this would be a hack to chatch this coercion troubles, e.g.

def __call__(self,n,x):

    try:
        super(Func_chebyshev_T,self).__call__(n,x)
    except TypeError:
        return self._clenshaw_method_(n,x)

would work for that problem and some others. But how to catch such things correctly? Use a try catch, or make several is_instance checks? This could get a quite long list on coercions...

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

new fixes and patches

[jdemeyer]

comment:69

Attachment: trac_9706_new_clenshaw.patch.gz

Replying to @sagetrac-maldun:

would work for that problem and some others. But how to catch such things correctly?

I think there are only two cases: the "symbolic" case and the "algebraic" case. The latter means that we really consider the polynomial as a polynomial, not a symbolic function. In chebyshev(n,x), if either n or x is symbolic, we are in the symbolic case, otherwise we're in the algebraic case. In the algebraic case, the index n must be a concrete integer and we use the iterative algorithm.

[jdemeyer]

comment:70

maldun: I don't like all the log's in your approach (I don't think you need them), but otherwise I'm happy with either the recursive or iterative approach.

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:71

Replying to @jdemeyer:

Replying to @sagetrac-maldun:

would work for that problem and some others. But how to catch such things correctly?

I think there are only two cases: the "symbolic" case and the "algebraic" case. The latter means that we really consider the polynomial as a polynomial, not a symbolic function. In chebyshev(n,x), if either n or x is symbolic, we are in the symbolic case, otherwise we're in the algebraic case. In the algebraic case, the index n must be a concrete integer and we use the iterative algorithm.

Thank you for the input, that's a great idea!

Based on that I propose the following switch:

def __call__(self,n,x):
    if n in ZZ: # check if n is integer -> consider polynomial as algebraic structure
        self._eval_(n,x) # Let eval methode decide which is best
    else:
        super(OrthogonalPolynomial,self).__call__(n,x)

[jdemeyer]

comment:72

I think that chebyshev_T(1/2, 2) should raise a ValueError (or can we make sense of this?). So, in your code there should really be 3 cases: integer, symbolic and "something else" which is always an error.

So, I would do something like

def __call__(self,n,x):
    if is_Expression(n):
        return super(OrthogonalPolynomial, self).__call__(n,x)
    # We consider the polynomial really as a polynomial,
    # not a symbolic expression.
    try:
        n = ZZ(n)
    except StandardError:
        raise ValueError("Index for symbolic polynomials must be an integer")
    return self._eval_polynomial(n, x)

[jdemeyer]

Description changed:

--- 
+++ 
@@ -15,3 +15,4 @@

 * [attachment: trac_9706_chebyshev.patch](https://github.com/sagemath/sage-prod/files/10650399/trac_9706_chebyshev.patch.gz)
 * [attachment: trac_9706-review-ts.patch](https://github.com/sagemath/sage-prod/files/10650409/trac_9706-review-ts.patch.gz)
+* [attachment: trac_9706_new_clenshaw.patch​](https://github.com/sagemath/sage/files/ticket9706/92d8417838742f4fc31f28f66fb06843.gz)

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:74

Replying to @jdemeyer:

maldun: I don't like all the log's in your approach (I don't think you need them), but otherwise I'm happy with either the recursive or iterative approach.

I removed the logs, it is much faster now, but still slower than the recursive implementation about a constant factor of 2 (which is not much considering that we are talking about µs, but still)

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

Description changed:

--- 
+++ 
@@ -15,4 +15,3 @@

 * [attachment: trac_9706_chebyshev.patch](https://github.com/sagemath/sage-prod/files/10650399/trac_9706_chebyshev.patch.gz)
 * [attachment: trac_9706-review-ts.patch](https://github.com/sagemath/sage-prod/files/10650409/trac_9706-review-ts.patch.gz)
-* [attachment: trac_9706_new_clenshaw.patch​](https://github.com/sagemath/sage/files/ticket9706/92d8417838742f4fc31f28f66fb06843.gz)

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:75

Replying to @jdemeyer:

I think that chebyshev_T(1/2, 2) should raise a ValueError (or can we make sense of this?). So, in your code there should really be 3 cases: integer, symbolic and "something else" which is always an error.

So, I would do something like

def __call__(self,n,x):
    if is_Expression(n):
        return super(OrthogonalPolynomial, self).__call__(n,x)
    # We consider the polynomial really as a polynomial,
    # not a symbolic expression.
    try:
        n = ZZ(n)
    except StandardError:
        raise ValueError("Index for symbolic polynomials must be an integer")
    return self._eval_polynomial(n, x)

No there are several other cases to consider: You can also have complex and real input values if you consider a chebyshev polynomial as extension of the Hypergeometric function 1F2 like mpmath does: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/orthogonal.html#chebyt That's the reason for the message in the _derive_ method, since it would be theoretically possible to differentiate a chebyshev polynomial with respect to the index.

So the way I proposed makes sense, since this could be important for symbolic computation purposes where hypergeometric functions play an important role (eg. Zeilberger algorithm), and of course for analytical considerations.

[jdemeyer]

comment:76

OK, I agree.

What remains to do:

• support negative indices
• add many doctests (essentially, all the examples I mentioned on this ticket should become doctests)
• for ease of reviewing, fold all the patches into one patch.

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

Patch in a whole with corrections etc.

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:77

Attachment: trac_9706_chebyshev_patch_new.patch.gz

New patch attached. I incorporated all changes discussed.

@Pari I tried the evaluation with pari, but with Fredericks recursion, there is no gain in speed, due to the type checks beforehand. And the recursion in sage avoid type conversions.

[jdemeyer]

Description changed:

--- 
+++ 
@@ -13,5 +13,4 @@

 Apply:

-* [attachment: trac_9706_chebyshev.patch](https://github.com/sagemath/sage-prod/files/10650399/trac_9706_chebyshev.patch.gz)
-* [attachment: trac_9706-review-ts.patch](https://github.com/sagemath/sage-prod/files/10650409/trac_9706-review-ts.patch.gz)
+* [attachment: trac_9706_chebyshev_patch_new.patch](https://github.com/sagemath/sage-prod/files/10650398/trac_9706_chebyshev_patch_new.patch.gz)

[jdemeyer]

comment:78

meldun: please adjust the "apply" section in the ticket description so it's clear which patch(es) should be applied.

[jdemeyer]

Description changed:

--- 
+++ 
@@ -13,4 +13,7 @@

 Apply:

+* [attachment: trac_9706_chebyshev.patch](https://github.com/sagemath/sage-prod/files/10650399/trac_9706_chebyshev.patch.gz)
+* [attachment: trac_9706-review-ts.patch](https://github.com/sagemath/sage-prod/files/10650409/trac_9706-review-ts.patch.gz)
+* [attachment: trac_9706_new_clenshaw.patch​](https://github.com/sagemath/sage/files/ticket9706/92d8417838742f4fc31f28f66fb06843.gz)
 * [attachment: trac_9706_chebyshev_patch_new.patch](https://github.com/sagemath/sage-prod/files/10650398/trac_9706_chebyshev_patch_new.patch.gz)

[jdemeyer]

Attachment: trac_9706_chebyshev.patch.gz

cheby changes

[jdemeyer]

Description changed:

--- 
+++ 
@@ -14,6 +14,3 @@
 Apply:

 * [attachment: trac_9706_chebyshev.patch](https://github.com/sagemath/sage-prod/files/10650399/trac_9706_chebyshev.patch.gz)
-* [attachment: trac_9706-review-ts.patch](https://github.com/sagemath/sage-prod/files/10650409/trac_9706-review-ts.patch.gz)
-* [attachment: trac_9706_new_clenshaw.patch​](https://github.com/sagemath/sage/files/ticket9706/92d8417838742f4fc31f28f66fb06843.gz)
-* [attachment: trac_9706_chebyshev_patch_new.patch](https://github.com/sagemath/sage-prod/files/10650398/trac_9706_chebyshev_patch_new.patch.gz)

[jdemeyer]

comment:80

Folded the 4 patches.

[jdemeyer]

comment:82

What's the advantage of the _clenshaw_method_ over the recursive method? I see no need for the two different implementations and suggest to remove _clenshaw_method_.

[jdemeyer]

comment:83

More doctests are needed, this still doesn't work:

sage: chebyshev_T(1/2, 0)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-2-c142cc68c50b> in <module>()
----> 1 chebyshev_T(Integer(1)/Integer(2), Integer(0))

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in __call__(self, n, x)
    554             return self._eval_(n,x) # Let eval methode decide which is best
    555         else: # Consider OrthogonalPolynomial as symbol
--> 556             return super(OrthogonalPolynomial,self).__call__(n,x)
    557 
    558 

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.BuiltinFunction.__call__ (sage/symbolic/function.cpp:8126)()

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:5531)()

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _eval_(self, *args)
    493                 if not is_Expression(args[-1]):
    494                     try:
--> 495                         return self._evalf_(*args)
    496                     except AttributeError:
    497                         pass

/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _evalf_(self, *args, **kwds)
    641             precision = step_parent.prec()
    642         except AttributeError:
--> 643             precision = RR.prec()
    644 
    645         from sage.libs.mpmath.all import call as mpcall

NameError: global name 'RR' is not defined

[jdemeyer]

comment:84

Why do we have __call__ and _eval_ given they do essentially the same thing? Please keep in mind [comment:69] also in the _eval_() function (I really don't understand why it needs to be so complicated).

This is also still broken:

sage: parent(chebyshev_T(4, RIF(5)))
Real Field with 53 bits of precision

Please fix all issues that I mentioned and turn them into doctests.

[jdemeyer]

comment:85

For the REFERENCES, see http://sagemath.org/doc/developer/conventions.html#docstring-markup-with-rest-and-sphinx for the right syntax.

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

comment:86

@clenshaw_method: there is a difference. clenshaw method also applies a direct formula for small n and calls the recursive method else. The difference is that the recursive evaluation does not give an expanded representation of the polynomial, which is wanted for small n, because that was the standard till now and people expect this, especially if you are used to mathematica or maxima. Expanding huge expressions costs a lot of time, and this approach is much faster in that situation. Of course it is a matter of naming. But the reason why I have 2 methods, is to avoid too long code segments, and splitting them apart is better for readability. It also is important concerning other orthogonal polynomials.

@call & _eval_ : This convention is part of the BuiltinFunction structure. __call__ does all the stuff like coercions, transforming into a symbolic expression (e.g. if n is a symbolic value don't return a polynomial but hold the closed form.) _eval cares about evaluating the polynomial (e.g return a number if x is a number etc.) Look into the symbolic.function module for more details And eval is so complicated because there are several cases to consider: correct evaluation of symbolic expressions, numerical expressions and numpy arrays etc. This is also part of the BuiltinFunction structure. And you also have to keep in mind that this method should work for all ortho polys.

@bugs Sorry, during the patch merging process I had forgotten to apply a patch, which I'm now missing, since I work on different machines ... I will correct this tomorrow. It's annoying since I already had fixed it ...

[ac9ad401-3030-4fb0-957e-6c14f81e046a]

Attachment: trac_9706_bugfixes_and_doctests.patch.gz

incorporated things that were already been done ...