use FLINT to speed up Chebyshev T polynomial creation

[rws]

In #16670 the superiority of FLINT for creation of big Chebyshev-T-polynomials was confirmed. At T_10000 the speedup is already about 50x versus the present implementation. The issue is outsourced to this ticket.

[rws]

New commits:

​ba4e34d	16670: use FLINT for chebyshev_T polynomial
​df781ab	16670: fix previous patch: use argument parent gen
​e8dd233	16812: final fix and doctest

[jdemeyer]

Documentation is missing:

def chebyshev_T(unsigned long n, var='x'):
    """
    """

[jdemeyer]

Also the obvious question: can this be made to work also for chebyshev_U?

[jdemeyer]

I think the conditions n>10 and is_PolynomialRing(P) and P.base() == ZZ and P.ngens() == 1 and x == P.gen() should be generalized:

1. It should work for any n, at least negative n should be supported.

2. It should work for any ring of characteristic 0, not just ZZ.

3. Ideally, it should work for any kind of non-constant polynomial, not just generators of univariate polynomial rings.

Point 2 and 3 could be supported by first computing the polynomial in ZZ[x] and then substituting and/or changing the base ring as needed. I'm not saying this is needed on this ticket, but it would be the right thing to do.

[fredrik.johansson]

Why does flint_arith.chebyshev_T create an fmpz_poly, convert it to an array of ZZs, and then convert it back to a polynomial? Over ZZ[x], the roundtrip is probably more expensive than the actual computation. It should be enough to create a new Polynomial_integer_dense_flint instance and apply arith_chebyshev_t to its .__poly (the n > 10 condition can certainly be dropped).

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​7b7b8d4	16812: add documentationin libs.flint.arith
​d89d9ad	16812: same for chebyshev_U

[rws]

Replying to fredrik.johansson:

...It should be enough to create a new Polynomial_integer_dense_flint instance and apply arith_chebyshev_t to its .__poly

Well, I tried that but needed to cast polynomial.__poly to fmpz_poly_t which is not possible in Cython because the type is an array. OTOH cimporting the whole class Polynomial_integer_dense_flint declaration needs instances of the member functions defined. If I don't do that and just declare

cdef class Polynomial_integer_dense_flint(Polynomial):
    cdef fmpz_poly_t __poly

the function libs.flint.arith.chebyshev_T will finally compile but the result is not a rings.polynomial.Polynomial_integer_dense_flint and returning it and assigning it as such will bomb Sage.

I will now add a member function chebyshev_T to rings...Polynomial_integer_dense_flint to avoid all that hassle. My first forage into Cython which certainly presents itself as a fickle customer.

[fredrik.johansson]

Putting a method on Polynomial_integer_dense_flint sounds like the easiest way to do it.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​85a141b	16812: move chebyshev_t into class Polynomial_integer_dense_flint
​70214f5	16812: remove restriction on n

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​f5837ef	16812: generalize chebyshev_T() by substituting the argument and changing the base ring
​e5993da	16812: fix "except Exception: pass"
​6f3786b	16812: _eval_numpy_() was not called on numpy objects
​100622a	16812: padics doctests sent Sage into coma

[rws]

Replying to jdemeyer:

I think the conditions n>10 and is_PolynomialRing(P) and P.base() == ZZ and P.ngens() == 1 and x == P.gen() should be generalized:

1. It should work for any n, at least negative n should be supported.

2. It should work for any ring of characteristic 0, not just ZZ.

3. Ideally, it should work for any kind of non-constant polynomial, not just generators of univariate polynomial rings.

Point 2 and 3 could be supported by first computing the polynomial in ZZ[x] and then substituting and/or changing the base ring as needed. I'm not saying this is needed on this ticket, but it would be the right thing to do.

I have implemented this now and see no reason for the restriction in (2) to characteristic 0, as it seems to me that changing the base ring to, e.g., a finite field works nicely too.

In the process I uncovered at least two bugs.

I have removed the padics doctest for the moment because it was evaluated using the new FLINT code and the subsequent substitution was a bit too much. In this case the recursive evaulation is actually faster. But I have no idea at the moment when to select which algorithm, in order to not to fall in this trap.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​b032291	Merge branch 'develop' into t/16812/use_flint_to_speed_up_chebyshev_t_polynomial_creation
​61972b5	16812: fix doctest

[jdemeyer]

Replying to rws:

I have removed the padics doctest for the moment because it was evaluated using the new FLINT code and the subsequent substitution was a bit too much. In this case the recursive evaulation is actually faster. But I have no idea at the moment when to select which algorithm, in order to not to fall in this trap.

needs_work for exactly this reason.

