Bionic Orthopoly

# 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...