OrthoPoly

# 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  XY 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/