Post Reply 
(33S) Legendre Polynomials
05-18-2022, 02:59 AM
Post: #1
(33S) Legendre Polynomials
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.
Visit this user's website Find all posts by this user
Quote this message in a reply
05-18-2022, 04:02 PM
Post: #2
RE: (33S) Legendre Polynomials
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
Find all posts by this user
Quote this message in a reply
05-18-2022, 05:11 PM
Post: #3
RE: (33S) Legendre Polynomials
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.
Find all posts by this user
Quote this message in a reply
05-18-2022, 07:34 PM
Post: #4
RE: (33S) Legendre Polynomials
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.
Find all posts by this user
Quote this message in a reply
05-18-2022, 08:34 PM (This post was last modified: 05-18-2022 09:08 PM by Albert Chan.)
Post: #5
RE: (33S) Legendre Polynomials
(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 Pn 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.

P100(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.

[Image: 7b08cb339066f2105ff011f0f47d2630a481c27f]

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
Find all posts by this user
Quote this message in a reply
05-18-2022, 11:13 PM
Post: #6
RE: (33S) Legendre Polynomials
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
Find all posts by this user
Quote this message in a reply
05-20-2022, 04:10 PM
Post: #7
RE: (33S) Legendre Polynomials
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.
Find all posts by this user
Quote this message in a reply
05-20-2022, 05:51 PM (This post was last modified: 05-20-2022 06:39 PM by Albert Chan.)
Post: #8
RE: (33S) Legendre Polynomials
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)
Find all posts by this user
Quote this message in a reply
05-20-2022, 08:08 PM
Post: #9
RE: (33S) Legendre Polynomials
(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.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)