Opened 5 years ago

Closed 4 years ago

#17151 closed enhancement (fixed)

symbolic Laguerre / associated Laguerre polynomials

Reported by: rws Owned by:
Priority: major Milestone: sage-6.8
Component: symbolics Keywords: special, function, holonomic, orthogonal
Cc: Merged in:
Authors: Ralf Stephan Reviewers: Marc Mezzarobba
Report Upstream: N/A Work issues:
Branch: f0d809f (Commits) Commit: f0d809ffdc8d1342aabbc12c95bb5d5c15b1648b
Dependencies: #17953 Stopgaps:

Description (last modified by 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)])

Change History (36)

comment:1 Changed 5 years ago by rws

  • Description modified (diff)

comment:2 Changed 5 years ago by rws

  • Branch set to u/rws/symbolic_laguerre___associated_laguerre_polynomials

comment:3 Changed 5 years ago by rws

  • Commit set to 1626a8b208a1c2641d05cf8c76d9df86714fb4f0
  • Status changed from new to needs_review

New commits:

6d10729Simplify numerical evaluation of BuiltinFunctions
b6e1ed4Merge remote-tracking branches 'origin/u/jdemeyer/ticket/17131' and 'origin/u/jdemeyer/ticket/17133' into ticket/17130
04b2b37Merge branch 'u/jdemeyer/ticket/17130' of trac.sagemath.org:sage into t/17151/symbolic_laguerre___associated_laguerre_polynomials
3d3d1fa17151: skeleton impl.
b3edc5117151: implement symbolic Laguerre pol.
c8beac217151: implement laguerre(n,x)
1626a8b17151: implement gen_laguerre

comment:4 Changed 5 years ago by rws

  • Authors set to Ralf Stephan

comment:5 Changed 5 years ago by git

  • Commit changed from 1626a8b208a1c2641d05cf8c76d9df86714fb4f0 to f727de55a3d6ea92be9c70606f57f81795b4c73a

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

f727de517151: reduce evalf logic further; fix merge conflict

comment:6 Changed 5 years ago by kcrisman

  • Status changed from needs_review to needs_work

For some reason this branch is red.

comment:7 Changed 5 years ago by git

  • Commit changed from f727de55a3d6ea92be9c70606f57f81795b4c73a to 6d05ac6987afb4423a6978888e5566b9730fa793

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

6d05ac6Merge branch 'develop' into t/17151/symbolic_laguerre___associated_laguerre_polynomials

comment:8 Changed 5 years ago by rws

  • Status changed from needs_work to needs_review

comment:9 Changed 5 years ago by rws

  • Status changed from needs_review to needs_work
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

comment:10 Changed 5 years ago by rws

  • Branch changed from u/rws/symbolic_laguerre___associated_laguerre_polynomials to u/rws/17151

comment:11 Changed 5 years ago by rws

  • Commit changed from 6d05ac6987afb4423a6978888e5566b9730fa793 to ea2ff9b11a520edd58b1db527ef41eed6b0e7f85
  • Milestone changed from sage-6.4 to sage-6.6
  • Status changed from needs_work to needs_review

Fixed. Squashed it all into one commit.

EDIT: Or not.


New commits:

80f136bMerge branch 'develop' into t/17151/symbolic_laguerre___associated_laguerre_polynomials
ea2ff9b17151: symbolic Laguerre / associated Laguerre polynomials
Last edited 5 years ago by rws (previous) (diff)

comment:12 Changed 5 years ago by mmezzarobba

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

comment:13 follow-up: Changed 5 years ago by 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?

comment:14 follow-up: Changed 5 years ago by 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'

comment:15 Changed 5 years ago by mmezzarobba

  • Status changed from needs_review to needs_work

Otherwise looks good to me overall...

comment:16 in reply to: ↑ 14 ; follow-up: Changed 5 years ago by 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).

comment:17 in reply to: ↑ 16 Changed 5 years ago by 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)).

Last edited 5 years ago by mmezzarobba (previous) (diff)

comment:18 in reply to: ↑ 13 ; follow-ups: Changed 5 years ago by 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.

comment:19 Changed 5 years ago by git

  • Commit changed from ea2ff9b11a520edd58b1db527ef41eed6b0e7f85 to 009c0289def83e6e3c2d0911bd44b275cf0d635a

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

ae4c8e4Merge branch 'develop' into t/17151/17151
438123e17151: work around mpmath issue 307
009c02817151: improve _pol_laguerre

comment:20 in reply to: ↑ 18 ; follow-up: Changed 5 years ago by rws

  • Status changed from needs_work to needs_review

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.

comment:21 in reply to: ↑ 20 Changed 5 years ago by rws

  • Dependencies set to #17953

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.

comment:22 Changed 5 years ago by rws

  • Branch changed from u/rws/17151 to u/rws/17151-1

comment:23 Changed 5 years ago by rws

  • Commit changed from 009c0289def83e6e3c2d0911bd44b275cf0d635a to 351b5c69625489c7c12139436b193a3735fd429f

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:

a99ad5517151: symbolic Laguerre / associated Laguerre polynomials
f0fe6f417953: any symbolic function arg prevents forced result conversion to numeric
351b5c6Merge branch 'u/rws/inconsistency_in_returned_type_of_builtinfunction_result' of trac.sagemath.org:sage into tmp1

comment:24 Changed 5 years ago by rws

  • Milestone changed from sage-6.6 to sage-6.7

Passes all patchbot tests.

comment:25 in reply to: ↑ 18 ; follow-ups: Changed 5 years ago by 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?

comment:26 follow-up: Changed 5 years ago by 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

comment:27 Changed 5 years ago by 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).

comment:28 in reply to: ↑ 25 Changed 5 years ago by rws

  • Status changed from needs_review to needs_work

comment:29 in reply to: ↑ 26 Changed 4 years ago by 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.

Last edited 4 years ago by rws (previous) (diff)

comment:30 Changed 4 years ago by rws

  • Branch changed from u/rws/17151-1 to u/rws/17151-2

comment:31 Changed 4 years ago by rws

  • Commit changed from 351b5c69625489c7c12139436b193a3735fd429f to 6631c77809f8e1025105508232ef2615f363e358
  • Milestone changed from sage-6.7 to sage-6.8
  • Status changed from needs_work to needs_review

With the new branch all other comments were addressed.


New commits:

6631c7717151: symbolic Laguerre / associated Laguerre polynomials

comment:32 Changed 4 years ago by chapoton

  • Status changed from needs_review to needs_work

needs rebase, does not apply

comment:33 Changed 4 years ago by chapoton

  • Branch changed from u/rws/17151-2 to public/ticket/17151
  • Commit changed from 6631c77809f8e1025105508232ef2615f363e358 to f0d809ffdc8d1342aabbc12c95bb5d5c15b1648b
  • Status changed from needs_work to needs_review

rebased


New commits:

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

comment:34 in reply to: ↑ 25 Changed 4 years ago by mmezzarobba

  • Reviewers set to Marc Mezzarobba
  • Status changed from needs_review to positive_review

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!

comment:35 Changed 4 years ago by rws

Thanks for the review.

comment:36 Changed 4 years ago by vbraun

  • Branch changed from public/ticket/17151 to f0d809ffdc8d1342aabbc12c95bb5d5c15b1648b
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.