[jdemeyer]

For the p-adics, how about this: compute the leading term of the resulting polynomial, which is 2^(n-1). If that result is zero, use the recursive algorithm. Otherwise, use the FLINT algorithm.

As a general rule, one should never remove doctests to "fix" something, unless you can show that there is a problem with the doctest itself.

[jdemeyer]

Replying to rws:

I have implemented this now and see no reason for the restriction in (2) to characteristic 0, as it seems to me that changing the base ring to, e.g., a finite field works nicely too.

It would certainly work in theory, but I just don't know if the FLINT algorithm is also faster in this case (I have not tried it).

[rws]

Replying to jdemeyer:

As a general rule, one should never remove doctests to "fix" something, unless you can show that there is a problem with the doctest itself.

I think you haven't noticed that I fixed a different doctest, not the p-adic one. Actually I was aware that the p-adic doctest should go in again, that's why I wrote "for the moment" earlier. I just wanted to have a version that passes all doctests (I failed to note that another doctest in a different file was bad).

[jdemeyer]

Replying to rws:

I think you haven't noticed that I fixed a different doctest, not the p-adic one.

I don't understand what you're trying to say. Fixing one doctest is not an excuse for breaking a different doctest.

Actually I was aware that the p-adic doctest should go in again, that's why I wrote "for the moment" earlier. I just wanted to have a version that passes all doctests.

OK fine, I understand that. But as long as not all those temporarily-removed doctests pass, a ticket should remain needs_work.

[rws]

On another note, the behaviour of the function up to this is unexpectedly complicated and not documented. Given integer n, if the second argument is symbolic then if n<32 a formula is applied else the recursive algorithm. With this ticket, additionally polynomials are handled with FLINT (modulo the refinements discussed).

Now, I have a patch ready for an algorithm keyword (flint/recursive). Also, I want to remove that n<32 bit together with the eval_formula method that really is no longer necessary now we have FLINT. It also leads to such inconsistencies as

            sage: var('n,x')
            (n, x)
            sage: chebyshev_T(5,x)
            16*x^5 - 20*x^3 + 5*x
            sage: chebyshev_T(64, x)
            2*(2*(2*(2*(2*(2*x^2 - 1)^2 - 1)^2 - 1)^2 - 1)^2 - 1)^2 - 1

which I would like to have either so or so, not both.

I think I propose to handle a symbolic argument exactly as a polynomial. No one expects such a nested result except if explicitly requested. It's also only faster because it's not expanded.

Do you have any objections?

[jdemeyer]

Replying to rws:

I think I propose to handle a symbolic argument exactly as a polynomial.

Do you have any objections?

Not at all, your suggestion sounds like a good idea.

[rws]

Another bug: special values apply only to integer n, as T(1/2, 0) shows, where the formula for x==0 in _eval_special_values_() gives complex values, while the correct result is 1/sqrt(2). So, assumptions have to be checked or the function restricted to n in ZZ.

Update: changed n==0 to x==0 in the above.

[rws]

I would also say there is a bug in mpmath because it cannot evaluate chebyshev_T(10^6, 0.1) without throwing a mpmath.NoConvergence while Sage's recursive algorithm is happy:

sage: chebyshev_T(10^6, 0.1, algorithm="recursive")
0.636384327171504
sage: timeit('chebyshev_T(10^6, 0.1, algorithm="recursive")')
625 loops, best of 3: 229 µs per loop

The result agrees with Wolfram to the 9th place which has 0.636384327247008190135495327481582817644932547700768253402752...

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​d21e70c	16812: implement algorithm keyword for eval()
​c4b8c17	16812: move examples to top class; doc adaptations
​a37b50e	16812: add links from top to classes/methods
​e5b19e4	Merge branch 'develop' into t/16812/use_flint_to_speed_up_chebyshev_t_polynomial_creation
​9fc5592	16812: we do not want _evalf_() called >100 times...
​830ff87	16812: we do not like "except Exception:"
​cc2e34f	16812: revamp Chebyshev polys

[rws]

This still needs some doctest cleanup.

[fredrik.johansson]

Replying to rws:

I would also say there is a bug in mpmath because it cannot evaluate chebyshev_T(10^6, 0.1) without throwing a mpmath.NoConvergence while Sage's recursive algorithm is happy:

Yep, mpmath currently doesn't optimize Chebyshev polynomials in any way for integer n (it always just calls 2F1). So there's no reason for Sage to use it for numerical evaluation.

[maldun]

Replying to rws:

On another note, the behaviour of the function up to this is unexpectedly complicated and not documented. Given integer n, if the second argument is symbolic then if n<32 a formula is applied else the recursive algorithm. With this ticket, additionally polynomials are handled with FLINT (modulo the refinements discussed).

