+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (33S) Legendre Polynomials (/thread-18385.html)
(33S) Legendre Polynomials - Eddie W. Shore - 05-18-2022 02:59 AM
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 RTNExamples:
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.
RE: (33S) Legendre Polynomials - Thomas Klemm - 05-18-2022 04:02 PM
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 * sExample
\( 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 ENDHP-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 RTNSome 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 RTNFor 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 RTNRE: (33S) Legendre Polynomials - John Keith - 05-18-2022 05:11 PM
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.
RE: (33S) Legendre Polynomials - Thomas Klemm - 05-18-2022 07:34 PM
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.
RE: (33S) Legendre Polynomials - Albert Chan - 05-18-2022 08:34 PM
(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) endNumerical 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
endReccurence relation formula is faster, and *much* more accurate.
![[Image: 7b08cb339066f2105ff011f0f47d2630a481c27f]](https://wikimedia.org/api/rest_v1/media/math/render/svg/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
RE: (33S) Legendre Polynomials - Thomas Klemm - 05-18-2022 11:13 PM
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 00Registers
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
RE: (33S) Legendre Polynomials - Thomas Klemm - 05-20-2022 04:10 PM
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 00Examples
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.
RE: (33S) Legendre Polynomials - Albert Chan - 05-20-2022 05:51 PM
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) endlua> 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)
RE: (33S) Legendre Polynomials - Thomas Klemm - 05-20-2022 08:08 PM
(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 21I 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.