# Ticket #13720(needs_review enhancement)

Opened 6 months ago

## Scale legendre_P to [a,b]

Reported by: Owned by: mjo burcin major sage-5.10 symbolics N/A Michael Orlitzky

### Description

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

## Change History

### comment:1 Changed 6 months ago by mjo

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

### comment:2 follow-up: ↓ 3 Changed 6 months 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 6 months 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.

### Changed 6 months ago by mjo

Fix variable quoting in exceptions and allow symbolic endpoints

### comment:4 follow-up: ↓ 5 Changed 6 months 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 6 months 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 6 months 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 4 weeks 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 4 weeks ago by mjo

• Status changed from needs_work to needs_review