Now, I have a patch ready for an algorithm keyword (flint/recursive). Also, I want to remove that n<32 bit together with the eval_formula method that really is no longer necessary now we have FLINT. It also leads to such inconsistencies as

            sage: var('n,x')
            (n, x)
            sage: chebyshev_T(5,x)
            16*x^5 - 20*x^3 + 5*x
            sage: chebyshev_T(64, x)
            2*(2*(2*(2*(2*(2*x^2 - 1)^2 - 1)^2 - 1)^2 - 1)^2 - 1)^2 - 1

which I would like to have either so or so, not both.

I think I propose to handle a symbolic argument exactly as a polynomial. No one expects such a nested result except if explicitly requested. It's also only faster because it's not expanded.

Do you have any objections?

Yes I have. There are users (like me) who want the possibility to switch between this representations, especially for large values of n. So please still allow to call the method explicitely.

[rws]

Replying to maldun

Yes I have. There are users (like me) who want the possibility to switch between this representations, especially for large values of n. So please still allow to call the method explicitely.

Exactly this is implemented. You can call the method, or you can say algorithm=recursive.

[maldun]

Replying to rws:

Replying to maldun

Yes I have. There are users (like me) who want the possibility to switch between this representations, especially for large values of n. So please still allow to call the method explicitely.

Exactly this is implemented. You can call the method, or you can say algorithm=recursive.

Good to hear! And thank you for your work! As mentioned in the other ticket: Don't hesitate to ask for help.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​90090fd

16812: doctest every execution path in chebyshev_T

[rws]

The FLINT routines are now called on polynomials with exact base ring and used for symbolic input. Everything else gets the recursive treatment which is really fast if you don't have to expand afterwards.

Please review.

[maldun]

Typing

chebyshev_T(1000000,x)

in the old version gives a proper result in 7.52 The exact same input gives me full memory (8GB on my machine) and about 1GB swap + a frozen system in the new version. Higher input gives a an error of flint (Memory allocation error)

Some switch would be nice. I guess the best measure would be the break even point on memory consumption.

Edit: I think a good point would be about order 5000-10000: It is reasonable from view of speed/memory and an explicit representation would be too cluttered anyway

[maldun]

New commits:

​59d04f5

Trac 16812: added a switch for evaluation between recursive and flint

[maldun]

I proposed a patch for that

[rws]

Replying to maldun:

Typing

chebyshev_T(1000000,x)

in the old version gives a proper result in 7.52 The exact same input gives me full memory (8GB on my machine) and about 1GB swap + a frozen system in the new version.

That result may be proper for a symbolic x but the times are completely different for ring elements because expansion, which is avoided with symbolic x, happens automatically. This speed difference for polynomial ring elements actually is the reason for this ticket (see description).

Higher input gives a an error of flint (Memory allocation error)

Some switch would be nice. I guess the best measure would be the break even point on memory consumption.

Edit: I think a good point would be about order 5000-10000: It is reasonable from view of speed/memory and an explicit representation would be too cluttered anyway

I don't think Sage has to anticipate user memory size, I can go happily to 150.000 with 8GB:

sage: R.<x> = ZZ[]
sage: timeit('chebyshev_T(10^5,x)')
5 loops, best of 3: 454 ms per loop
sage: timeit('chebyshev_T(150000,x)')
5 loops, best of 3: 1.01 s per loop

Who knows what will be in ten years? Note also that noone displays such results, and that's good because the output process takes near all memory and CPU, not FLINT.

Now, it might have been a bad idea by me to use FLINT on symbolic input but I thought users might appreciate unconvoluted polynomials as output. OTOH, they might get another motivation from this to use rings if they want straight output of polys. So I'll roll back the symbolic involvement, no limit is needed, and Bob's your uncle.

[maldun]

Replying to rws:

Who knows what will be in ten years? Note also that noone displays such results, and that's good because the output process takes near all memory and CPU, not FLINT.

Now, it might have been a bad idea by me to use FLINT on symbolic input but I thought users might appreciate unconvoluted polynomials as output. OTOH, they might get another motivation from this to use rings if they want straight output of polys. So I'll roll back the symbolic involvement, no limit is needed, and Bob's your uncle.

Than the solution is clear: Make a type check on the input type: If symbolic + higher order use recursion, if Polynomial Ring use flint

It's not about availabilty of RAM: no one wants that a Chebyshev polynomal consumes 1GB of RAM, and no one wants a freezing system, especially if there is a compact represention which only needs the logarithmical amount of time and memory. And also: You don't have much use with an expanded symbolic Cheby poly of order 1e+5, only a cluttered screen, whereas the unexpanded style still is human readable.

