(33S) Legendre Polynomials

[Eddie W. Shore]

Legendre Polynomials

The value of Legendre Polynomials can be calculated using a closed formula from Rodrigues' formula:

P_n(x) = Σ( comb(n, k) * comb(n+k, k) * ((x - 1)/2)^k, k= 0, n)

Program:
HP 33S:
LBL L: Size: LN = 30, CK = 14EC
LBL M: Size: LN = 105, CK = AA05

Run XEQ L.

Code:


L0001  LBL L
L0002  0
L0003  STO P
L0004  INPUT X
L0005  INPUT N
L0006  STO K
M0001  LBL M
M0002  RCL N
M0003  RCL K
M0004  nCr
M0005  RCL N
M0006  RCL+ K
M0007  LASTx
M0008  nCr
M0009  ×
M0010  -1
M0011  RCL+ X
M0012  2
M0013  ÷
M0014  RCL K
M0015  y^x
M0016  ×
M0017  STO+ P
M0018  DSE K
M0019  GTO M
M0020  1
M0021  STO+ P
M0022  RCL P
M0023  RTN

Examples:
N = 2, X = 0.25; Result: -0.406250
N = 3, X = -0.46; Result: 0.446660
N = 4, X = 0.73; Result: -0.380952

Source:
"Legendre polynomials" Wikipedia. https://en.wikipedia.org/wiki/Legendre_polynomials Updated April 6, 2022. Last Accessed April 29, 2022.

[Thomas Klemm]

Or then we use Horner's method to avoid exponentiation, e.g. for \(n = 5\):

\(
\begin{align}
u &= \frac{x - 1}{2} \\
\\
P_5(x) &= 1 + u \cdot \left(30 + u \cdot \left(210 + u \cdot \left(560 + u \cdot \left(630 + u \cdot 252 \right) \right) \right) \right) \\
\end{align}
\)

It's easy to see that:

\(
\begin{align}
\binom{n}{k-1} &= \binom{n}{k} \frac{k}{n-k+1} \\
\\
\binom {n+k-1}{k-1} &= \binom {n+k}{k} \frac{k}{n + k} \\
\end{align}
\)

This leads to a recursive formula to calculate the coefficient from \(k \to k-1\):

\(
\begin{align}
\binom{n}{k-1} \binom {n+k-1}{k-1} = \binom {n}{k}\binom {n+k}{k} \frac{k}{n-k+1} \frac{k}{n + k}
\end{align}
\)

It is used to calculate the next coefficient from the previous one.

Note: All of these coefficients are integers.
Therefore, the order of operation is important.
We also try to keep the intermediate results small.

Python

Code:

from math import comb

def legendre(n, x):
    u = (x - 1) / 2
    a = comb(2*n, n)
    s = 0
    for k in range(n, 0, -1):
        s = a + u * s
        a = a * k // (n + k) * k // (n - k + 1)
    return a + u * s

Example

\( P_5(0.3) \)

5
ENTER
0.3
XEQ "LP"

0.34538625


Programs

In case of the HP-42S and HP-15C (and other models as well) RCL-arithmetic could be used.
Also I tried to avoid stack acrobatics to keep the code similar.
Feel free to optimize the program for your favourite model.

HP-42S

Code:

00 { 52-Byte Prgm }
01▸LBL "LP"
02 1
03 -
04 2
05 ÷
06 STO 00
07 R↓
08 STO 01
09 STO 02
10 RCL 01
11 +
12 RCL 01
13 COMB
14 0
15▸LBL 00
16 RCL 00
17 ×
18 X<>Y
19 +
20 LASTX
21 RCL 02
22 ×
23 RCL 01
24 RCL 02
25 +
26 ÷
27 RCL 02
28 ×
29 RCL 01
30 RCL 02
31 -
32 1
33 +
34 ÷
35 X<>Y
36 DSE 02
37 GTO 00
38 RCL 00
39 ×
40 +
41 END

HP-15C

