Opened 3 years ago

# Scale legendre_P to [a,b]

Reported by: Owned by: mjo burcin major sage-6.10 symbolics Michael Orlitzky, Frédéric Chapoton Francis Clarke, Karl-Dieter Crisman N/A u/chapoton/13720 (Commits) 0843342b41d1ce501cab9b84c89e2ef1d84e5dba

The Legendre polynomials, returned by legendre_P(), of the first kind are orthogonal over [-1,1] and are normalized to have value +-1 at the endpoints.

When solving least-squares problems, it's convenient to be able to construct them over an arbitrary interval [a,b].

### comment:1 Changed 3 years ago by mjo

• Authors set to Michael Orlitzky
• Status changed from new to needs_review

### comment:2 follow-up: ↓ 3 Changed 3 years ago by fwclarke

• Status changed from needs_review to needs_work
• The code seems over-complicated to me. Replacing the last few lines with
```    R = PolynomialRing(ZZ, 't')
f = R([(-1)^(n - k)*binomial(n, k)*binomial(n + k, k) for k in range(n + 1)])
return f((x - a)/(b - a))
```
is much simpler and is significantly faster. Moreover, this makes expressions such as
```legendre_P(4, Zmod(5), 3, 6)
```

evaluate correctly; at it stands the patch yields an

```ArithmeticError: 0^0 is undefined.
```
• I don't see why a and b need to be real numbers. Actually they could be polynomial variables, giving:
```sage: R.<r,s,t> = QQ[]
sage: legendre_P(2, t, r, s)
6*((r - t)/(r - s) - 1)*(r - t)/(r - s) + 1
```
• I don't see the need to use bool in doctests such as
```bool(legendre_P(0, x, 0, 1) == p0)
```
• It's not clear to me why back-quotes are used in some of the error strings:
```        raise TypeError('`n` must be a natural number')
```

but

```        raise ValueError('n must be nonnegative')
```

### comment:3 in reply to: ↑ 2 Changed 3 years ago by mjo

• The code seems over-complicated to me. Replacing the last few lines with
```    R = PolynomialRing(ZZ, 't')
f = R([(-1)^(n - k)*binomial(n, k)*binomial(n + k, k) for k in range(n + 1)])
return f((x - a)/(b - a))
```
is much simpler and is significantly faster. Moreover, this makes expressions such as
```legendre_P(4, Zmod(5), 3, 6)
```

evaluate correctly; at it stands the patch yields an

```ArithmeticError: 0^0 is undefined.
```

I tried this, but it's producing the wrong results. For example,

```sage: legendre_P(1, x)
-1/2*x - 5/2
```

That said, there are two reasons it might seem over-complicated.

The first is that I was very careful not to break any existing doctests, even though they're doing some crazy things. I tried to comment each of these hurdles.

The second is that I wanted to be clear about what was happening. The c(m) and g(m) functions are as defined in A&S. So the final return statement is exactly the form given in the reference. The use of the affine transform phi() allows me to retain that form, and also provides some rationale for the result. I think this makes it (more) clear that we're composing the standard P() over [-1,1] with an affine transform, which will intuitively preserve orthogonality.

• I don't see why a and b need to be real numbers. Actually they could be polynomial variables, giving:
```sage: R.<r,s,t> = QQ[]
sage: legendre_P(2, t, r, s)
6*((r - t)/(r - s) - 1)*(r - t)/(r - s) + 1
```

I'm happy to change this if it can be done without breaking the existing doctests. I'll wait to see what you say about the first thing.

• I don't see the need to use bool in doctests such as
```bool(legendre_P(0, x, 0, 1) == p0)
```

Often, the == operator will return a symbolic equality rather than True/False?. This is true even in this simple case:

```sage: p0 = 1
sage: legendre_P(0, x, 0, 1) == p0
1 == 1
```

If you cast the symbolic equality to bool, you get,

• True, if Sage can convince itself that equality holds.
• False, if Sage can "prove" inequality.
• False, if it's not sure.

So it should be safe to test for True.