Always consider: The main input method should provide some more or less smart default behaviour

On another note: Encourage people to use polynomial ring by putting proper examples in the docstring

[maldun]

sage: T = chebyshev_T(100000000000,x)
Exception (FLINT memory_manager). Unable to allocate memory.
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-1-bef10d0acafb> in <module>()
----> 1 T = chebyshev_T(Integer(100000000000),x)

/home/maldun/sage/sage-6.3/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in __call__(self, n, *args, **kwds)
    624         # If n is an integer: consider the polynomial as an algebraic (not symbolic) object
    625         if n in ZZ and not kwds.get('hold', False):
--> 626             ret = self._eval_(n, *args, **kwds)
    627             if ret is not None:
    628                 return ret

/home/maldun/sage/sage-6.3/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _eval_(self, n, x, *args, **kwds)
    679             if ((is_Polynomial(x) and x.parent().base_ring().is_exact()) or
    680                 (is_Expression(x) and (not x.is_numeric()) and ((n > 10000) or (n < -10000)))):
--> 681                 return self.eval_algebraic(n, x)
    682             else:
    683                 return self.eval_recursive(n, x)

/home/maldun/sage/sage-6.3/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in eval_algebraic(self, n, arg)
    918                 n = -n
    919             R = PolynomialRing(ZZ, 'x')
--> 920             pol = Polynomial_integer_dense_flint.chebyshev_T(n, R, R.gen())
    921             if is_PolynomialRing(P):
    922                 pol = pol.subs({pol.parent().gen():arg})

/home/maldun/sage/sage-6.3/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_integer_dense_flint.so in sage.rings.polynomial.polynomial_integer_dense_flint.Polynomial_integer_dense_flint.chebyshev_T (build/cythonized/sage/rings/polynomial/polynomial_integer_dense_flint.cpp:14398)()

/home/maldun/sage/sage-6.3/local/lib/python2.7/site-packages/sage/ext/c_lib.so in sage.ext.c_lib.sig_raise_exception (build/cythonized/sage/ext/c_lib.c:1040)()

RuntimeError: Aborted

Recursive version works without abort

This kills my system:

sage: R.<x> = ZZ[]
sage: timeit('chebyshev_T(10^7,x)')

so it is not only the use of symbolic expressions

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​e5c0a83	Trac 16812: Make switch in default evaluation method
​9e6001c	Merge branch 'cheby3' into cheby
​f02d699	Trac 16812: removed non default switch mechanism
​6165cbb	Trac 16812: Make another proposal for a switch mechanism

[maldun]

I propose the following switch mechanism: If the input is symbolic and deg > 10000 return the recursive result else use the regular algorithm

For other input the user has to take responsibility Of course polynomials don't work good with the python recursion, flint does a better job. For symbolic expressions for very high degree the recursion algorithm does a quite good job

[rws]

I fixed your doctest. Let's see what the reviewer/s say.

---

New commits:

​2ee34da

16812: fix doctest

[rws]

sage -t --long src/sage/graphs/bipartite_graph.py  # 1 doctest failed
sage -t --long src/sage/functions/orthogonal_polys.py  # 11 doctests failed

[rws]

It maybe better to just make the Chebys GinacFunctions and use Flint from within Pynac.

[rws]

While other orthogonal polynomial functions (as well as almost all other symbolic functions) are simply BuiltinFunctions or GinacFunctions this ticket suffers from the fact that the Cheby functions in Sage allow the algorithm keyword (by using a __call__ method in the superclass). The way such polymorphisms are resolved in other functions is by having the interface (the callable global function chebyshev_T()) dispatching on the keyword. This refactoring step (#24554) should be done before all items of this ticket can be implemented.

[rws]

We will convert the Chebyshevs to GinacFunctions and call Flint from Pynac, which also will be the fastest way of creating the expression polynomial.

[rws]

New commits:

​907aff6	24497: pkg version/chksum
​d085587	24497: doctest fixes
​810f238	Merge commit 'd085587a0f4319f03a830c39eb5b3d83836f18dd' of git://trac.sagemath.org/sage into HEAD
​5c343ab	24554: refactor chebyshev_T
​23190a6	24554: refactor chebyshev_U
​2582330	24554: remove ChebyshevFunction
​96e4e39	Merge branch 'u/rws/remove___call_____usage_in_chebyshev_functions' of git://trac.sagemath.org/sage into t16812
​d7a055d	16812: make Chebyshev T/U GinacFunctions

[rws]

Branch sems corrupt.