Code:

   000 {             } 
   001 {    42 21 11 } f LBL A
   002 {           1 } 1
   003 {          30 } −
   004 {           2 } 2
   005 {          10 } ÷
   006 {       44  0 } STO 0
   007 {          33 } R⬇
   008 {       44  1 } STO 1
   009 {       44  2 } STO 2
   010 {       45  1 } RCL 1
   011 {          40 } + 
   012 {       43 36 } g LSTΧ
   013 {       43 40 } g Cy,x
   014 {           0 } 0
   015 {    42 21  0 } f LBL 0
   016 {       45  0 } RCL 0
   017 {          20 } ×
   018 {          34 } x↔y
   019 {          40 } +
   020 {       43 36 } g LSTΧ
   021 {       45  2 } RCL 2
   022 {          20 } ×
   023 {       45  1 } RCL 1
   024 {       45  2 } RCL 2
   025 {          40 } +
   026 {          10 } ÷
   027 {       45  2 } RCL 2
   028 {          20 } ×
   029 {       45  1 } RCL 1
   030 {       45  2 } RCL 2
   031 {          30 } −
   032 {           1 } 1
   033 {          40 } +
   034 {          10 } ÷
   035 {          34 } x↔y
   036 {    42  5  2 } f DSE 2
   037 {       22  0 } GTO 0
   038 {       45  0 } RCL 0
   039 {          20 } ×
   040 {          40 } +
   041 {       43 32 } g RTN

Some models use a specific register I for loops.

HP-11C

Code:

 01-42,21,11     # f LBL A
 02-       1     # 1
 03-      30     # −
 04-       2     # 2
 05-      10     # ÷
 06-    44 0     # STO 0
 07-      33     # R⬇
 08-    44 1     # STO 1
 09-   44 25     # STO I
 10-    45 1     # RCL 1
 11-      40     # + 
 12-   43 36     # g LSTΧ
 13-    43 1     # g Cy,x
 14-       0     # 0
 15-42,21, 0     # f LBL 0
 16-    45 0     # RCL 0
 17-      20     # ×
 18-      34     # x↔y
 19-      40     # +
 20-   43 36     # g LSTΧ
 21-   45 25     # RCL I
 22-      20     # ×
 23-    45 1     # RCL 1
 24-   45 25     # RCL I
 25-      40     # +
 26-      10     # ÷
 27-   45 25     # RCL I
 28-      20     # ×
 29-    45 1     # RCL 1
 30-   45 25     # RCL I
 31-      30     # −
 32-       1     # 1
 33-      40     # +
 34-      10     # ÷
 35-      34     # x↔y
 36-    42 5     # f DSE
 37-    22 0     # GTO 0
 38-    45 0     # RCL 0
 39-      20     # ×
 40-      40     # +
 41-   43 32     # g RTN

For models that don't have a function to compute binomial coefficients, we can use factorial instead:

\(
\begin{align}
\binom{2n}{n} = \frac{(2n)!}{(n!)^2}
\end{align}
\)

HP-67

Code:

001: 31 25 11           # f LBL A 
002: 01                 # 1
003: 51                 # −
004: 02                 # 2
005: 81                 # ÷
006: 33 00              # STO 0
007: 35 53              # h R⬇
008: 33 01              # STO 1
009: 35 33              # h ST I
010: 34 01              # RCL 1
011: 61                 # + 
012: 35 81              # h N!
013: 34 01              # RCL 1
014: 35 81              # h N!
015: 32 54              # g x^2
016: 81                 # ÷
017: 00                 # 0
018: 31 25 00           # f LBL 0
019: 34 00              # RCL 0
020: 71                 # ×
021: 35 52              # h x↔y
022: 61                 # +
023: 35 82              # h LST Χ
024: 35 34              # h RC I
025: 71                 # ×
026: 34 01              # RCL 1
027: 35 34              # h RC I
028: 61                 # +
029: 81                 # ÷
030: 35 34              # h RC I
031: 71                 # ×
032: 34 01              # RCL 1
033: 35 34              # h RC I
034: 51                 # −
035: 01                 # 1
036: 61                 # +
037: 81                 # ÷
038: 35 52              # h x↔y
039: 31 33              # f DSZ
040: 22 00              # GTO 0
041: 34 00              # RCL 0
042: 71                 # ×
043: 61                 # +
044: 35 22              # h RTN

[John Keith]

Not sure this is applicable here, but this section of the Wikipedia article shows a way to compute the coefficients without using COMB. This should be faster since it uses only simple arithmetic.

[Thomas Klemm]

From the aforementioned link:

\(
(n+1)P_{n+1}(x)=(2n+1)xP_{n}(x)-nP_{n-1}(x)
\)

Initial values are:

\(
\begin{align}
P_0(x) &= 1 \\
P_1(x) &= x \\
\end{align}
\)

I knew that we have an expert here: (49g 50g) Gegenbauer and Jacobi Polynomials

From the Wikipedia article Gegenbauer polynomials:

Quote:When α = 1/2, the equation reduces to the Legendre equation, and the Gegenbauer polynomials reduce to the Legendre polynomials.

[Albert Chan]

(05-18-2022 02:59 AM)Eddie W. Shore Wrote:  Legendre Polynomials

The value of Legendre Polynomials can be calculated using a closed formula from Rodrigues' formula:

P_n(x) = Σ( comb(n, k) * comb(n+k, k) * ((x - 1)/2)^k, k= 0, n)

I think P_n as polynomial of (x-1)/2, is meant for symbolic calculations.

BTW, mpmath use this form: hyp2f1(-n, n+1, 1, (1-x)/2).
For integer n, both forms are equivalent (hyp2f1 version can work with arbitrary n)

Code:

function hyp2f1(a,b,c,z)
    local k, s, t = 1, 1, 1
    repeat
        t = t*a*b*z/(c*k)
        s = s + t
        a, b, c, k = a+1, b+1, c+1, k+1
    until s == s-t
    return s
end

function P_bad(n,x) return hyp2f1(-n,n+1,1,(1-x)/2) end

Numerical calculations with big n, coefs blow-up, with alternate signs, we have massive cancellations.
Below example, terms grow as big as ±1E48, before tappering down, making sum garbage.

P₁₀₀(0.25) = 1 - 3787.5 + 3585578.90625 - 1508034728.3203125 + ...

lua> P_bad(100, 0.25)
9.184297382421194e+031

Code:

function P(n,x) -- legendre recurrence relation formula
    local p0, p1 = 1, x
    for k=2,n do k0=k-1; k1=k+k0; p0, p1 = p1, (k1*p1*x-k0*p0)/k end
    return p1
end

Reccurence relation formula is faster, and *much* more accurate.



It may seems we have the same cancellation issues, but half the subtractions are really sums.
The other half, most of subtractions does not suffer catastrophic cancellations.

lua> P(100, 0.25) -- error = -2 ULP
0.07812465474298284

[Thomas Klemm]

Here we go with the recursive formula:

\(
\begin{align}
P_0(x) &= 1 \\
\\
P_1(x) &= x \\
\\
P_{n+1}(x) &= \frac{(2n+1)xP_{n}(x)-nP_{n-1}(x)}{n+1}
\end{align}
\)

HP-25

Code:

01: 23 01    : STO 1
02: 01       : 1
03: 23 02    : STO 2
04: 23 03    : STO 3
05: 22       : Rv
06: 24 02    : RCL 2
07: 24 03    : RCL 3
08: 61       : *
09: 21       : x<->y
10: 23 02    : STO 2
11: 01       : 1
12: 14 73    : f LASTx
13: 51       : +
14: 23 03    : STO 3
15: 14 73    : f LASTx
16: 51       : +
17: 61       : *
18: 24 01    : RCL 1
19: 61       : *
20: 21       : x<->y
21: 41       : -
22: 24 03    : RCL 3
23: 24 00    : RCL 0
24: 14 71    : f x=y
25: 13 29    : GTO 29
26: 22       : Rv
27: 71       : /
28: 13 06    : GTO 06
29: 22       : Rv
30: 71       : /
31: 13 00    : GTO 00

Registers

0: \(n\) … the order of the Legendre polynomial
1: \(x\) … where to evaluate the polynomial at
2: \(P_{k-1}\) … the previous value
3: \(k\) … the counter


Examples

I've changed the usage slightly since entering the order \( n \) again and again annoyed me.

Initialisation

5
STO 0
CLEAR PRGM

Calculation

0.3
R/S

0.345386250

0.6
R/S

-0.152640000


This leads to a blinkenlicht festival.

100
STO 0
0.25
R/S

0.07812465477

[Thomas Klemm]

Or then we use program L4 (Legendre Polynomials) of the HP-25 Library of the PPC Journal:

Code:

01: 23 01    : STO 1
02: 01       : 1
03: 23 02    : STO 2
04: 23 04    : STO 4
05: 41       : -
06: 24 01    : RCL 1
07: 15 51    : g x>=0
08: 13 21    : GTO 21
09: 01       : 1
10: 51       : +
11: 71       : /
12: 15 03    : g ABS
13: 14 02    : f SQRT
14: 14 07    : f LN
15: 23 02    : STO 2
16: 24 01    : RCL 1
17: 61       : *
18: 01       : 1
19: 51       : +
20: 32       : CHS
21: 23 03    : STO 3
22: 24 00    : RCL 0
23: 24 04    : RCL 4
24: 14 51    : f x>=y
25: 13 44    : GTO 44
26: 24 02    : RCL 2
27: 61       : *
28: 24 04    : RCL 4
29: 02       : 2
30: 61       : *
31: 01       : 1
32: 23 51 04 : STO + 4
33: 51       : +
34: 24 03    : RCL 3
35: 23 02    : STO 2
36: 24 01    : RCL 1
37: 15 03    : g ABS
38: 61       : *
39: 61       : *
40: 41       : -
41: 24 04    : RCL 4
42: 71       : /
43: 13 20    : GTO 20
44: 14 71    : f x=y
45: 13 48    : GTO 48
46: 24 02    : RCL 2
47: 13 00    : GTO 00
48: 24 03    : RCL 3
49: 13 00    : GTO 00

Examples

CLEAR RPGM

4
STO 0
R/S

1060.375

9
STO 0
0.5
R/S

-0.267898559


The program also implements the Legendre functions of the second kind (Qn)

Use a negative value for \( x \) in this case:

0
STO 0
-0.5
R/S

0.549306145

10
STO 0
-0.4
R/S

0.373991229


Resources

(05-17-2018 02:31 PM)rprosperi Wrote:  Here's a link: HP-25 Library Programs-PPC Jrnl

(05-18-2018 04:07 AM)teenix Wrote:  I have extracted about 63 of the HP-25 programs and have them in a zip file at

http://www.teenix.org/HP25ppc.zip

They are in a format that will load into the HP emulator and MultiCalc and they are also in text format.

[Albert Chan]

From https://solitaryroad.com/c679.html, P(n,x) and Q(n,x) are based on same recurrence relation formula.
Note: this version extended range to all non-negative integers n. (i.e. n=0 included)

Code:

function legendre_hlp(n,x,p0,k0)
    local p1,k1 = x*p0-k0
    for k=2,n do k0=k-1; k1=k+k0; p0, p1 = p1, (k1*p1*x-k0*p0)/k end
    return (n==0 and p0) or p1
end
    
function P(n,x) return legendre_hlp(n,x,1,0) end
function Q(n,x) return legendre_hlp(n,x,atanh(x),1) end

lua> P(4,4)
1060.375
lua> P(9,0.5)
-0.2678985595703125

lua> Q(0, 0.5)
0.5493061443340549
lua> Q(10, 0.4)
0.37399122844670774

Numbers matched L4 (Legendre Polynomials) numbers (previous post)

(05-20-2022 04:10 PM)Thomas Klemm Wrote:  The program also implements the Legendre functions of the second kind (Qn)

Use a negative value for \( x \) in this case ...

What happens if x = 0 ? Which legendre kind get returned ?
I would assume first kind (x = +0), since signed zero may not be supported.

This may be confusing, since P(n,0) ≠ Q(n,0)

[Thomas Klemm]

(05-20-2022 05:51 PM)Albert Chan Wrote:  What happens if x = 0 ? Which legendre kind get returned ?

If only there was an online emulator for the HP-25, with instructions on how to load a program:

Quote:As with my other HTML5 microcode emulators, tap the display (LED digits) to bring up more options.

(05-20-2022 05:51 PM)Albert Chan Wrote:  I would assume first kind (x = +0), since signed zero may not be supported.

If only polynomials were continuous and we could enter tiny numbers like \( -10^{-99} \) or \( 10^{-99} \), that would be great.

And if we only could understand the meaning of these instructions:

Code:

06: 24 01    : RCL 1
07: 15 51    : g x>=0
08: 13 21    : GTO 21

I highly recommend to consult the original contribution by James J. Davidson (547) in his neat handwriting.
It contains stack diagrams and comments that are helpful to understand the program.