• It's not clear to me why back-quotes are used in some of the error strings:
```        raise TypeError('`n` must be a natural number')
```

but

```        raise ValueError('n must be nonnegative')
```

I was just careless with that one, I'll fix it. I think backticks or single quotes are needed around 'a' and 'b', otherwise it reads weird. But 'n' at the beginning is fine without it.

### comment:4 follow-up: ↓ 5 Changed 3 years ago by mjo

• Status changed from needs_work to needs_review

I've fixed the variable quoting in the new patch, and allowed the endpoints a,b to be symbolic. There's also a new test to make sure that the symbolic endpoints work as expected.

I haven't heard back, so maybe you bought my reasoning for the affine transform =)

The bool() casts do have to stay, otherwise the tests just don't work.

Thanks for taking a look.

### comment:5 in reply to: ↑ 4 ; follow-up: ↓ 6 Changed 3 years ago by fwclarke

• Status changed from needs_review to needs_work

I've fixed the variable quoting in the new patch, and allowed the endpoints a,b to be symbolic. There's also a new test to make sure that the symbolic endpoints work as expected.

Why force them to be symbolic? Better, surely, to allow anything that will evaluate and let coercion handle it.

I haven't heard back, so maybe you bought my reasoning for the affine transform =)

No, just busy with other things. There is an affine transform in what I wrote, but it's a transform of the interval [a, b] onto [0, 1] rather than [-1, 1]. This is simpler, and the Legendre polynomials on the interval [0, 1] have a simple form.

I tried this, but it's producing the wrong results. For example,

```sage: legendre_P(1, x)
-1/2*x - 5/2
```

This is what results from using ^ rather than ** in python. I'm sorry my proposed code had only been testing in Sage, and therefore was being prepared. This dealt with, my version does pass all the doctests, is faster when if the variable is a polynomial generator, and the error

```sage: legendre_P(4, Zmod(5)(1), 3, 6)
...
ArithmeticError: 0^0 is undefined.
```

is eliminated. Apologies for mistyping this before.

The bool() casts do have to stay, otherwise the tests just don't work.

Another way of doing this is to introduce a polynomial variable. Thus

```sage: R.<x> = QQ[]
sage: legendre_P(1, x) == x
True
```

I know very little about the symbolic side of Sage, but I do find things like

```sage: SR(1) == SR(1)
1 == 1
```

very strange, and confusing when compared with

```sage: 1 == 1
True
```

### comment:6 in reply to: ↑ 5 Changed 3 years ago by mjo

Why force them to be symbolic? Better, surely, to allow anything that will evaluate and let coercion handle it.

I removed the RR check, so you can actually pass in anything you want now. I'm not convinced this isn't a foot-gun, but it's still a clear improvement over what we have now so I'm happy with it.

This is what results from using ^ rather than ** in python.

Ah, yes, I should have caught that, sorry.

I'm sorry my proposed code had only been testing in Sage, and therefore was being prepared. This dealt with, my version does pass all the doctests, is faster when if the variable is a polynomial generator, and the error

```sage: legendre_P(4, Zmod(5)(1), 3, 6)
...
ArithmeticError: 0^0 is undefined.
```

is eliminated. Apologies for mistyping this before.

I believe this is another bug, elsewhere. I've set out to fix it and am trapped in a rabbit hole at the moment, about four bugs from here. Once that's sorted out, I'll do some experiments.

The bool() casts do have to stay, otherwise the tests just don't work.

Another way of doing this is to introduce a polynomial variable. Thus

```sage: R.<x> = QQ[]
sage: legendre_P(1, x) == x
True
```

I know very little about the symbolic side of Sage, but I do find things like

```sage: SR(1) == SR(1)
1 == 1
```

very strange, and confusing when compared with

```sage: 1 == 1
True
```

You'll get no argument from me there. I do prefer to test with the default symbol x though just because that's what most people will use.

### comment:7 follow-up: ↓ 8 Changed 3 years ago by kcrisman

