(33S) Legendre Polynomials

[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