use FLINT to speed up Chebyshev T polynomial creation

[rwst]

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.

Depends on #24554
Depends on #24668

Component: symbolics

Keywords: flint, speedup

Author: Ralf Stephan

Branch/Commit: u/rws/16812 @ d7a055d

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

[rwst]

Branch: u/rws/use_flint_to_speed_up_chebyshev_t_polynomial_creation

[rwst]

Commit: e8dd233

[rwst]

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

[rwst]

Changed keywords from none to flint, speedup

[jdemeyer]

comment:3

Documentation is missing:

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

[jdemeyer]

comment:4

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

[jdemeyer]

comment:5

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]

comment:6

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).

[sagetrac-git]

Changed commit from e8dd233 to d89d9ad

[sagetrac-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

[rwst]

comment:8

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]

comment:9

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

[sagetrac-git]

Changed commit from d89d9ad to 70214f5

[sagetrac-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

[sagetrac-git]

Changed commit from 70214f5 to 100622a

[sagetrac-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

[rwst]

comment:12

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.

[sagetrac-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

[sagetrac-git]

Changed commit from 100622a to 61972b5

[jdemeyer]

comment:15

Replying to @rwst:

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.

[rwst]

Changed commit from 6165cbb to 2ee34da

[rwst]

comment:41

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

---

New commits:

2ee34da

16812: fix doctest

[rwst]

comment:42

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

[rwst]

Dependencies: #17531

[rwst]

comment:44

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

[rwst]

comment:45

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.

[rwst]

Changed dependencies from #17531 to #24554

[rwst]

Changed commit from 2ee34da to none

[rwst]

comment:46

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

[rwst]

Changed branch from u/rws/use_flint_to_speed_up_chebyshev_t_polynomial_creation to none

[rwst]

Branch: u/rws/16812

[rwst]

Commit: d7a055d

[rwst]

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

[rwst]

Changed dependencies from #24554 to #24554, pynac-0.7.16

[rwst]

Changed dependencies from #24554, pynac-0.7.16 to #24554, #24668

[rwst]

comment:50

Branch sems corrupt.