Orthogonal Polynomials and Related Functions for the HP-41
Overview
1°) Legendre Polynomials
2°) Laguerre Polynomials
3°) Generalized Laguerre Polynomials
4°) Hermite Polynomials
5°) Chebyshev Polynomials & the
"Connaissanse des Temps"
a) Chebyshev Polynomials of the
first and second kind
b) The "Connaissance des Temps"
c) Chebyshev Functions
6°) Ultraspherical Polynomials &
Functions
7°) Jacobi's Polynomials & Functions
8°) Fibonacci Polynomials & Functions
9°) Coefficients of Orthogonal Polynomials
-Many of the orthogonal polynomials may be expressed as hypergeometric
functions.
-Thus, they can be generalized to non-integer n-values
-Several programs have been added to compute these functions
-Execution time may of course be saved if you have loaded the M-code
functions "HGF" & "GAM" instead of the corresponding focal programs.
1°) Legendre Polynomials
Formulae: n.Pn(x)
= (2n-1).x.Pn-1(x) - (n-1).Pn-2(x) ;
P0(x) = 1 ; P1(x) = x
Data Registers: /
Flags: /
Subroutines: /
-Synthetic register M may be replaced by any unused data register.
| 01 LBL "LEG"
02 X<>Y 03 E3 04 / 05 1 06 + 07 STO M 08 LASTX 09 LBL 01 10 STO N 11 RCL Z 12 * 13 STO O 14 - 15 RCL M 16 INT 17 1/X 18 1 19 - 20 * 21 RCL O 22 + 23 RCL N 24 X<>Y 25 ISG M 26 GTO 01 27 CLA 28 END |
( 46 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| T | / | x |
| Z | / | x |
| Y | n | Pn-1(x) |
| X | x | Pn(x) |
( 0 < n < 999 )
Example: Calculate P7(4.9)
7 ENTER^
4.9 XEQ "LEG" gives P7(4.9)
= 1698444.019 ( in 4 seconds )
and P6(4.9) = 188641.3852 is
in Y-register.
Note:
-If n is not an integer, Pn(x) may be calculated
by Pn(x) = F(-n,1+n;1,(1-x)/2) where F = the
hypergeometric function
-It is actually a special case of the associated Legendre funcions
Pnm(x) with m = 0
Data Registers: R00 thru R03: temp
Flags: /
Subroutine: "HGF" ( cf "Hypergeometric
Functions for the HP-41" )
| 01 LBL "LGF"
02 X<>Y 03 CHS 04 STO 01 05 1 06 STO 03 07 RCL 01 08 - 09 STO 02 10 1 11 R^ 12 - 13 2 14 / 15 XEQ "HGF" 16 END |
( 25 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | x | Pn(x) |
Example:
0.7 ENTER^
PI 1/X XEQ "LGF" >>>> P0.7(1/PI)
= 0.559708571
---Execution time = 14s---
2°) Laguerre Polynomials
Formulae: Ln(x)
= (2n-1-x).Ln-1(x) - (n-1)2.Ln-2(x)
; L0(x) = 1 ; L1(x) = 1 - x
Data Registers: /
Flags: /
Subroutines: /
-Synthetic register M may be replaced by any unused data register.
| 01 LBL "LAG"
02 X<>Y 03 .1 04 % 05 1 06 + |
07 STO M
08 SIGN 09 LBL 01 10 X<>Y 11 RCL M 12 INT |
13 DSE X
14 "" 15 X^2 16 * 17 RCL M 18 INT |
19 ST+ X
20 DSE X 21 R^ 22 ST- Y 23 X<> T 24 ST* Y |
25 X<> Z
26 - 27 ISG M 28 GTO 01 29 CLA 30 END |
( 51 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| T | / | x |
| Z | / | x |
| Y | n | Ln-1(x) |
| X | x | Ln (x) |
( 0 < n < 1000 )
Example: Calculate L7 (3.14)
7 ENTER^
3.14 XEQ "LAG" >>>> L7 (3.14)
= -4932.43995 (
in 5 seconds )
X<>Y >>>> L6
(3.14) = -189.9784735
Note:
-Some authors divide Ln (x) by n!
In this case simply add RCL M 1 - INT
FACT / after line 28
( but Y Z and T-registers now contain the former Ln-1
(x) )
-With this definition, L7 (3.14) = -0.978658720
3°) Generalized Laguerre Polynomials
Formulae: L0(a)
(x) = 1 ; L1(a) (x) = a+1-x
; n.Ln(a) (x) = (2.n+a-1-x).Ln-1(a)
(x) - (n+a-1).Ln-2(a) (x)
Data Registers: /
Flags: /
Subroutines: /
| 01 LBL "LANX"
02 STO M 03 CLX 04 E3 05 / 06 1 07 ST- Z 08 + |
09 STO N
10 X<>Y 11 STO O 12 CLST 13 SIGN 14 LBL 01 15 X<>Y 16 CHS |
17 RCL N
18 INT 19 RCL O 20 + 21 * 22 RCL N 23 INT 24 ST+ X |
25 RCL O
26 + 27 RCL M 28 - 29 R^ 30 * 31 + 32 RCL N |
33 INT
34 / 35 ISG N 36 GTO 01 37 CLA 38 END |
( 61 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Z | a | Ln-1(a)(x) |
| Y | n | Ln-1(a)(x) |
| X | x | Ln(a)(x) |
( 0 < n < 1000 , n integer )
Example:
1.4 ENTER^
7 ENTER^
PI XEQ "LANX" >>>> L7(1.4)(Pi)
= 1.688893513 ( 5 seconds )
X<>Y L6(1.4)(Pi) = 2.271353727
Notes:
Ln(0)(x) = (1/n!).Ln(x) where Ln(x) = Laguerre's Polynomial
-For non integer n, we can use Ln(a)(x) = [ Gam(n+a+1) / Gam(n+1) / Gam(a+1) ] M(-n,a+1,x)
where Gam = Gamma function
and M
= Kummer's function
Data Registers: R00 thru R02: temp
Flags: /
Subroutines: "GAM" ( cf "Gamma function
for the HP-41" ) , "KUM" ( cf "Kummer's function
for the HP-41" )
| 01 LBL "LANX"
02 RDN 03 CHS 04 STO 01 05 CLX 06 SIGN 07 + 08 STO 02 09 R^ 10 XEQ "KUM" 11 STO 00 12 RCL 02 13 RCL 01 14 - 15 XEQ "GAM" 16 ST* 00 17 1 18 RCL 01 19 - 20 XEQ "GAM" 21 ST/ 00 22 RCL 02 23 XEQ "GAM" 24 ST/ 00 25 RCL 00 26 END |
( 54 bytes / SIZE 003 )
| STACK | INPUTS | OUTPUTS |
| Z | a | / |
| Y | n | / |
| X | x | Ln(a)(x) |
Example:
2 SQRT
7 SQRT
PI 1/X XEQ "LANX" >>>> Lsqrt(7)(sqrt(2))(1/PI)
= 3.625609187
---Execution time = 15s---
4°) Hermite Polynomials
Formulae: Hn(x)
= 2x.Hn-1(x) - 2(n-1).Hn-2(x) ; H0(x)
= 1 ; H1(x) = 2x
Data Registers: /
Flags: /
Subroutines: /
-Synthetic register M may be replaced by any unused data register.
| 01 LBL "HMT"
02 X<>Y 03 .1 04 % 05 1 |
06 +
07 STO M 08 SIGN 09 LBL 01 10 X<>Y |
11 RCL M
12 INT 13 DSE X 14 "" 15 * |
16 CHS
17 X<>Y 18 ST* Z 19 X<> Z 20 + |
21 ST+ X
22 ISG M 23 GTO 01 24 CLA 25 END |
( 42 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| T | / | x |
| Z | / | x |
| Y | n | Hn-1(x) |
| X | x | Hn(x) |
( 0 < n < 1000 )
Example: Calculate H7 (3.14)
7 ENTER^
3.14 XEQ "HMT" >>>> H7 (3.14)
= 73726.24332 (
in 4 seconds )
X<>Y >>>> H6
(3.14) = 21659.28040
Note:
-For fractional n-values, we have Hn(x) = 2n sqrt(PI) [ (1/Gam((1-n)/2)) M(-n/2,1/2,x2) - ( 2.x / Gam(-n/2) ) M((1-n)/2,3/2,x2) ]
where Gam = Gamma function
and
M = Kummer's function
Data Registers: R00 thru R04: temp
Flags: /
Subroutines: "GAM" or "GAM+"
... ( cf "Gamma function for the HP-41" ) , "KUM"
( cf "Kummer's function for the HP-41" )
| 01 LBL "HNX"
02 STO 03 03 X<>Y 04 STO 04 05 .5 06 STO 02 07 * 08 CHS |
09 STO 01
10 RCL 03 11 X^2 12 XEQ "KUM" 13 STO 05 14 RCL 01 15 XEQ "GAM" 16 ST/ 03 |
17 .5
18 ST+ 01 19 ISG 02 20 RCL 01 21 XEQ "GAM" 22 ST/ 05 23 RCL 00 24 XEQ "KUM" |
25 RCL 03
26 ST+ X 27 * 28 RCL 05 29 X<>Y 30 - 31 2 32 RCL 04 |
33 Y^X
34 * 35 PI 36 SQRT 37 * 38 END |
( 69 bytes / SIZE 005 )
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | x | Hn(x) |
Example:
0.4 ENTER^
PI 1/X XEQ "HNX" >>>>
H0.4(1/PI) = 1.010321996
---Execution time = 17s---
Notes:
The results are not very accurate for large x-values.
Hn(x) may also be computed via parabolic cylinder
functions.
5°) Chebyshev Polynomials & the "Connaissance
des Temps"
a) Chebyshev
Polynomials of the first and second kind
Formulae: Tn(x) = 2x.Tn-1(x)
- Tn-2(x) ; T0(x) = 1 ; T1(x)
= x ( first kind ) & Un(x)
= 2x.Un-1(x) - Un-2(x) ; U0(x)
= 1 ; U1(x) = 2x ( second kind
)
Data Registers: /
Flag: F02
if F02 is clear, "CHB" calculates the Chebyshev polynomials of the first
kind: Tn(x)
if F02 is set, "CHB" calculates the Chebyshev polynomials
of the second kind: Un(x)
Subroutines: /
-Synthetic register M may be replaced by any unused data register.
| 01 LBL "CHB"
02 STO M 03 ST+ M 04 FS? 02 05 CLX 06 X<>Y 07 E3 08 / 09 1 10 + 11 STO Z 12 SIGN 13 LBL 01 14 RCL M 15 X<>Y 16 ST* Y 17 X<> Z 18 - 19 ISG Z 20 GTO 01 21 CLA 22 END |
( 40 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Y | n | Tn-1(x)orUn-1(x) |
| X | x | Tn (x)orUn(x) |
( 0 < n < 1000 )
Example: Compute T7 (0.314) and U7(0.314)
CF 02
7 ENTER^
0.314 XEQ "CHB" >>>> T7
(0.314) = -0.786900700 and
T6 (0.314) = 0.338782777 in Y-register
( in 2 seconds )
SF 02
7 ENTER^
0.314 R/S
>>>> U7 (0.314) = -0.582815681
and U6 (0.314) = 0.649952293
in Y-register
b) The "Connaissance
des Temps"
-The "Connaissance des Temps" is a French ephemeris which contains coefficients
of Chebyshev polynomials
for very accurate coordinates of the Sun , the Moon and the
planets.
-If y(t) is a coordinate of a celestial body at a given instant t (
t belonging to the interval [ t0 ; t0 + DT
] ) , y is given by the formula:
y = a0 + a1.T1(x) + ....... + an.Tn(x)
where x = -1 + 2(t-t0)/DT
( -1 <= x <= +1 )
and a0 ; a1
; ....... ; an are published in the Connaissance des Temps
Data Registers: only the (n+1) registers
used to store the coefficients a0 , a1
, ....... , an
Flags: /
Subroutines: /
-Synthetic registers M N O may be replaced by any unused data registers.
| 01 LBL "CdT"
02 X<>Y 03 STO M 04 SIGN 05 ST+ M |
06 STO N
07 RCL Y 08 STO O 09 RCL IND L 10 LBL 01 |
11 R^
12 ST+ X 13 RCL N 14 ST* Y 15 X<> O |
16 -
17 STO N 18 RCL IND M 19 * 20 + |
21 ISG M
22 GTO 01 23 CLA 24 END |
( 46 bytes )
| STACK | INPUTS | OUTPUTS |
| T | / | x |
| Z | / | x |
| Y | bbb.eee | x |
| X | x | y(x) |
where bbb.eee is the control number of the registers which contain the (n+1) coefficients a0 ; a1 ; ....... ; an
Example: Compute the radius vector of Saturn for 2000 March 12 at 0h TT ( TT = TAI+32.184s )
-We find the following coefficients in the "Connaissance des Temps 2000" ( page B37 ) valid for 2000/01/00 0h to 2000/12/33 0h ( DT = 368 days )
a0 = 9.14765315 ; a1 = -0.03544281 ; a2 = 0.00109597 ; a3 = 0.00002140 ; a4 = 0.00000039 ; a5 = -0.00000083 ( in AU )
-For instance, we store these 6 numbers in registers R10 to R15
( bbb.eee = 10.015 )
-In this example t-t0 = 72 days therefore x = -1 +2*72/368
= -0.6086956522
10.015
ENTER^
-0.6086956522 XEQ "CdT" >>>>
radius vector of Saturn = 9.16896253 AU ( in 2
seconds )
c) Chebyshev
Functions
-The Chebyshev function of the 1st kind is defined
by Tn(x) = Cos ( n Acos x )
-But the following routine also uses:
Tn(x) = 2n-1 (x-1)
-n [ 1 + ((x+1)/(x-1))1/2 ] -2n + 2 -n-1
(x-1)n [ 1 + ((x+1)/(x-1))1/2 ] 2n
which is preferable if x > 1
Data Registers: /
Flags: /
Subroutines: /
| 01 LBL "TNX"
02 ABS 03 1 04 X<Y? 05 GTO 00 06 RDN 07 ACOS |
08 *
09 COS 10 RTN 11 LBL 00 12 LASTX 13 STO Z 14 X<>Y |
15 ST+ Z
16 - 17 STO Z 18 / 19 SQRT 20 1 21 + |
22 X^2
23 * 24 2 25 / 26 X<>Y 27 Y^X 28 ENTER^ |
29 1/X
30 + 31 2 32 / 33 END |
( 45 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | x | Tn(x) |
Example:
PI ENTER^
2 XEQ "TNX" >>>> Tpi(2)
= 31.32614116
-Likewise, Tpi(0.4) = -0.877394921
Note:
-If n is not an integer and x < -1 , there will be a DATA ERROR ( the exact result is complex in this case )
-The Chebyshev function of the 2nd kind is defined
by Un(x) = Sin [ (n+1) Acos x ] (1-x2)
-1/2
-But the following routine also uses:
Un(x) = { [ x + ( x2 - 1 )1/2
]n+1 - [ x + ( x2 - 1 )1/2 ] -n-1
} / [ 2 ( x2 - 1 )1/2 ]
if x > 1
Data Register: R00: temp
Flags: /
Subroutines: /
| 01 LBL "UNX"
02 ENTER^ 03 X^2 04 1 05 ST+ T 06 - 07 X>0? |
08 GTO 00
09 RDN 10 ACOS 11 * 12 SIN 13 R^ 14 CHS |
15 SQRT
16 / 17 RTN 18 LBL 00 19 SQRT 20 STO 00 21 RCL Y |
22 SIGN
23 ST* T 24 * 25 + 26 X<>Y 27 Y^X 28 ENTER^ |
29 1/X
30 - 31 RCL 00 32 ST+ X 33 / 34 END |
( 47 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | x | Un(x) |
with x # +/-1
Example:
PI ENTER^
2 XEQ "UNX" >>>> Upi(2)
= 67.48001814
-Likewise Upi(0.4) = -1.086783215
-If n is not an integer, x must be > -1
6°) Ultraspherical Polynomials & Functions
Formulae: C0(a) (x) = 1 ; C1(a) (x) = 2.a.x ; (n+1).Cn+1(a) (x) = 2.(n+a).x.Cn(a) (x) - (n+2a-1).Cn-1(a) (x) if a # 0
Cn(0) (x) = (2/n).Tn(x)
Data Registers: /
Flags:
F02 if a = 0 , none if a # 0
Subroutines: CHB - Chebyshev Polynomials
if a = 0 , none if a # 0
| 01 LBL "USP"
02 1 03 R^ 04 X#0? 05 GTO 00 06 CF 02 07 R^ 08 R^ 09 XEQ "CHB" 10 ST+ X 11 R^ |
12 INT
13 / 14 RTN 15 LBL 00 16 - 17 ST+ X 18 STO O 19 RDN 20 STO M 21 CLX 22 E3 |
23 /
24 1 25 + 26 STO N 27 CLST 28 SIGN 29 LBL 01 30 X<>Y 31 RCL O 32 RCL N 33 INT |
34 -
35 ST* Y 36 X<> L 37 ST+ X 38 RCL O 39 - 40 RCL M 41 * 42 R^ 43 * 44 + |
45 RCL N
46 INT 47 / 48 ISG N 49 GTO 01 50 CLA 51 END |
( 81 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Z | a | / |
| Y | n | Cn-1(a)(x) |
| X | x | Cn(a)(x) |
( 0 < n < 1000 , n integer )
Cn-1(a)(x) is in Y-register only if a # 0
Example:
1.5 ENTER^
7 PI 1/X XEQ
"USP" >>>> C7(1.5)(1/Pi)
= -0.989046386
X<>Y C6(1.5)(1/Pi) =
1.768780932
Notes:
-If n is not an integer, Ultraspherical functions are defined as:
Cn(0)(x) = (2/n) Tn(x)
where Tn(x) = Chebyshev function, first kind
and if a#0:
Cn(a)(x) = [ Gam(n+2a) / Gam(n+1)
/ Gam(2a) ] F( -n , n+2a , n+1/2 , (1-x)/2 ) where
F is the hypergeometric function and Gam = Gamma function
Data Registers: R00 thru R04: temp
Flags: /
Subroutines: "GAM" or "GAM+"
( cf "Gamma function for the HP-41" )
"HGF" ( cf "Hypergeometric functions for the HP-41" )
"TNX" ( cf paragraph 5 ) only if a = 0
| 01 LBL "USP"
02 STO 00 03 1 04 X<>Y 05 - 06 .5 07 STO 03 08 * 09 STO 02 |
10 X<> T
11 X#0? 12 GTO 00 13 X<>Y 14 STO 01 15 RCL 00 16 XEQ "TNX" 17 ST+ X 18 RCL 01 |
19 /
20 RTN 21 LBL 00 22 ST+ 03 23 ST+ X 24 STO 04 25 X<>Y 26 CHS 27 STO 01 |
28 -
29 X<> 02 30 XEQ "HGF" 31 STO 00 32 RCL 02 33 XEQ "GAM" 34 ST* 00 35 1 36 RCL 01 |
37 -
38 XEQ "GAM" 39 ST/ 00 40 RCL 04 41 XEQ "GAM" 42 ST/ 00 43 RCL 00 44 END |
( 82 bytes / SIZE 005 )
| STACK | INPUTS | OUTPUTS |
| Z | a | / |
| Y | n | / |
| X | x | Cn(a)(x) |
Example:
2 SQRT
7 SQRT
PI 1/X XEQ "USP" >>>> Csqrt(7)(sqrt(2))(1/PI)
= -1.575466394 ---Execution
time = 25s---
-Likewise "USP" returns Csqrt(7)(0)(1/PI)
= -0.746600516 in 2 seconds.
7°) Jacobi's Polynomials & Functions
Formulae: P0(a;b) (x) = 1 ; P1(a;b) (x) = (a-b)/2 + x.(a+b+2)/2
2n(n+a+b)(2n+a+b-2).Pn(a;b)
(x) = [ (2n+a+b-1).(a2-b2) + x.(2n+a+b-2)(2n+a+b-1)(2n+a+b)
] Pn-1(a;b) (x)
- 2(n+a-1)(n+b-1)(2n+a+b) Pn-2(a;b) (x)
Data Registers:
R00 thru R05: temp ( R01 = x )
Flags: /
Subroutines: /
| 01 LBL"JCP"
02 STO 01 03 CLX 04 E3 05 / 06 1 07 STO 04 08 + 09 STO 00 10 RDN 11 STO 03 12 X<>Y 13 STO 02 14 + 15 2 16 + 17 RCL 01 18 * |
19 RCL 02
20 RCL 03 21 - 22 + 23 2 24 LBL 01 25 / 26 ISG 00 27 FS? 30 28 GTO 02 29 X<> 04 30 CHS 31 RCL 00 32 INT 33 RCL 02 34 + 35 STO 05 36 1 |
37 -
38 * 39 RCL 00 40 INT 41 RCL 03 42 + 43 ST+ 05 44 1 45 - 46 * 47 RCL 05 48 ENTER^ 49 ST* Z 50 1 51 - 52 STO 05 53 ST* Y 54 LASTX |
55 -
56 * 57 RCL 01 58 * 59 RCL 02 60 X^2 61 RCL 03 62 X^2 63 - 64 RCL 05 65 * 66 + 67 RCL 04 68 * 69 2 70 / 71 + 72 RCL 05 |
73 1
74 - 75 RCL 02 76 RCL 03 77 + 78 RCL 00 79 INT 80 ST+ Y 81 * 82 * 83 GTO 01 84 LBL 02 85 RCL 04 86 X<>Y 87 END |
( 105 bytes / SIZE 006 )
| STACK | INPUTS | OUTPUTS |
| T | a | / |
| Z | b | / |
| Y | n | Pn-1(a;b)(x) |
| X | x | Pn(a;b)(x) |
( 0 < n < 1000 , n integer )
Example:
1.4 ENTER^
1.7 ENTER^
7 ENTER^
PI 1/X XEQ "JCP" >>>> P7(1.4;1.7)(1/Pi)
= -0.322234421 ( in 13 seconds )
X<>Y P6(1.4;1.7)(1/Pi) =
0.538220923
Note:
-The hypergeometric function also allows to evaluate Jacobi's functions even when n is not an integer:
Pn(a;b) (x) = [ Gam(a+n+1)
/ Gam(a+1) / Gam(n+1) ] F ( -n , a+b+n+1 , a+1 , (1-x)/2 )
Data Registers: R00 thru R03: temp
Flag: F00
Subroutines: "GAM" or "GAM+"
... ( cf "Gamma function for the HP-41" ) , "HGF" ( cf
"Hypergeometric functions for the HP-41" )
| 01 LBL "PABNX"
02 STO 00 03 RDN 04 STO 02 05 CHS 06 STO 01 07 X<> Z 08 STO 03 |
09 +
10 1 11 ST+ 03 12 + 13 ST+ 02 14 LASTX 15 RCL 00 16 - |
17 2
18 / 19 XEQ "HGF" 20 STO 02 21 RCL 03 22 RCL 01 23 - 24 XEQ "GAM" |
25 ST* 02
26 1 27 RCL 01 28 - 29 XEQ "GAM" 30 ST/ 02 31 RCL 03 32 FC? 00 |
33 XEQ "GAM"
34 FC?C 00 35 ST/ 02 36 RCL 02 37 END |
( 71 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| T | a | / |
| Z | b | / |
| Y | n | / |
| X | x | Pn(a,b)(x) |
Example:
2 SQRT
3 SQRT
7 SQRT
PI 1/X XEQ "PABNX" >>>> Psqrt(7)(sqrt(2),sqrt(3))(1/PI)
= -0.887945297
---Execution time = 26s---
8°) Fibonacci Polynomials & Functions
-They are defined by F0(x) = 0
F1(x) = 1 and Fn(x)
= x Fn-1(x) + Fn-2(x) if
n > 1
-In fact, the recurrence relation also works for n = 0 and n
= 1 if we define F-1(x) = 1 and F-2(x)
= -x
Data Registers: R00 = x , R01 = counter
Flags: /
Subroutines: /
| 01 LBL "FNX"
02 STO 00 03 CHS 04 1 05 RCL Z 06 E3 07 / 08 STO 01 09 RDN 10 LBL 01 11 STO Z 12 RCL 00 13 * 14 + 15 ISG 01 16 GTO 01 17 END |
( 30 bytes / SIZE 002 )
| STACK | INPUTS | OUTPUTS |
| Y | n | Fn-1(x) |
| X | x | Fn(x) |
Example:
8
ENTER^
PI 1/X XEQ "FNX" >>>> F8(1/PI)
= 1.615692565
---Execution time = 2s---
RDN F7(1/PI) = 1.660297175
Note:
Fn(x) may also be defined if n is not an integer by the formula:
Fn(x) = ( x2 + 4 ) -1/2
[ ( Xn - cos ( n.PI ) X -n ]
where X = [ x + ( x2 + 4 )1/2 ]n
2 -n
Data Register: R00: temp
Flags: /
Subroutines: /
| 01 LBL "FNX"
02 ENTER^ 03 X^2 04 4 05 + 06 SQRT |
07 STO 00
08 RCL Y 09 ABS 10 + 11 2 12 / |
13 X<>Y
14 SIGN 15 R^ 16 * 17 Y^X 18 ENTER^ |
19 1/X
20 1 21 CHS 22 ACOS 23 R^ 24 * |
25 COS
26 * 27 - 28 RCL 00 29 / 30 END |
( 39 bytes / SIZE 001 )
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | x | Fn(x) |
Example:
0.4 ENTER^
PI 1/X XEQ "FNX" >>>> F0.4(1/PI)
= 0.382888178
---Execution time = 2s---
-Of course, this version also works for integer n.
9°) Coefficients of Orthogonal Polynomials
-Most of the orthogonal polynomials may be computed by a recurrence relation of the type:
a(n) Pn(x) = [ b(n) + c(n) x ] Pn-1(x) - d(n) Pn-2(x)
-This program computes the coefficients of the orthogonal
polynomials Pn(x) & Pn-1(x) for a given
n > 0 ( if n = 0 , P = 1 )
-You must specify the 1st register Rbb to store the 1st coefficient
of Pn(x)
Data Registers: R00 = k.nnn
R01 = control number of Pn(x)
R02 = control number of Pn-1(x)
R03 = a R04 = b R05 = c R06 = d ( and R07: temp for "LANX" & "USP" - and R07-R08-R09-R10: temp for "JCP" )
Rbb ..... Ree = Pn(x) Rbb' ..... Ree' = Pn-1(x) with bb' = ee + 1
Flag: F02 ( Chebyshev polynomials
only ) CF 02 = Chebyshev polynomials of the 1st kind
SF 02 = Chebyshev polynomials of the 2nd kind
Subroutines: /
-Lines 270-272-300 are synthetic TEXT0
instructions ( NOPs )
-They may be replaced by STO X ....
| 01 LBL "LEG"
02 XEQ 13 03 LBL 01 04 ISG 00 05 FS? 30 06 GTO 16 07 RCL 00 08 INT 09 ENTER^ 10 STO 03 11 DSE X 12 STO 06 13 + 14 STO 05 15 CLX 16 STO 04 17 XEQ 14 18 GTO 01 19 LBL "LANX" 20 RCL Z 21 STO 07 22 RDN 23 XEQ 13 24 RCL 01 25 RCL IND X 26 CHS 27 STO IND Y 28 RCL 07 29 X<>Y 30 - 31 ISG Y 32 STO IND Y 33 LBL 02 34 ISG 00 35 FS? 30 36 GTO 16 37 RCL 00 38 INT 39 ENTER^ 40 STO 03 41 RCL 07 42 + 43 1 44 CHS 45 STO 05 46 + |
47 STO 06
48 + 49 STO 04 50 XEQ 14 51 GTO 02 52 LBL "CHB" 53 XEQ 13 54 FS? 02 55 ST+ IND 01 56 LBL 03 57 ISG 00 58 FS? 30 59 GTO 16 60 CLX 61 STO 04 62 SIGN 63 STO 03 64 STO 06 65 2 66 STO 05 67 XEQ 14 68 GTO 03 69 LBL "USP" 70 RCL Z 71 STO 07 72 RDN 73 XEQ 13 74 RCL 07 75 ST+ X 76 STO IND 01 77 LBL 04 78 ISG 00 79 FS? 30 80 GTO 16 81 RCL 00 82 INT 83 STO 03 84 STO 06 85 RCL 07 86 1 87 - 88 ST+ 06 89 ST+ 06 90 + 91 ST+ X 92 STO 05 |
93 CLX
94 STO 04 95 XEQ 14 96 GTO 04 97 LBL "HMT" 98 XEQ 13 99 ST+ IND 01 100 LBL 05 101 ISG 00 102 FS? 30 103 GTO 16 104 CLX 105 STO 04 106 SIGN 107 STO 03 108 RCL 00 109 INT 110 ST+ X 111 2 112 STO 05 113 - 114 STO 06 115 XEQ 14 116 GTO 05 117 LBL "JCP" 118 R^ 119 STO 07 120 R^ 121 STO 08 122 R^ 123 R^ 124 XEQ 13 125 RCL 01 126 RCL 07 127 RCL 08 128 + 129 2 130 ST+ Y 131 / 132 STO IND Y 133 RCL 07 134 RCL 08 135 - 136 2 137 / 138 ISG Z |
139 STO IND Z
140 LBL 06 141 ISG 00 142 FS? 30 143 GTO 16 144 RCL 00 145 INT 146 STO 09 147 RCL 07 148 + 149 RCL 08 150 + 151 STO Y 152 RCL 09 153 + 154 STO 10 155 2 156 - 157 * 158 RCL 09 159 ST+ X 160 * 161 STO 03 162 RCL 10 163 1 164 - 165 RCL 07 166 X^2 167 RCL 08 168 X^2 169 - 170 * 171 STO 04 172 RCL 10 173 RCL 10 174 1 175 - 176 ST* Y 177 LASTX 178 - 179 * 180 STO 05 181 RCL 07 182 RCL 08 183 RCL 09 184 1 |
185 -
186 ST+ Z 187 + 188 * 189 RCL 10 190 * 191 ST+ X 192 STO 06 193 XEQ 14 194 GTO 06 195 LBL 13 196 E3 197 / 198 ISG X 199 STO 00 200 X<>Y 201 RCL X 202 1 203 + 204 0 205 STO IND Y 206 CLX 207 E3 208 / 209 + 210 STO 01 211 2.001 212 + 213 STO 02 214 SIGN 215 STO IND 01 216 STO IND 02 217 RTN 218 LBL 14 219 RCL 02 220 FRC 221 E3 222 * 223 RCL 02 224 INT 225 DSE X 226 - 227 E3 228 / 229 RCL 02 230 INT |
231 STO Z
232 + 233 2 234 + 235 E3 236 / 237 + 238 REGMOVE 239 RCL 02 240 1.002 241 + 242 STO 02 243 FRC 244 LASTX 245 INT 246 E3 247 ST* Z 248 / 249 + 250 RCL 01 251 FRC 252 LASTX 253 INT 254 DSE X 255 E3 256 ST* Z 257 / 258 + 259 SIGN 260 LASTX 261 + 262 LBL 08 263 RCL IND Y 264 STO IND Y 265 CLX 266 RCL IND L 267 STO IND Z 268 RDN 269 DSE X 270 "" 271 DSE L 272 "" 273 DSE Y 274 GTO 08 275 RCL IND L 276 STO IND Z |
277 RCL 01
278 E-3 279 + 280 STO 01 281 ISG X 282 0 283 STO IND 01 284 STO IND Y 285 RCL 01 286 RCL 06 287 CHS 288 LBL 09 289 ST* IND Y 290 ISG Y 291 GTO 09 292 RCL 01 293 RCL 02 294 LBL 10 295 RCL 05 296 RCL IND Y 297 * 298 ST+ IND Z 299 ISG Z 300 "" 301 CLX 302 RCL 04 303 LASTX 304 * 305 ST+ IND Z 306 RDN 307 ISG X 308 GTO 10 309 RCL 01 310 RCL 03 311 LBL 11 312 ST/ IND Y 313 ISG Y 314 GTO 11 315 RTN 316 LBL 16 317 RCL 02 318 RCL 01 319 END |
( 491 bytes
/ SIZE ??? )
Legendre Polynomials
| STACK | INPUTS | OUTPUTS |
| Y | bbb | bbb.eee(Ln-1) |
| X | n > 0 | bbb.eee(Ln) |
Example: Find Legendre polynomial of order
n = 6
-If you choose bb = 11
11 ENTER^
6
XEQ "LEG" >>>> 11.017
---Execution time = 40s---
X<>Y 18.023
R11-R12-R13-R14-R15-R16-R17 = [ 231/16 0 -315/16 0 105/16 0 -5/16 ]
-Whence L6(x) = ( 231 x6 - 315 x4 + 105 x2 - 5 ) / 16
R18-R19-R20-R21-R22-R23 = [ 63/8 0 -70/8 0 15/8 0 ]
-Whence L5(x) = ( 63 x5 - 70 x3 + 15 x ) / 8
Generalized Laguerre Polynomials
| STACK | INPUTS | OUTPUTS |
| Z | a | / |
| Y | bbb | bbb.eee(Ln-1(a)) |
| X | n > 0 | bbb.eee(Ln(a)) |
Example: Find Ln(a)(x)
with a = 3 , n = 6
-If you choose bb = 11
3
ENTER^
11 ENTER^
6
XEQ "LANX" >>>> 11.017
---Execution time = 41s---
X<>Y 18.023
R11 ..... R17 give
L6(3)(x) = x6 / 720 - 3 x5
/ 40 + 3 x4 / 2 - 14 x3 + 63 x2 - 126
x + 84
and R18 ... R23
L5(3)(x) = - x5 / 120 + x4
/ 3 - 14 x3 / 3 + 28 x2 - 70 x + 56
Chebyshev Polynomials
| STACK | INPUTS | OUTPUTS |
| Y | bbb | bbb.eee(Pn-1) |
| X | n > 0 | bbb.eee(Pn) |
CF 02 = polynomials of the 1st kind
SF 02 = polynomials of the 2nd kind
Example1: Polynomials of the 1st kind Find T6(x)
-Again with bb = 11
CF 02
11 ENTER^
6 XEQ "CHB" >>>> 11.017
---Execution time = 39s---
X<>Y 18.023
R11 ..... R17 give T6(x) = 32 x6
- 48 x4 + 18 x2 - 1
and R18 ... R23 T5(x)
= 16 x5 - 20 x3 + 5 x
Example2: Polynomials of the 2nd kind Find U6(x)
-Again with bb = 11
SF 02
11 ENTER^
6 XEQ "CHB" >>>> 11.017
---Execution time = 39s---
X<>Y 18.023
R11 ..... R17 give U6(x) = 64 x6
- 80 x4 + 24 x2 - 1
and R18 ... R23 U5(x)
= 32 x5 - 32 x3 + 6 x
Ultraspherical Polynomials
| STACK | INPUTS | OUTPUTS |
| Z | a # 0 | / |
| Y | bbb | bbb.eee(Cn-1(a)) |
| X | n > 0 | bbb.eee(Cn(a)) |
Example: If a = 3 & n = 6 ,
Cn(a)(x) = ?
-If you choose again bb = 11
3
ENTER^
11 ENTER^
6
XEQ "USP" >>>> 11.017
---Execution time = 41s---
X<>Y 18.023
R11 ..... R17 give
C6(3)(x) = 1792 x6 - 1680 x4
+ 360 x2 - 10
and R18 ... R23
C5(3)(x) = 672 x5 - 480 x3
+ 60 x
>>> "USP" does not work if a = 0 but in this case, you can use the formula: Cn(0) (x) = (2/n).Tn(x)
Hermite Polynomials
| STACK | INPUTS | OUTPUTS |
| Y | bbb | bbb.eee(Hn-1) |
| X | n > 0 | bbb.eee(Hn) |
Example: Find Hermite polynomial of order
n = 6
-If you choose bb = 11
11 ENTER^
6
XEQ "HMT" >>>> 11.017
---Execution time = 40s---
X<>Y 18.023
R11 ..... R17 give H6(x) = 64 x6
- 480 x4 + 720 x2 - 120
and R18 ... R23 H5(x)
= 32 x5 - 160 x3 + 120 x
Jacobi Polynomials
| STACK | INPUTS | OUTPUTS |
| T | a | / |
| Z | b | / |
| Y | bbb | Pn-1(a;b)(x) |
| X | n > 0 | Pn(a;b)(x) |
Example: Find Pn(a;b)(x)
with a = 3 , b = 4 , n = 6
-If you choose bb = 11
3
ENTER^
4
ENTER^
11 ENTER^
6
XEQ "JCP" >>>> 11.017
---Execution time = 47s---
X<>Y 18.023
R11 ..... R17 give
P6(3;4)(x) = 423.9375 x6 - 133.875
x5 - 334.6875 x4 + 78.75 x3
+ 59.0625 x2 - 7.875 x - 1.3125
and R18 ... R23
P5(3;4)(x) = 193.375 x5 - 56.875
x4 - 113.75 x3 + 22.75 x2
+ 11.375 x - 0.875
Notes:
-Of course, when the coefficients are not integers, the HP-41 may give
approximate values, not the fractions directly.
-Though mathematically equivalent, evaluating these polynomials would
often produce p(x) with a lower precision because of cancellation of leading
digits,
especially for large n-values: the signs of the coefficients
alternate.
-If you only want to calculate p(x), the programs listed in paragraphs
1°) to 8°) are both faster and more accurate
-They can also be used for larger n-values.
References:
[1] Abramowitz and Stegun , "Handbook of Mathematical Functions"
- Dover Publications - ISBN 0-486-61272-4
[2] http://functions.wolfram.com/