symbolic Laguerre / associated Laguerre polynomials

[rws]

Not only is laguerre(n,x) not symbolic, a naive implementation is already 2x as fast at n=100 and 10x at n=1000:

R.<x> = PolynomialRing(QQ, 'x')
def lag(n):
    return R([binomial(n,k)*(-1)^k/factorial(k) for k in xrange(n+1)])

• ​https://en.wikipedia.org/wiki/Laguerre_polynomials
• ​http://mathworld.wolfram.com/LaguerrePolynomial.html
• ​http://mathworld.wolfram.com/AssociatedLaguerrePolynomial.html

[rws]

New commits:

​6d10729	Simplify numerical evaluation of BuiltinFunctions
​b6e1ed4	Merge remote-tracking branches 'origin/u/jdemeyer/ticket/17131' and 'origin/u/jdemeyer/ticket/17133' into ticket/17130
​04b2b37	Merge branch 'u/jdemeyer/ticket/17130' of trac.sagemath.org:sage into t/17151/symbolic_laguerre___associated_laguerre_polynomials
​3d3d1fa	17151: skeleton impl.
​b3edc51	17151: implement symbolic Laguerre pol.
​c8beac2	17151: implement laguerre(n,x)
​1626a8b	17151: implement gen_laguerre

[git]

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

​f727de5

17151: reduce evalf logic further; fix merge conflict

[kcrisman]

For some reason this branch is red.

[git]

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

​6d05ac6

Merge branch 'develop' into t/17151/symbolic_laguerre___associated_laguerre_polynomials

[rws]

sage -t --long src/sage/symbolic/expression.pyx  # 1 doctest failed
sage -t --long src/sage/matrix/matrix2.pyx  # 1 doctest failed
sage -t --long src/sage/graphs/bipartite_graph.py  # 1 doctest failed

[rws]

Fixed. Squashed it all into one commit.

EDIT: Or not.

---

New commits:

​80f136b	Merge branch 'develop' into t/17151/symbolic_laguerre___associated_laguerre_polynomials
​ea2ff9b	17151: symbolic Laguerre / associated Laguerre polynomials

[mmezzarobba]

sage: laguerre(-9, 2.)
0.000000000000000

sage: laguerre(-9, 2).n()
1566.22186244286

[mmezzarobba]

Are the three separate cases in _pol_laguerre() necessary? Couldn't one just do something like

R = PolynomialRing(QQ, 'x')
pol = R([binomial(n,k)*(-1)**k/factorial(k) for k in xrange(n+1)])
return pol(x)

no matter what x is?

[mmezzarobba]

sage: gen_laguerre(3, 1, x+1).expand()
-1/6*x^3 + 3/2*x^2 - 5/2*x - 1/6
sage: gen_laguerre(0, 1, x+1).expand()
...
AttributeError: 'sage.rings.integer.Integer' object has no attribute 'expand'

[mmezzarobba]

Otherwise looks good to me overall...

[rws]

Replying to mmezzarobba:

sage: gen_laguerre(0, 1, x+1).expand()
...
AttributeError: 'sage.rings.integer.Integer' object has no attribute 'expand'

Well that is ZZ(1) and cannot be expanded. I don't think we should return SR(1).

[mmezzarobba]

Replying to rws:

Replying to mmezzarobba:

sage: gen_laguerre(0, 1, x+1).expand()
...
AttributeError: 'sage.rings.integer.Integer' object has no attribute 'expand'

Well that is ZZ(1) and cannot be expanded. I don't think we should return SR(1).

Why? Having similar calls to the same function return elements of different parents (or even values of different types) is rarely a good idea... That being said, I don't think we should return SR(1) in all cases, but perhaps something like pushout(QQ, parent(x)).

[rws]

Replying to mmezzarobba:

Are the three separate cases in _pol_laguerre() necessary? Couldn't one just do something like

R = PolynomialRing(QQ, 'x')
pol = R([binomial(n,k)*(-1)**k/factorial(k) for k in xrange(n+1)])
return pol(x)

no matter what x is?