I think that people on this ticket may find #9706 interesting, though I imagine it has bitrotted somewhat and Burcin had a very long list of comments for it.

### comment:8 in reply to: ↑ 7 Changed 3 years ago by mjo

• Status changed from needs_work to needs_review

I think that people on this ticket may find #9706 interesting, though I imagine it has bitrotted somewhat and Burcin had a very long list of comments for it.

That's too big to ever get merged... I would break out the functions into individual tickets where possible.

Anyway, I was waiting here on #13786 which fixes the ArithmeticError: 0^0 is undefined mentioned above. Unfortunately, both implementations now fail due to some other error on legendre_P(4, Zmod(5)(1), 3, 6), but what can you do.

Now the only difference between my implementation and the one in comment 2 is speed: mine is slightly faster unscaled, and that one is slightly faster on a nonstandard interval. But both deltas are negligible, and I have an aesthetic preference for the algorithm that's described in the reference.

This still passes all the tests and is a big improvement over what we have so I'm going to put it up for review again.

### comment:9 Changed 2 years ago by kcrisman

Brief questions:

• When special-casing n=0, you have
```return 1
```
Will this be a Python int?
• What if n is a Python int?
```We can accept Python integers for ``x`` as well::
```
• How's speed? I can imagine the sum thing getting slow.
• Am I misreading this? Should it be beta = 1?
```# From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
```
Oh, I see. It's actually (22.3.1) and maybe this should be in the code with a reference section.

### comment:10 Changed 2 years ago by kcrisman

• Reviewers set to Francis Clarke, Karl-Dieter Crisman

### comment:11 Changed 2 years ago by mjo

When n=0, I try to return unity from the same ring as the variable:

```if n == 0:
# Easy base case, save time. Attempt to return a value in the
# same field/ring as x.
try:
return x.parent()(1)
except AttributeError:
# In case something without a parent was given for x.
return 1
```

The exceptional case is just there as a failsafe, and it would probably give a python int. I guess ZZ(1) is better, so I've changed it.

If n is a python int, we're happy. Early in the function we check that n can be coerced into ZZ. Then later,

```# Ensure that (2**n) is an element of ZZ. This is used later --
# we can divide finite field elements by integers but we can't
# multiply them by rationals at the moment.
n = ZZ(n)
```

For speed, the current implementation has,

```sage: timeit('legendre_P(100,x)')
5 loops, best of 3: 57.1 ms per loop
```

The posted patch has,

```sage: timeit('legendre_P(100,x)')
5 loops, best of 3: 47.3 ms per loop
```

So, even for the unscaled version it's an improvement. There's a slowdown for the scaled version, of course:

```sage: timeit('legendre_P(100,x,-5,10)')
5 loops, best of 3: 70.4 ms per loop
```

But it's still very fast.

You're right about the reference, I must have typo'd it. It's fixed, and I've added a REFERENCES section to the docs.

The new patch is up, I forgot to check the box, so use the one with the "2" at the end.

### comment:12 Changed 2 years ago by kcrisman

• Status changed from needs_review to needs_work

```sage: legendre_P(0,float(1))
1
sage: type(_)
<type 'int'>
```

but now

```sage: legendre_P(0,float(1))
1
sage: type(_)
<type 'sage.rings.integer.Integer'>
```

Work:

• So to my surprise, plotting these things works great. It would be nice to add one or two examples.
• I wouldn't ask for this on this ticket, but it probably should have one or two examples of the errors actually catching silly input. I find it worthwhile to test those branches.

Is that okay? Otherwise I think this is nice.

Last edited 2 years ago by kcrisman (previous) (diff)

### comment:13 Changed 2 years ago by kcrisman

Oh, I should ask the dumb question of whether this scaling is standard in terminology, or whether we need any other disclaimers.

### comment:14 Changed 2 years ago by mjo

• Description modified (diff)
• Status changed from needs_work to needs_review

Ok, I added a plotting example, and three tests for nonsense input. They're at the end of their respective sections.

What might not be standard terminology? "Scaling"? In any case I don't think I'm qualified to say, but it made sense to me at the time.

