hp41programs

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