Bionic Orthogonal Polynomials for the HP-41
Overview
0°) Extra Registers & M-Code Routines
1°) Legendre Polynomials
2°) Generalized Laguerre's Polynomials
3°) Hermite Polynomials
4°) Chebyshev Polynomials
5°) UltraSpherical Polynomials
6°) Jacobi Polynomials
7°) Associated Legrendre Functions
of the 1st kind - Integer Order and Degree
0°) Extra Registers & M-Code Routines
-The following programs use 1 to 3 extra-registers that I've called:
U V W
-STO W & RCL W are listed in "A few M-Code Routines
for the HP-41"
-Registers U and V are coded in the same way.
-Of course, if n is not too large, you can replace these registers by R97 R98 R99
-Other M-Code routines are employed:
X/E3 X+1 X-1 X/2 cf "A few M-Code Routines for the HP-41"
-The product of 2 bi-ons with proportional "imaginary parts" B*B1 cf "Bi-ons for the HP-41"
-And M-Code routines that are useful for anions may also be used for bi-ons without any modification:
A+A A-A ST*A ST/A AMOVE ASWAP
cf "Anions" and "Anionic Functions for the HP-41"
1°) Legendre Polynomials
Formulae:
m.Pm(b) = (2m-1).b.Pm-1(b)
- (m-1).Pm-2(b) , P0(b) = 1 ,
P1(b) = b
Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "LEGB" )
• R01 ...... • R2n = the n components of the bion b which are saved in registers R2n+1 to R4n
R4n+1 ........ R6n = Pm(b)
R6n+1 ........ R8n = Pm-1(b)
>>> When the program stops: R01 ...... R2n = the n components of the result
Flags: /
Subroutines: X/E3 X-1
AMOVE ASWAP ST*A ST/A A-A B*B1 STO
W RCL W
| 01 LBL "LEGB"
02 X/E3 03 STO W 04 12 05 AMOVE 06 CLX 07 ST*A 08 14 09 AMOVE 10 SIGN 11 STO 01 |
12 LBL 01
13 13 14 AMOVE 15 RCL W 16 ENTER^ 17 ISG X 18 X=0? 19 GTO 00 20 STO W 21 41 22 AMOVE |
23 X<> Z
24 INT 25 ST*A 26 14 27 AMOVE 28 31 29 AMOVE 30 B*B1 31 RCL W 32 INT 33 ST+ X |
34 X-1
35 ST*A 36 34 37 ASWAP 38 23 39 ASWAP 40 A-A 41 RCL W 42 INT 43 ST/A 44 32 |
45 AMOVE
46 GTO 01 47 LBL 00 48 RCL 00 49 X/E3 50 ISG X 51 END |
( 98 bytes
/ SIZE 8n+1 )
| STACK | INPUT | OUTPUT |
| X | m | 1.2n |
where m is an integer
( 0 <= m < 1000 )
and 1.2n
is the control number of the result
Example: m = 7 b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3 ( biquaternion -> 8 STO 00 )
0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08
7 XEQ "LEGB" >>>> 1.008 ---Execution time = 37s---
And we get in registers R01 thru R08
P7(b) = ( -240.8232644
+ 180.2220151 i ) + ( -8.4974508 + 82.1615056
i ) e1
+ ( -6.6831028 + 128.8517448 i ) e2 + ( -4.8687548
+ 175.541984 i ) e3
2°) Generalized Laguerre's Polynomials
Formulae:
L0(a) (b) = 1
, L1(a) (b) = a+1-b ,
m Lm(a) (b) = (2.m+a-1-b) Lm-1(a)
(b) - (m+a-1) Lm-2(a) (b)
Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "LANB" )
• R01 ...... • R2n = the n components of the bion b which are saved in registers R2n+1 to R4n
R4n+1 ........ R6n = Lm(b)
R6n+1 ........ R8n = Lm-1(b)
>>> When the program stops: R01 ...... R2n = the n components of the result
Flags: /
Subroutines: X/E3 X-1
AMOVE ASWAP ST*A ST/A A-A B*B1
STO V RCL V STO W RCL W
| 01 LBL "LANB"
02 X/E3 03 STO W 04 X<>Y 05 STO V 06 12 07 AMOVE 08 CLX 09 ST*A 10 14 11 AMOVE 12 SIGN |
13 STO 01
14 LBL 01 15 13 16 AMOVE 17 RCL W 18 ISG X 19 X=0? 20 GTO 00 21 STO W 22 23 23 AMOVE 24 12 |
25 ASWAP
26 X<> Z 27 INT 28 ST+ X 29 X-1 30 RCL V 31 + 32 ST- 01 33 B*B1 34 14 35 ASWAP 36 SIGN |
37 RCL W
38 INT 39 - 40 RCL V 41 - 42 ST*A 43 24 44 ASWAP 45 A-A 46 RCL W 47 INT 48 ST/A |
49 32
50 AMOVE 51 GTO 01 52 LBL 00 53 RCL 00 54 X/E3 55 ISG X 56 END |
( 104 bytes / SIZE 8n+1 )
| STACK | INPUT | OUTPUT |
| Y | a | / |
| X | m | 1.2n |
where m is an integer
( 0 <= m < 1000 ) , a is a real number
and
1.2n is the control number of the result
Example: a = sqrt(2) m = 7 b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3 ( biquaternion -> 8 STO 00 )
0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08
2 SQRT
7 XEQ "LANB"
>>>> 1.008
---Execution time = 39s---
And we get in registers R01 thru R08
L7sqrt(2) (b) =
( -5.035303539 - 80.93462041 i ) + ( -27.58137526
- 1.567289524 i ) e1
+ ( -43.15232859 - 0.238461640 i ) e2 + ( -58.72328181
+ 1.090366269 i ) e3
3°) Hermite Polynomials
Formulae:
H0(b) = 1
, H1(b) = 2 b
Hm(b) = 2.b Hm-1(b)
- 2 (m-1) Hm-2(b)
Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "HMTB" )
• R01 ...... • R2n = the n components of the bion b which are saved in registers R2n+1 to R4n
R4n+1 ........ R6n = Hm(b)
R6n+1 ........ R8n = Hm-1(b)
>>> When the program stops: R01 ...... R2n = the n components of the result
Flags: /
Subroutines: X/E3 AMOVE
ASWAP ST*A ST/A A-A B*B1
STO W RCL W
| 01 LBL "HMTB"
02 X/E3 03 STO W 04 12 05 AMOVE 06 CLX 07 ST*A 08 14 09 AMOVE |
10 SIGN
11 STO 01 12 LBL 01 13 13 14 AMOVE 15 RCL W 16 ENTER^ 17 ISG X 18 X=0? |
19 GTO 00
20 STO W 21 41 22 AMOVE 23 X<> Z 24 INT 25 ST*A 26 14 27 AMOVE |
28 31
29 AMOVE 30 B*B1 31 34 32 ASWAP 33 23 34 ASWAP 35 A-A 36 2 |
37 ST*A
38 32 39 AMOVE 40 GTO 01 41 LBL 00 42 RCL 00 43 X/E3 44 ISG X 45 END |
( 87 bytes / SIZE 8n+1 )
| STACK | INPUT | OUTPUT |
| X | m | 1.2n |
where m is an integer
( 0 <= m < 1000 )
and 1.2n
is the control number of the result
Example: m = 7 b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3 ( biquaternion -> 8 STO 00 )
0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08
7 XEQ "HMTB" >>>> 1.008 ---Execution time = 34s---
And we get in registers R01 thru R08
H7(b) = ( 3548.365108
+ 3989.424742 i ) + ( 2872.096870 - 177.1388928
i ) e1
+ ( 4466.300006 - 506.1044224 i ) e2 + ( 6060.503142
- 835.0699520 i ) e3
4°) Chebyshev Polynomials
Formulae:
CF 02
Tm(b) = 2b.Tm-1(b) - Tm-2(b)
; T0(b) = 1 ; T1(b) = b
( first kind )
SF 02
Um(b) = 2b.Um-1(b) - Um-2(b) ;
U0(b) = 1 ; U1(b) = 2b
( second kind )
Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "CHBB" )
• R01 ...... • R2n = the n components of the bion b
R4n+1 ........ R6n = Tm(b) or Um(b)
R6n+1 ........ R8n = Tm-1(b) or Um-1(b)
>>> When the program stops: R01 ...... R2n = the n components of the result
Flag: F02
if F02 is clear, "CHBB" calculates the Chebyshev polynomials of the first
kind: Tm(b)
if F02 is set, "CHBB" calculates the Chebyshev polynomials
of the second kind: Um(b)
Subroutines: X/E3 X-1 AMOVE
ASWAP ST*A A-A B*B1
STO W RCL W
| 01 LBL "CHBB"
02 X/E3 03 STO W 04 12 05 AMOVE 06 14 07 AMOVE |
08 2
09 ST*A 10 12 11 ASWAP 12 CLX 13 ST*A 14 14 |
15 FS? 02
16 AMOVE 17 SIGN 18 STO 01 19 LBL 01 20 13 21 AMOVE |
22 RCL W
23 ISG X 24 X=0? 25 GTO 00 26 STO W 27 B*B1 28 34 |
29 ASWAP
30 23 31 ASWAP 32 A-A 33 23 34 ASWAP 35 GTO 01 |
36 LBL 00
37 INT 38 X-1 39 RCL 00 40 X/E3 41 ISG X 42 END |
( 82 bytes / SIZE 8n+1 )
| STACK | INPUT | OUTPUT |
| Y | / | m |
| X | m | 1.2n |
where m is an integer
( 0 <= m < 1000 )
and 1.2n
is the control number of the result
Example: m = 7 b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3 ( biquaternion -> 8 STO 00 )
0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08
• Chebyshev Polynomial 1st kind
CF 02
7 XEQ "CHBB"
>>>> 1.008
---Execution time = 25s---
And we get in registers R01 thru R08
T7(b) = ( -558.0438464
+ 462.2275712 i ) + ( -7.670764800 + 203.3249536
i ) e1
+ ( 4.299603200 + 317.8005888 i ) e2 + ( 16.26997120
+ 432.2762240 i ) e3
• Chebyshev Polynomial 2nd kind
SF 02
7 XEQ "CHBB"
>>>> 1.008
---Execution time = 25s---
And we get in registers R01 thru R08
U7(b) = ( -1176.666563
+ 804.7309824 i ) + ( -61.19336960 + 378.9647872
i ) e1
+ ( -65.14447360 + 596.0805376 i ) e2 + ( -69.09557760
+ 813.1962880 i ) e3
5°) UltraSpherical Polynomials
Formulae:
C0(a) (b) = 1 ; C1(a) (b) = 2.a.b ; (m+1).Cm+1(a) (b) = 2.(m+a).b.Cm(a) (b) - (m+2a-1).Cm-1(a) (b) if a # 0
Cm(0) (b)
= (2/m).Tm(b)
Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "USPB" )
• R01 ...... • R2n = the n components of the bion b which are saved in registers R2n+1 to R4n
R4n+1 ........ R6n = Cm(a) (b)
R6n+1 ........ R8n = Cm-1(a) (b)
>>> When the program stops: R01 ...... R2n = the n components of the result
Flags: /
Subroutines: X/E3 X-1
AMOVE ASWAP ST*A ST/A A+A B*B1
STO W RCL W STO V RCL V
"CHBB" if a = 0 ( cf paragraph above )
| 01 LBL "USPB"
02 X<>Y 03 X#0? 04 GTO 00 05 X<>Y 06 CF 02 07 XEQ "CHBB" 08 2 09 RCL Z 10 / 11 ST*A 12 GTO 02 13 LBL 00 14 STO V |
15 X<>Y
16 X/E3 17 STO W 18 12 19 AMOVE 20 CLX 21 ST*A 22 14 23 AMOVE 24 SIGN 25 STO 01 26 LBL 01 27 13 28 AMOVE |
29 RCL W
30 ISG X 31 X=0? 32 GTO 02 33 STO W 34 B*B1 35 RCL W 36 INT 37 X-1 38 RCL V 39 + 40 ST+ X 41 ST*A 42 14 |
43 ASWAP
44 2 45 RCL W 46 INT 47 - 48 RCL V 49 ST+ X 50 - 51 ST*A 52 34 53 ASWAP 54 23 55 ASWAP 56 A+A |
57 RCL W
58 INT 59 ST/A 60 23 61 ASWAP 62 GTO 01 63 LBL 02 64 RCL 00 65 X/E3 66 ISG X 67 END |
(
124 bytes / SIZE 8n+1 )
| STACK | INPUT | OUTPUT |
| Y | a | / |
| X | m | 1.2n |
where m is an integer
( 0 <= m < 1000 ) , a is a real number
and
1.2n is the control number of the result
Example: a = sqrt(2) m = 7 b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3 ( biquaternion -> 8 STO 00 )
0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08
2 SQRT
7 XEQ "USPB"
>>>> 1.008
---Execution time = 40s---
And we get in registers R01 thru R08
C7sqrt(2) (b) =
( -3142.090579 + 2008.230127 i ) + ( -198.7572026
+ 969.5883167 i ) e1
+ ( -232.4941704 + 1528.458350 i ) e2 + ( -266.2311383
+ 2087.328384 i ) e3
-Likewise, you will find:
C7(0) (b) = ( -159.4410990
+ 132.0650203 i ) + ( -2.191647086 + 58.09284388
i ) e1
+ ( 1.228458057 + 90.80016822 i ) e2 + ( 4.648563200
+ 123.5074926 i ) e3
6°) Jacobi Polynomials
Formulae:
P0(a;c) (b) = 1 ; P1(a;c) (b) = (a-c)/2 + b (a+c+2)/2
2m(m+a+c)(2m+a+c-2) Pm(a;c)
(b) = [ (2m+a+c-1).(a2-c2) + b (2m+a+c-2)(2m+a+c-1)(2m+a+c)
] Pm-1(a;c) (b)
- 2(m+a-1)(m+c-1)(2m+a+c) Pm-2(a;c) (b)
Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "JCPB" )
• R01 ...... • R2n = the n components of the bion b which are saved in registers R2n+1 to R4n
R4n+1 ........ R6n = Pm(a;c) (b)
R6n+1 ........ R8n = Pm-1(a;c) (b)
>>> When the program stops: R01 ...... R2n = the n components of the result
Flags: /
Subroutines: X/E3 X-1
AMOVE ASWAP ST*A ST/A A+A B*B1
STO W RCL W STO V RCL V STO U RCL U
-Lines 22 & 96 are three-byte GTO's
| 01 LBL "JCPB"
02 X/E3 03 STO W 04 RDN 05 STO V 06 X<>Y 07 STO U 08 12 09 AMOVE 10 CLX 11 ST*A 12 14 13 AMOVE 14 SIGN 15 STO 01 16 LBL 01 17 13 |
18 AMOVE
19 RCL W 20 ISG X 21 X=0? 22 GTO 00 23 STO W 24 23 25 AMOVE 26 12 27 ASWAP 28 X<> Z 29 INT 30 ST+ X 31 RCL U 32 + 33 RCL V 34 + |
35 ENTER^
36 STO M 37 X-1 38 ENTER^ 39 STO N 40 X-1 41 * 42 * 43 ST*A 44 RCL N 45 RCL U 46 X^2 47 RCL V 48 X^2 49 - 50 * 51 ST+ 01 |
52 B*B1
53 14 54 ASWAP 55 SIGN 56 RCL W 57 INT 58 STO M 59 - 60 RCL U 61 - 62 STO N 63 ST+ X 64 RCL M 65 RCL V 66 + 67 X-1 68 STO O |
69 *
70 RCL M 71 ST+ X 72 RCL U 73 + 74 RCL V 75 + 76 * 77 ST*A 78 24 79 ASWAP 80 A+A 81 RCL O 82 RCL N 83 - 84 RCL M 85 RCL U |
86 +
87 RCL V 88 + 89 * 90 RCL M 91 ST+ X 92 * 93 ST/A 94 32 95 AMOVE 96 GTO 01 97 LBL 00 98 RCL 00 99 X/E3 100 ISG X 101 CLA 102 END |
( 178 bytes / SIZE 8n+1 )
| STACK | INPUT | OUTPUT |
| Z | a | / |
| Y | c | / |
| X | m | 1.2n |
where m is an integer
( 0 <= m < 1000 ) , a & c are real numbers
and
1.2n is the control number of the result
Example: a = sqrt(2) c = sqrt(3) m = 7 b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3 ( biquaternion -> 8 STO 00 )
0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08
2 SQRT
3 SQRT
7 XEQ "JCPB"
>>>> 1.008
---Execution time = 55s---
And we get in registers R01 thru R08
P7sqrt(2) , sqrt(3) (b)
= ( -1604.402811 + 779.134567 i ) + ( -142.7952083
+ 499.3596976 i ) e1
+ ( -182.8117487 + 790.4247449 i ) e2 + ( -222.8282902
+ 1081.489792 i ) e3
7°) Associated Legrendre Functions of the 1st
kind - Integer Order and Degree
"PMNB" compute the functions of type 2 - if CF 03 - or the functions
of type 3 - if SF 03.
Formulae: (n-m) Pnm(a) = a (2n-1) Pn-1m(a) - (n+m-1) Pn-2m(a)
Type 2
Pmm(a) = (-1)m
(2m-1)!! ( 1-a2 )m/2
where (2m-1)!! = (2m-1)(2m-3)(2m-5).......5.3.1
Type 3
Pmm(a) = (2m-1)!!
( a2-1 )m/2
Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "PMNB" )
• R01 ...... • R2n = the n components of the bion b which are saved in registers R2n+1 to R4n
R4n+1 ........ R6n = Pnm (b)
R6n+1 ........ R8n = Pn-1m (b)
>>> When the program stops: R01 ...... R2n = the n components of the result
Flags: /
Subroutines: X/E3 X-1
X/2 AMOVE ASWAP ST*A ST/A A+A B*B1
"B^X"
STO W RCL W STO V RCL V STO U RCL U
| 01 LBL "PMNB"
02 STO V 03 X<>Y 04 STO U 05 12 06 AMOVE 07 B*B1 08 SIGN 09 ST- 01 10 CHS 11 FC? 03 12 ST*A 13 RCL U 14 X/2 15 XEQ "B^X" 16 RCL U |
17 ST+ X
18 1 19 ST+ Y 20 X<>Y 21 LBL 01 22 2 23 - 24 X>0? 25 ST* Y 26 X>0? 27 GTO 01 28 CLX 29 SIGN 30 FC? 03 31 CHS 32 RCL U |
33 Y^X
34 * 35 ST*A 36 RCL 00 37 4 38 * 39 .1 40 % 41 ISG X 42 + 43 RCL 00 44 - 45 CLRGX 46 RCL V 47 X/E3 48 RCL U |
49 +
50 STO W 51 LBL 02 52 13 53 AMOVE 54 RCL W 55 ISG X 56 X=0? 57 GTO 00 58 STO W 59 INT 60 ST+ X 61 X-1 62 ST*A 63 B*B1 64 4 |
65 RCL 00
66 ST* Y 67 ENTER^ 68 SIGN 69 RCL U 70 - 71 RCL W 72 INT 73 - 74 SIGN 75 LBL 03 76 CLX 77 X<> IND Z 78 LASTX 79 * 80 ST+ IND Y |
81 DSE Z
82 DSE Y 83 GTO 03 84 RCL W 85 INT 86 RCL U 87 - 88 ST/A 89 34 90 AMOVE 91 GTO 02 92 LBL 00 93 RCL 00 94 X/E3 95 ISG X 96 END |
( 160 bytes / SIZE 8n+1 )
| STACK | INPUTS | OUTPUTS |
| Y | m | / |
| X | n | 1.2n |
where m
& n are integers with 0 <=
m
<= n
and 1.2n
is the control number of the result.
Example: m = 3 , n = 7 , b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3 ( biquaternion -> 8 STO 00 )
0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08
• CF 03 Legendre Function of type 2
3 ENTER^
7 XEQ "PMNB" >>>> 1.008
---Execution time = 48s---
And we have in R01 thru R08
P37 ( b ) = ( 6071.566263
+ 53115.31993 i ) + ( 17513.35572 - 15006.54081
i ) e1
+ ( 26120.31168 - 24811.27212 i ) e2 + ( 34727.26763
- 34016.00343 i ) e3
• SF 03 Legendre Function of type 3
-Store again 0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08
3 ENTER^
7 XEQ "PMNB" >>>> 1.008
---Execution time = 48s---
And we have in R01 thru R08
P37 ( b ) = ( -46823.11240
+ 45127.54828 i ) + ( 1048.341587 + 18931.32096
i ) e1
+ ( -3149.918553 + 29448.99338 i ) e2 + ( 5251.495520
+ 39966.66578 i ) e3
Notes:
-This program does not check if m & n
are integers that satisfy 0 <= m <= n
-Register W may be replaced by register V.
-These functions are not always polynomials, but the method is similar
to the other routines of this page...