### comment:15 Changed 2 years ago by mjo

Quick update to use tmp_filename() per Jeroen's comment on #12852.

### comment:16 Changed 2 years ago by jdemeyer

• Milestone changed from sage-5.11 to sage-5.12

### comment:17 Changed 2 years ago by chapoton

apply only sage-trac_13720.patch

### comment:18 follow-up: ↓ 19 Changed 2 years ago by fwclarke

• Status changed from needs_review to needs_work

I really don't understand why a and b get converted into symbolic expressions.  This has some very strange consequences, e.g.,

```sage: legendre_P(3, 7, -1/2, 1/2).parent()
Symbolic Ring
```

It is my understanding that coercion should arrange that the parent of an expression is as close as possible to the parents of the constituant parts.  It is thus wrong to force (almost) everything into the Symbolic Ring.

But anyway, I'm afraid I still prefer my version of the code: (1) for its simplicity (having 3 local functions, two of them only used once, seems far too over-elaborate); (2) my code is significantly faster.

If you don't like that then something like

```    if a == -1 and b == 1:
_init()
return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)), locals={'x':x})
else:
return legendre_P(n, (2*x - a - b)/(b - a), -1, 1)
```

would be a very simple change to the existing code which avoids 'reinventing the wheel'.  This would also have the advantage that an almost identical change would provide scaled versions for the the Legendre functions legendre_Q of the second kind.  However the maxima code is very slow.

Incidentally, in all cases there needs to be a check to see if a == b.

### comment:19 in reply to: ↑ 18 Changed 2 years ago by mjo

I really don't understand why a and b get converted into symbolic expressions.  This has some very strange consequences, e.g.,

```sage: legendre_P(3, 7, -1/2, 1/2).parent()
Symbolic Ring
```

It is my understanding that coercion should arrange that the parent of an expression is as close as possible to the parents of the constituant parts.  It is thus wrong to force (almost) everything into the Symbolic Ring.

I had a comment in the code about this, but I've forgotten the details. In any case -- whatever the issue was -- it seems to be fixed. I removed the manual coercions of n,a,b and all of the tests still pass.

So this is fixed:

```sage: legendre_P(3, 7, -1/2, 1/2).parent()
Rational Field
```

Incidentally, in all cases there needs to be a check to see if a == b.

This will throw a divide-by-zero. What's the alternative?

### comment:20 Changed 2 years ago by mjo

This will throw a divide-by-zero. What's the alternative?

Nevermind, stupid question =)

### comment:21 Changed 2 years ago by mjo

• Branch set to u/mjo/ticket/13720
• Commit set to ce0a4578af9efdef9356d318d698daed1ef27e9b
• Description modified (diff)
• Status changed from needs_work to needs_review

This should address the a==b issue and avoid the coercions.

New commits:

 ​ce0a457 Trac #13720: Allow scaling of legendre_P to an arbitrary interval [a,b].

### comment:22 Changed 2 years ago by mjo

• Authors changed from Michael Orlitzky to Michael Orlitzky, Frédéric Chapoton

Adding an author, as I merged the reviewer patch.

### comment:23 Changed 22 months ago by vbraun_spam

• Milestone changed from sage-6.1 to sage-6.2

### comment:24 Changed 20 months ago by chapoton

• Status changed from needs_review to needs_work

This needs to be rebased.

### comment:25 Changed 20 months ago by git

• Commit changed from ce0a4578af9efdef9356d318d698daed1ef27e9b to 36b4c29920bc4cfbc09b8b7a1fd9de757f8179d3

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

 ​36b4c29 Trac #13720: Allow scaling of legendre_P to an arbitrary interval [a,b].

### comment:26 Changed 20 months ago by mjo

• Status changed from needs_work to needs_review

I had to change the _init() doctest to use the laguerre polynomials instead, since the implementation of legendre_P() no longer uses maxima.

### comment:27 follow-up: ↓ 28 Changed 20 months ago by kcrisman

I'm not sure if they're related, but the Maxima