No, because (apart from speed reasons) the R([...]) construction is not available in all parents. However, when investigating this I found that the argument is coerced into SR by the superclasses BuiltinFunction/Function, so it must be unwrapped for the real parent. OTOH, for simplicity the branch with sum would be a catch-all. Let's see if going polynomial really is faster and, if not, do only summing for all parents.

[git]

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

​ae4c8e4	Merge branch 'develop' into t/17151/17151
​438123e	17151: work around mpmath issue 307
​009c028	17151: improve _pol_laguerre

[rws]

Let's see if going polynomial really is faster and, if not, do only summing for all parents.

Here's the catch: Using R(...) for polynomial creation makes for 5x speedup with L(1000,x) and 10x with L(1000,x^2+x+1) if x is a polynomial generator vs. a symbol. OTOH, summing over monomials is faster with sum and this is the only way to get a result for other rings.

I don't thing we should return SR(1) in all cases, but perhaps something like pushout(QQ, parent(x)).

There is a problem: even if I return SR(ZZ(1)), it gets converted back somewhere to Integer(1). This will have to be addressed in another ticket.

[rws]

There is a problem: even if I return SR(ZZ(1)), it gets converted back somewhere to Integer(1). This will have to be addressed in another ticket.

See #17953.

[rws]

Squashed it all into one commit. Please note I'll rename OrthogonalPolynomial to OrthogonalFunction in the other ticket, so I depend on this being merged first. In the end the type problem required a one liner in #17953 and didn't require any forced result conversion to SR.

---

New commits:

​a99ad55	17151: symbolic Laguerre / associated Laguerre polynomials
​f0fe6f4	17953: any symbolic function arg prevents forced result conversion to numeric
​351b5c6	Merge branch 'u/rws/inconsistency_in_returned_type_of_builtinfunction_result' of trac.sagemath.org:sage into tmp1

[rws]

Passes all patchbot tests.

[mmezzarobba]

Replying to rws:

Replying to mmezzarobba:

Are the three separate cases in _pol_laguerre() necessary? Couldn't one just do something like

R = PolynomialRing(QQ, 'x')
pol = R([binomial(n,k)*(-1)**k/factorial(k) for k in xrange(n+1)])
return pol(x)

no matter what x is?

No, because (apart from speed reasons) the R([...]) construction is not available in all parents.

I don't understand what you mean, sorry: I just tried the exact code that I suggested, and it works on all examples in the docstring, except that the result comes out in Horner form (which we could perhaps change if nothing else depends on it or somehow get or get around otherwise). What parents do you have in mind?

[mmezzarobba]

I also don't like that

sage: laguerre(3, polygen(QQ)).parent()
Symbolic Ring

while for instance

sage: chebyshev_T(3, polygen(QQ)).parent()
Univariate Polynomial Ring in x over Rational Field
sage: pow(polygen(QQ), 3).parent()
Univariate Polynomial Ring in x over Rational Field

[mmezzarobba]

Not a bug, but would be nice to have: the implementation doesn't seem to know that L(n, 0) = 1 and more generally L(n, α, 0) = binomial(n+α,n).

[rws]

Replying to mmezzarobba:

I also don't like that

sage: laguerre(3, polygen(QQ)).parent()
Symbolic Ring

Equivalently to #16813 this behaviour can only be changed when we have #18832. EDIT: pow is not a BuiltinFunction and chebyshev_T uses its own __call__ which is frowned upon. I would advise to not make this ticket dependent on the feature you suggest but to open a followup ticket.

[rws]

With the new branch all other comments were addressed.

---

New commits:

​6631c77

17151: symbolic Laguerre / associated Laguerre polynomials

[chapoton]

needs rebase, does not apply

[chapoton]

rebased

---

New commits:

​f0d809f

Merge branch 'u/rws/17151-2' into 6.9.b4

[mmezzarobba]

As explained in one of my previous comments, I still don't understand why you don't like the idea of using QQ['x'] to compute the coefficients in _pol_laguerre() (especially since #18282 has been merged and the problem with symbolic results being computed in Horner form no longer exists). But that's something that can be fixed later and this ticket has languished long enough already.

Thanks for working on it!

[rws]

Thanks for the review.