```assoc_legendre_p[v,u] (z)      Legendre function of degree v and order u
assoc_legendre_q[v,u] (z)      Legendre function, 2nd kind
```

from the symbolic wiki on Trac could be relevant here. If so, we should have conversions - if not, sorry for the noise.

### comment:28 in reply to: ↑ 27 ; follow-up: ↓ 29 Changed 20 months ago by mjo

I'm not sure if they're related, ...

Those two functions are different than our legendre_{P,Q}. Maxima has matching legendre_p(n,x) and legendre_q(n,x) which we're currently using for our implementations.

I updated the note on the wiki page.

### comment:29 in reply to: ↑ 28 Changed 20 months ago by kcrisman

I'm not sure if they're related, ...

Those two functions are different than our legendre_{P,Q}. Maxima has matching legendre_p(n,x) and legendre_q(n,x) which we're currently using for our implementations.

I updated the note on the wiki page.

Thanks, and I updated based on that update, very helpful.

### comment:30 Changed 19 months ago by vbraun_spam

• Milestone changed from sage-6.2 to sage-6.3

### comment:31 Changed 19 months ago by rws

• Status changed from needs_review to needs_work

patchbot:

```sage -t --long src/sage/functions/orthogonal_polys.py  # 33 doctests failed
```

### comment:32 Changed 19 months ago by chapoton

• Branch changed from u/mjo/ticket/13720 to u/chapoton/13720
• Commit changed from 36b4c29920bc4cfbc09b8b7a1fd9de757f8179d3 to 1de3cf47ac8d5d696d2ccf8f1e46281eba4564b4
• Status changed from needs_work to needs_review

I have corrected the problem (and also make a little cleanup). Back to needs review.

New commits:

 ​09ed77c Merge branch 'u/mjo/ticket/13720' of ssh://trac.sagemath.org:22/sage into 13720 ​1de3cf4 trac #13720 pep8 and pyflakes cleanup

### comment:33 Changed 16 months ago by vbraun_spam

• Milestone changed from sage-6.3 to sage-6.4

### comment:34 Changed 12 months ago by git

• Commit changed from 1de3cf47ac8d5d696d2ccf8f1e46281eba4564b4 to 0843342b41d1ce501cab9b84c89e2ef1d84e5dba

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

 ​0843342 Merge branch 'u/chapoton/13720' into 6.5.b2

### comment:36 Changed 9 months ago by rws

Conflicts with #16813. We might want to merge efforts.

### comment:37 Changed 5 months ago by chapoton

• Milestone changed from sage-6.4 to sage-6.8

### comment:38 Changed 7 weeks ago by kdilks

• Branch changed from u/chapoton/13720 to u/kdilks/13720

### comment:39 Changed 7 weeks ago by kdilks

• Branch changed from u/kdilks/13720 to u/chapoton/13720
• Milestone changed from sage-6.8 to sage-6.10

### comment:40 Changed 7 weeks ago by kdilks

• Branch changed from u/chapoton/13720 to u/kdilks/13720
• Commit changed from 0843342b41d1ce501cab9b84c89e2ef1d84e5dba to 8d1354cf63bcca9b01733e4f57aaa16556e3ba0f

New commits:

 ​8d1354c fixed merge conflict

### comment:41 Changed 7 weeks ago by rws

• Status changed from needs_review to needs_work

You introduced merge conflict markers in your branch. Better always eyeball your diff before committing.

### comment:42 Changed 7 weeks ago by kdilks

Yeesh, yeah, thought it was just the one thing, but a lot more needs to be cleaned up. I'll take care of that tomorrow.

### comment:43 Changed 7 weeks ago by kdilks

• Branch changed from u/kdilks/13720 to u/chapoton/13720
• Commit changed from 8d1354cf63bcca9b01733e4f57aaa16556e3ba0f to 0843342b41d1ce501cab9b84c89e2ef1d84e5dba

I'll just change it back and let Frédéric deal with rebasing.

Note: See TracTickets for help on using tickets.